numerical integration inf to inf

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;

Answer 1

[Torsten]

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

Answer 2

[Steven Lord]

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?

[John D'Errico]

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.

Answer 3

[John D'Errico]

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.

https://en.wikipedia.org/wiki/Gaussian_quadrature

https://en.wikipedia.org/wiki/Chebyshev–Gauss_quadrature

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.