Please help, I have no idea how to do this integral in matlab

[Walter Roberson]

No closed form solution. Runs the risk of having a singularity if a_0 > theta_0

[Laura]

If theta_0 were 30.5, could it be done?

Answer 1

[Torsten]

The integral seems to exist for those values of theta0 with cos(theta0) ~= 0, thus theta0 ~= (2*k+1)*pi/2 for k in Z.

In this case, the Taylor expansion of sin(theta0)-sin(theta) around theta0 shows a behaviour like 1/sqrt(x) around x=0.

If cos(theta0) = 0, the Taylor expansion of sin(theta0)-sin(theta) around theta0 shows a behaviour like 1/x around x = 0 which means that the integral diverges.

[Laura]

If theta_0 were 30.5, could it be done?

[Torsten]

Seems to work:

syms theta theta0
theta0 = 30.5*pi/180;
f = 1/sqrt(sin(theta0)-sin(theta));
sol = vpaintegral(f,theta,0,theta0)

sol = 

1.54036

[Laura]

thank youu. You saved me :)

Answer 2

[Walter Roberson]

syms a_0 positive
syms theta real; assumeAlso(theta >= 0);
theta_0 = sym(30.5);
F = 1/(sind(theta_0) - sind(theta))

F = 

a_0 = sym([5:5:25 26 27 28 29]).';
values = arrayfun(@(A) int(F, theta, 0, A), a_0);
values(1)    %for example

ans = 

valued = double(values);
plot(a_0, valued)

[Torsten]

@Walter Roberson

I think you forgot the sqrt.

[Laura]

thank you for help :)

[Walter Roberson]

@Torsten you are right, I did forget the sqrt() !

It looks like MATLAB is able to integrate even so, but with a more complicated formula.

The imaginary components of the results are so small that they are surely due to round-off error in the numeric calculations.

format long g
syms a_0 positive
syms theta real; assumeAlso(theta >= 0);
theta_0 = sym(30.5);
F = 1/sqrt(sind(theta_0) - sind(theta))

F = 

a_0 = sym([5:5:25 26 27 28 29 30:.1:30.4]).';
values = arrayfun(@(A) int(F, theta, 0, A), a_0);
values(1)    %for example

ans = 

valued = double(values);
plot(a_0, real(valued), a_0, imag(valued))
legend({'real', 'imaginary'})

isimag = imag(valued) ~= 0;
a_0(isimag)

ans = 

valued(isimag)

ans = 
           7.34891328498187 +  1.60441675001006e-67i
           15.5050816756634 +  1.20331256250754e-67i
           24.7962226985339 -  6.85758772181718e-68i
           35.8443076292128 -   3.0831418652267e-67i
           50.1801693418716 -  4.50826319501674e-67i
           53.7878830561832 -  4.49347594386457e-67i
           57.8335596505002 -  4.34745183873693e-67i
           62.5233721908558 -  4.02767753257133e-67i
            68.307102046422 -  3.52121418060963e-67i
           76.7289013342626 -  2.78739684218337e-67i
           77.9450026132779 -  2.70237014805841e-67i
            79.325691819995 -  2.60994982835737e-67i
           80.9638641965888 -   2.5230747278384e-67i
           83.0993245089765 -  2.43250281453138e-67i

[Torsten]

Your values seem to converge towards

syms theta theta0
theta0 = 30.5*pi/180;
f = 1/sqrt(sin(theta0)-sin(theta));
sol = vpaintegral(f,theta,0,theta0)

sol = 

1.54036

sol = sol * 180/pi

sol = 

vpa(sol)

ans = 

88.256285088573850083615798717037