matlab does not give any answer.

56 次查看（过去 30 天）
Anal
Anal about 23 hours 前
编辑： Walter Roberson about 7 hours 前
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 about 9 hours 前
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 Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by