Problem using Matlab to do Elliptic Integrals

Hello everyone, I have problem doing two elliptic integrals in the for loop. The basic equations are as shown in the attached picture.

I have 3 variables, x, y, and v0. I want to plot the relationship curve of x and y.

My code is list below:

-----------------------------------------------------

 syms fi;
v = [0:0.001:pi/2];
L = length(v); 
y = zeros(1,L);
k = zeros(1,L);
for i = 1:L  
    t1 = 1./sqrt(sin(v(i))-sin(fi)); 
    m1 = int(t1,0,v(i));
    k(i) = vpa(m1)
    t2 = sin(fi)/(sin(v(i))-sin(fi));
    m2 = 1/k(i)*int(t1,0,v(i));
    y(i) = vpa(m2)
    f(i) = k(i)^2./2;
end
plot(f,y,'b');

--------------------------------

The error massage shows I use wrongly of vpa,

-----------------------------------------------------

The following error occurred converting from sym to double: Error using mupadmex Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.

If the input expression contains a symbolic variable, use the VPA function instead.

Error in NonlinearStiffness (line 10) k(i) = vpa(m1)

-------------------------------------------------

However, when I just run the code for once, i.e. I give the exact value to v(i), it can give me a answer.

--------------------------------------------------

>> t1 = 1./sqrt(sin(0.2)-sin(fi));

   m1 = int(t1,0,0.2);

>> m1

m1 =

(7217745006463825^(1/2)*9007199254740992^(1/2)*ellipticF(pi/4, 18014398509481984/7217745006463825)*2*i)/7217745006463825 - (2*7217745006463825^(1/2)*ellipticF(pi/4 - 1/10, 18014398509481984/7217745006463825)*(9007199254740992*sin(1/5) - 1789454248277167)^(1/2))/(7217745006463825*(1789454248277167/9007199254740992 - sin(1/5))^(1/2))

>> vpa(m1)

ans =

0.901 + 1.145e-36*i

------------------------------------------------

Why the code is wrong in the for loop?

Why the result has a imaginary part "i"?

Answer 1

[Walter Roberson]

You should try again at 0.0 . Although your integral is logically 0 there as you would be integrating from 0 to 0, when you use int() the indefinite integral is done and then the bounds are substituted in, so if there is a singularity in the indefinite integral at the limit point then you might not get back 0.

The integral works out to be strictly real-valued over 0 to Pi, but the indefinite integral has to take into account that you may be integrating in places where the sin() is negative. Accordingly the indefinite integral is a formula in complex numbers which theoretically works out as real for the range you are interested in. but only theoretically: due to round-off, the imaginary component does not always completely vanish.

Answer 2

[Zhou Xi]

Hi Walter,

Thank you very much for your help.

I changed the 0 to 0.0 according to your suggestion but the problem is still there.

syms fi;
v = [0:0.001:pi/2];
L = length(v); 
y = zeros(1,L);
k = zeros(1,L);
for i = 1:L  
    t1 = 1./sqrt(sin(v(i))-sin(fi)); 
    m1 = int(t1,0.0,v(i));
%     m1 = int(1./sqrt(sin(v(i))-sin(fi)),fi,0,v(i));
    k(i) = vpa(m1)
    t2 = sin(fi)/(sin(v(i))-sin(fi));
    m2 = 1/k(i)*int(t1,0.0,v(i));
%     m2 = 1/k(i)*int(sin(fi)/(sin(v(i))-sin(fi)),fi,0,v(i));
    y(i) = vpa(m2)
    f(i) = k(i)^2./2;
end
plot(f,y,'b');

The error message is:

Warning: Explicit integral could not be found. 
The following error occurred converting from sym to double:
Error using mupadmex
Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.

If the input expression contains a symbolic variable, use the VPA function instead.

Error in NonlinearStiffness (line 10)
    k(i) = vpa(m1)

Before it went wrong, the command window show the elements of the k as below. I have no idea why it shows such information.

    k =

      Columns 1 through 18

         0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0

      Columns 19 through 36

         0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0

      Columns 37 through 54
  .................
    Columns 1561 through 1565

          0                  0                  0                  0                  0          

    Columns 1566 through 1570

          0                  0                  0                  0                  0          

    Column 1571

          0

[Walter Roberson]

Your command window is showing k because you have

   k(i) = vpa(m1)

with no semi-colon.

You need to give the command

dbstop if error

and then run the code. When it stops, examine the value of v(i) and tell us what it is.