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));
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?
When e is small, I would expect the result to be close to this:
int(1/sqrt(1-x^2),[-1,1])
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,
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))
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)
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
nodes = cheby1nodes(n)
nodes =
0.951056516295154 0.587785252292473 6.12323399573677e-17 -0.587785252292473 -0.951056516295154
dot(cheby1weights(n),f(nodes).*sqrt(1-nodes.^2))
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.
