numerical integration inf to inf
4 views (last 30 days)
Show older comments
PTK on 30 Jan 2023
I have to write a code to solve numerical integration using the trapezoidal method. The code is correct, however when I put the ranges from -1 to 1 the result is Inf Inf Inf Inf Inf Inf Inf ..
Range and function cannot be changed...
does anyone know how to solve this problem?
Thanks in advance to everyone who tries to help.
f= @ (x) 1/sqrt((1-x^2)+(e/2)*(1-x^4));
Err = abs (II(n)-II(n-1));
disp (['Result is:' num2str(II(n)) ', Accuracy is:' num2str(Err) ', Number of Interactions is:' num2str(n)])
xs = linspace(a,b,n+1);
xlabel ('variation of x');
Torsten on 30 Jan 2023
Edited: Torsten on 30 Jan 2023
e = 0.01;
f = 1/sqrt((1-x^2)+(e/2)*(1-x^4));
value = vpaintegral(f,-1,1)
f = @(x)1./sqrt((1-x.^2)+(e/2)*(1-x.^4));
x = linspace(-1+eps,1-eps,250000000);
fx = f(x);
value = trapz(x,fx)
Steven Lord on 30 Jan 2023
Look at your function over the region in question.
f= @(x) 1./sqrt((1-x.^2)+(e/2)*(1-x.^4));
fplot(f, [-2 2])
Your function goes to infinity at the endpoints.
So if the value of the integral is not inf, what do you expect it should be and why do you expect it to be that value?
Actually, the problem does have a finite solution. Note that this integral, which for small value of e, is quite similar,
does have a well known analytical solution.
You cannot use trapezoidal rule to do the integration exactly fully to the endpoints, since the shape of that function is too extreme near the endpoints. Singularities are a real pain in the neck. If the endpoints are FIXED AND immovable, AND you need to use trapezoidal rule, then you have a problem.
Instead, you try to do it symbolically, at least to see what the correct answer is.
K = 1/sqrt((1-x^2)+(e/2)*(1-x^4));
So int fails to find an analytical solution. That is not completely surprising. But vpaintegral does find a solution. Is that number something we can trust?
When e is small, I would expect the result to be close to this:
which does have a known result at pi. So for small e, indeed we did get a number close to pi. So I do trust the value predicted by vpaintegral.
Could you have solved it using trapezoidal integration?
You could try, if you were allowed to change the limits by even a tiny amount. For example,
format long g
f = @(x) 1./sqrt((1-x.^2)+(e/2)*(1-x.^4));
xt = linspace(-1 + eps,1 - eps,1000);
But this fails. The trapezoidal rule is simply unable to approximate the extreme shape of that function near the limits. That does not surprise me at all. In fact, I assume this is why you were assigned this problem as homework. My guess is that even a Simpson's quadrature (or any higher order Newton-Cotes rule) will fail here too. They are all based on polynomials, and polynomials abhor singularities.
One idea (In know you are probably supposed to do this using trapezoidal intgration. But that is your homework, not mine.) I might wonder if a transformation might help, of this general form:
x = sin(u);
If we try that,
I = subs(K,x,sin(u))
This changes the limits of integration from [-1,1], to [-pi/2,pi/2]
And that also fails. Well, it was a thought. ;-)
But there are alternatives.
integral(f,-1 + eps,1-eps)
So integral also survives. It is only trapezoidal rule that fails miserably. I might not be surprised if perhaps a Gauss-Chebyshev quadrature (the first kind) might also work. But that is way beyond anything you were asked to do.
In the latter link, we see the nodes and weights for a first kind Gauss-Chebyshev quadrature, are simple. For an n-node rule, they are
cheby1nodes = @(N) cos((2*(1:N) - 1)/(2*N)*pi);
cheby1weights = @(N) ones(1,N)*pi/N;
With 5 nodes, for example
n = 5;
nodes = cheby1nodes(n)
As you can see, this did very well.
Could you have used other methods? Of course. We could even have solve it using an ODE solver. Perhaps ODE15s, or something like that. But not trapezoidal rule.
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!