"Matthew " <mwk5v@virginia.edu> wrote in message <gl8cd6$4rj$1@fred.mathworks.com>...
> Hello everyone,
>
> I'm trying to evaluate a symbolic integral in Matlab whose solution is a linear combination of two elliptic integrals, but Matlab does not appear to evaluate the integral correctly. The integral that I am trying to evaluate is:
>
> E = sin^2(x)*sqrt(sin^2(x) + m*cos^2(x)), 0 <= m <= 1
>
> I want to integrate with respect to x, for x = 0..pi.
>
> At the limits of m, the solution is in closedform and becomes a trivial integral
>
> m=0: E = sin^2(x)*sqrt(sin^2(x)) = sin^3(x), whose integral from x=0..pi is simply 4/3.
>
> m=1: E = sin^2(x)*sqrt[sin^2(x)+cos^2(x)] = sin^2(x) since sin^2(x)+cos^2(x)=1, and the integral is simply pi/2.
>
> However, when I take the integral in Matlab, I get the following expression:
>
> int(E) = 2/3*[m*K(g)+(m2)*E(g)]/(m1), g=sqrt(1m)
>
> where K(g) is the complete integral of the first kind, and E(g) is the complete integral of the second kind. Substituting m=0 into this expression gives the correct value for int(E) = 4/3, but substituting m=0.99999999999 into this expression (since m1 would yield a division by zero error), gives a very large value. Indeed, as m>1, int(E) > Inf. I haven't been able to get anywhere on this problem for the last few days, and any assistance the community can provide would be greatly appreciated.
>
> Thanks,
> Matt
I believe you will find that the expression matlab gave is actually correct, even though as m approaches 1, it becomes increasingly difficult to calculate accurately. It is analogous to the integral of sin(x)/x as x approaches 0. That approaches a definite limit, namely 1, but becomes increasingly difficult to calculate directly. Below the point x = 0.00000001 the computation becomes increasingly inaccurate.
In your expression as m approaches 1, the quantity 1m in the denominator approaches zero, but K(g) and E(g) both approach pi/2 as g approaches zero which make the numerator also approach zero. One can determine this limit using power series. We have
K(g) = pi/2*(1+1/4*g^2 + higher powers of g)
E(g) = pi/2*(11/4*g^2 + higher powers of g)
(See http://en.wikipedia.org/wiki/Elliptic_integral)
Substituting just the terms up to g^2 in the expression gives
2/3*(m*K(g)+(m2)*E(g))/(m1) =
pi/3*((1g^2)*(1+1/4*g^2)(1+g^2)*(11/4*g^2))/(g^2) =
pi/3*(2*g^2+1/2*g^2)/(g^2) =
pi/3*3/2 = pi/2
If the higher powers of g were included, their terms would approach zero, so this is the true limit as m approaches 1, and it agrees with your earlier determination.
If you were to substitute m = .99999999, you should get a reasonably accurate approximation of pi/2 (provided you have accurate elliptic integrals.) Closer than that causes increasingly large errors which finally become wildly inaccurate, just as in the case of sin(x)/x.
Roger Stafford
