matlab does not give any answer.

Here is my code:

tic
clc; all clear;

a=18897.2598;   
w1=2*a;   % in a0
w2=3.78*a; % in a0
lambda=780*0.001*a;  %in a0
z_R1 =(pi/lambda).*(w1.^2);  % in micrometer
z_R2 =(pi/lambda).*(w2.^2);
R=1:1000:30001; % in a0
angle_CM=0:1:1;
[R_CM,Cos_theta_R] = meshgrid(R,angle_CM);
l_S=0; l_P=1;  s=1/100;
syms r
syms theta
for i=1
    for j=1:2
        
        Part_1_1(i,j)= vpa(1+((1./z_R1).*(R_CM(i,j).*Cos_theta_R(i,j)+r.*cos(theta))).^2);
        Part_2_1(i,j)= vpa(exp(((-1/(w1.^2)).*(R_CM(i,j).^2).*(1-(Cos_theta_R(i,j).^2)))./Part_1_1(i,j).^2));
        Part_3_1(i,j)= vpa(exp(((-1/(w1.^2)).*(r.^2).*((sin(theta)).^2))./Part_1_1(i,j).^2));
        Part_4_1(i,j) =vpa(exp(((-1/(w1.^2)).*2.*R_CM(i,j).*r.*sin(theta).*sqrt(1-Cos_theta_R(i,j).^2))./Part_1_1(i,j).^2));
        Total_1(i,j)=vpa(((Part_1_1(i,j).^(-0.5)).*(Part_2_1(i,j).*Part_3_1(i,j).*Part_4_1(i,j))));
        
        Part_1_2(i,j)= vpa(1+((1./z_R2).*(R_CM(i,j).*Cos_theta_R(i,j)+r.*cos(theta))).^2);
        Part_2_2(i,j)= vpa(exp(((-1/(w2.^2)).*(R_CM(i,j).^2).*(1-(Cos_theta_R(i,j).^2)))./Part_1_2(i,j).^2));
        Part_3_2(i,j)= vpa(exp(((-1/(w2.^2)).*(r.^2).*((sin(theta)).^2))./Part_1_2(i,j).^2));
        Part_4_2(i,j) =vpa(exp(((-1/(w2.^2)).*2.*R_CM(i,j).*r.*sin(theta).*sqrt(1-Cos_theta_R(i,j).^2))./Part_1_2(i,j).^2));
        Total_2(i,j)=vpa(((Part_1_2(i,j).^(-0.5)).*(Part_2_2(i,j).*Part_3_2(i,j).*Part_4_2(i,j))));
        
        Total(i,j) =vpa((abs((Total_1(i,j)-Total_2(i,j)))).^2);
        del_0_S=3.1311806; del_2_S=0.1786; del_4_S=0; del_6_S=0; del_8_S=0; % For n S_1/2 state of Cs, Quantum defect parameter
        n_S=10;
        n_star_S=double(n_S-del_0_S-(del_2_S./((n_S-del_0_S).^2))-(del_4_S./((n_S-del_0_S).^4))...
            -(del_6_S./((n_S-del_0_S).^6))-(del_8_S./((n_S-del_0_S).^8)));           % Effective n
        W_i= (whittakerW(n_star_S, l_S+0.5,(2.*(r))./(n_star_S)));
        Gii=  W_i./((sqrt(n_star_S.^2.*gamma(n_star_S+l_S+1).*gamma(n_star_S-l_S))));
        
        %include angular part of S1/2, i.e.,
        %Yj,MJ(theta,phi)=CG*Yl,ml(theta,phi)=((l+mj+0.5)/(2l+1))^0.5.*Yl,ml(theta,phi)
        %For S1/2 state anular part=1*Y0,0(theta,phi)=Y0,0(theta,phi)=0.5*sqrt(1/pi)
        
        fun_S_S(i,j)= 2*pi*(r.^2).*Gii.*Gii.*sin(theta).^2.*Total(i,j).*((0.5*sqrt(1/pi)).^2);
        r_min=s*n_S*n_S/(n_S+n_S);
        
        Matrix_element(i,j)=double(int(int(fun_S_S(i,j), 0, pi),r_min,1000));
        
    end
end
toc

Answer 1

[Suraj Kumar]

Hi @Anal,

Based on my understanding, the issue in the code is due to the symbolic integration of the function 'fun_S_S'.

The function fun_S_S involves several components, including trigonometric functions, polynomial terms, and potentially complex products. And symbolic integration attempts to find an exact analytical solution, which can be infeasible for complex expressions.

So as a workaround, you can use matlabFunction which transforms the symbolic expression into a MATLAB function handle, allowing for efficient numerical evaluation.This can be followed by numerical integration using the integral2 function in MATLAB.

You can refer to the attached code snippet for better understanding:

fun_S_S(i,j)= 2*pi*(r.^2).*Gii.*Gii.*sin(theta).^2.*Total(i,j).*((0.5*sqrt(1/pi)).^2);
r_min=s*n_S*n_S/(n_S+n_S);
fun_S_S_numeric = matlabFunction(fun_S_S(i,j), 'Vars', [r, theta]);
Matrix_element(i,j) = integral2(fun_S_S_numeric, r_min, 1000, 0, pi);

To learn more about matlabFunction and integral2 functions in MATLAB, you can refer to the following documentations:

• https://www.mathworks.com/help/symbolic/sym.matlabfunction.html
• https://www.mathworks.com/help/matlab/ref/integral2.html

Happy Coding!