numerical integration inf to inf

2 views (last 30 days)
PTK
PTK on 30 Jan 2023
Edited: John D'Errico on 30 Jan 2023
Hello guys
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.
clear all;
close all;
clc
e=0.01;
err=0.0001;
f= @ (x) 1/sqrt((1-x^2)+(e/2)*(1-x^4));
a=-1;
b=1;
for n=1:100
dx=(b-a)/n;
S1=f(a);
Sn=f(b);
Sab=0;
for i=a+dx:dx:b-dx
Sab=Sab+2*f(i);
end
I=(dx/2)*(S1+Sab+Sn);
II(n)=I;
if n>=2
Err = abs (II(n)-II(n-1));
if Err<=err
disp (['Result is:' num2str(II(n)) ', Accuracy is:' num2str(Err) ', Number of Interactions is:' num2str(n)])
break;
end
end
xs = linspace(a,b,n+1);
end
plot (II);
xlabel ('variation of x');
ylabel('Integration Value');
title('Trapezoid Method');
grid on;

Answers (3)

Torsten
Torsten on 30 Jan 2023
Edited: Torsten on 30 Jan 2023
Set
a=-1+eps;
b=1-eps;
instead of
a=-1;
b=1;
For comparison:
syms x
e = 0.01;
f = 1/sqrt((1-x^2)+(e/2)*(1-x^4));
value = vpaintegral(f,-1,1)
value = 
3.12988
Or:
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)
value = 3.5074

Steven Lord
Steven Lord on 30 Jan 2023
Look at your function over the region in question.
e=0.01;
f= @(x) 1./sqrt((1-x.^2)+(e/2)*(1-x.^4));
a=-1;
b=1;
fplot(f, [-2 2])
Your function goes to infinity at the endpoints.
f([-1 1])
ans = 1×2
Inf Inf
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?
  1 Comment
John D'Errico
John D'Errico on 30 Jan 2023
Edited: John D'Errico on 30 Jan 2023
Actually, the problem does have a finite solution. Note that this integral, which for small value of e, is quite similar,
syms x
int(1/sqrt(1-x^2),[-1,1])
ans = 
π
does have a well known analytical solution.

Sign in to comment.


John D'Errico
John D'Errico on 30 Jan 2023
Edited: John D'Errico on 30 Jan 2023
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.
syms x
e=0.01;
K = 1/sqrt((1-x^2)+(e/2)*(1-x^4));
int(K,[-1,1])
ans = 
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?
vpaintegral(K,[-1,1])
ans = 
3.12988
When e is small, I would expect the result to be close to this:
int(1/sqrt(1-x^2),[-1,1])
ans = 
π
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);
trapz(xt,f(xt))
ans =
94532.8315996745
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,
syms u
I = subs(K,x,sin(u))
I = 
This changes the limits of integration from [-1,1], to [-pi/2,pi/2]
int(I,u,[-pi/2,pi/2])
ans = 
And that also fails. Well, it was a thought. ;-)
But there are alternatives.
integral(f,-1 + eps,1-eps)
ans =
3.12988115106573
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)
nodes = 1×5
0.951056516295154 0.587785252292473 6.12323399573677e-17 -0.587785252292473 -0.951056516295154
dot(cheby1weights(n),f(nodes).*sqrt(1-nodes.^2))
ans =
3.1298811510667
As you can see, this did very well.
vpaintegral(K,[-1,1],'reltol',1e-30)
ans = 
3.12988115084129612486746309297
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!