matlab does not give any answer.

조회 수: 55 (최근 30일)
Anal
Anal 대략 19시간 전
편집: Walter Roberson 대략 2시간 전
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

답변 (1개)

Suraj Kumar
Suraj Kumar 대략 4시간 전
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:
Happy Coding!

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by