How to resolve the Plotting issue in my code?

Question

Athira T Das 2023-4-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1953814-how-to-resolve-the-plotting-issue-in-my-code

回答： Walter Roberson 2023-4-27

在 MATLAB Online 中打开

I done the partial differentiation of W and W3 and then substitute with an array of values and tried to plot.

clc;close all; clear all;

w=-3;T=38;S=28;Om =1; w0 =0.002; sigma0 = 0.04;muu0=4*pi*10^-7;

lambda=532*10^-9;k=2*pi/lambda;z=100;epsilon0=8.85*10^-12;

alphac=0.02;T=3.609699516774368e+03;

syms r1x r1y r2x r2y

J = 1i*k/(2*z) - T;

J1 = conj(J) ;

A = 1/(2*w0) + 1/(2*sigma0^2) + T;

Bx = (1i*k/2*z)*(r1x+r2x) - T*(r1x - r2x);By = (1i*k/2*z)*(r1y+r2y) - T*(r1y - r2y);

C = 2/(w0^2) + k^2/(4*A*z^2);

Dx = (1i*k/z)*(r1x-r2x) + 2*Om; Dy = (1i*k/z)*(r1y-r2y) + 2*Om;

Dx1 = (1i*k/z)*(r1x-r2x) - 2*Om; Dy1 = (1i*k/z)*(r1y-r2y) - 2*Om;

Ex = 0.5*(Dx - (1i*k*Bx)/(2*A*z)); Ey = 0.5*(Dy - (1i*k*By)/(2*A*z));

Ex1 = 0.5*(Dx1 - (1i*k*Bx)/(2*A*z)); Ey1 = 0.5*(Dy1 - (1i*k*By)/(2*A*z));

Fx=Bx + Om;Fy = By + Om;

Fx1=Bx - Om;Fy1 = By - Om;

Gx=0.5*((1i*k/z)*(r1x-r2x) - (1i*k*Fx/2*A*z) );Gy=0.5*((1i*k/z)*(r1y-r2y) - (1i*k*Fy/2*A*z) );

Gx1=0.5*((1i*k/z)*(r1x-r2x) - (1i*k*Fx1/2*A*z) );Gy1=0.5*((1i*k/z)*(r1y-r2y) - (1i*k*Fy1/2*A*z) );

H = A + (k^2*w0^2)/(8*z^2);

Ix = 0.5*(Bx - (1i*k*w0*w0*Dx)/(4*z));Iy = 0.5*(By - (1i*k*w0*w0*Dy)/(4*z));

Ix1 = 0.5*(Bx - (1i*k*w0*w0*Dx1)/(4*z));Iy1 = 0.5*(By - (1i*k*w0*w0*Dy1)/(4*z));

Jx = 0.5*(Fx + (k*k*w0*w0*(r1x-r2x))/(4*z*z));Jy = 0.5*(Fy + (k*k*w0*w0*(r1y-r2y))/(4*z*z));

Jx1 = 0.5*(Fx1 + (k*k*w0*w0*(r1x-r2x))/(4*z*z));Jy1 = 0.5*(Fy1 + (k*k*w0*w0*(r1y-r2y))/(4*z*z));

M1 =(( Ex^2 + Ey^2 )/(C^2) + (1/C) - (Ix^2 + Iy^2 )/(4*H^2) - (1/4*H) + (1i*(Ix*Ey - Ex*Iy)/(C*H)))*exp(((Bx^2 + By^2)/(4*A)) + ((Ex^2 + Ey^2)/C));

M2 =(( Ex1^2 + Ey1^2 )/(C^2) + (1/C) - (Ix1^2 + Iy1^2 )/(4*H^2) - (1/4*H) + (1i*(Ix1*Ey1 - Ex1*Iy1)/(C*H)))*exp(((Bx^2 + By^2)/(4*A)) + ((Ex1^2 + Ey1^2)/C));

M3 =(( Gx^2 + Gy^2 )/(C^2) + (1/C) - (Jx^2 + Jy^2 )/(4*H^2) - (1/4*H) + (1i*(Jx*Gy - Gx*Jy)/(C*H)))*exp(((Fx^2 + Fy^2)/(4*A)) + ((Gx^2 + Gy^2)/C));

M4 =(( Gx1^2 + Gy1^2 )/(C^2) + (1/C) - (Jx1^2 + Jy1^2 )/(4*H^2) - (1/4*H) + (1i*(Jx1*Gy1 - Gx1*Jy1)/(C*H)))*exp(((Fx1^2 + Fy1^2)/(4*A)) + ((Gx1^2 + Gy1^2)/C));

P = (k*k)/(16*A*C*z*z);

W=P*exp(conj(J)*(r1x^2+r1y^2) + (J)*(r2x^2+r2y^2) + 2*T*(r1x*r2x + r1y*r2y))*(M1 + M2 - M3 - M4);

dW1=diff(W,r2x);

dW2=diff(W,r2y);

W3 = r1y*dW1 - r1x*dW2;

r=linspace(-1,1,100);

q1 =(subs(W3,{r1x,r1y,r2x,r2y},{r,r,r,r}));

q2 =(subs(W,{r1x,r1y,r2x,r2y},{r,r,r,r}));

Lorb = imag(q1)*(-epsilon0/k);

S=(k/muu0*w0)*q2;

hw = 2*pi*3*10^8*6.62607015 * 10^-34/(2*pi*lambda);

lorb = hw.*Lorb./S;

plot(r,lorb)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2023-4-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1953814-how-to-resolve-the-plotting-issue-in-my-code#answer_1224179

your values all come out nan. You have terms with really silly exp like

exp(465047415018097955185175015559120439814522795787376405109289013546782463990466768163364960497178088168719688701/146734163107013360162621761657081706047762220035580747699952030138169819136)

which is like exp(3e36)