Please help, I have no idea how to do this integral in matlab
2 views (last 30 days)
Show older comments
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1528251/image.png)
2 Comments
Walter Roberson
on 2 Nov 2023
No closed form solution. Runs the risk of having a singularity if a_0 > theta_0
Accepted Answer
Torsten
on 2 Nov 2023
Edited: Torsten
on 3 Nov 2023
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.
More Answers (1)
Walter Roberson
on 2 Nov 2023
Edited: Walter Roberson
on 3 Nov 2023
syms a_0 positive
syms theta real; assumeAlso(theta >= 0);
theta_0 = sym(30.5);
F = 1/(sind(theta_0) - sind(theta))
a_0 = sym([5:5:25 26 27 28 29]).';
values = arrayfun(@(A) int(F, theta, 0, A), a_0);
values(1) %for example
valued = double(values);
plot(a_0, valued)
4 Comments
Walter Roberson
on 3 Nov 2023
@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))
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
valued = double(values);
plot(a_0, real(valued), a_0, imag(valued))
legend({'real', 'imaginary'})
isimag = imag(valued) ~= 0;
a_0(isimag)
valued(isimag)
Torsten
on 3 Nov 2023
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 = sol * 180/pi
vpa(sol)
See Also
Categories
Find more on Number Theory in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!