필터 지우기
필터 지우기

Not getting finite result

조회 수: 1 (최근 30일)
Athira T Das
Athira T Das 2022년 6월 9일
댓글: Athira T Das 2022년 6월 15일
clear all
clc
syms theta
syms z K
L=40
lambda=0.532*10^(-6)
a=10*(lambda*L)^0.5
A=1
n=0
m=0
k=2*pi/(lambda)
e =(k*a)/((k*a^2)-(1i*L))
ee =(k*a)/((k*a^2)+(1i*L))
G=A*ee*a*hermiteH(n,0)*hermiteH(m,0)
f=norm(G)
B=(1i*ee*(L-z)/(2*k^2*a))*((k*a^2)+(1i*z))*K^2
Hn=hermiteH(n,(z-L)/k*e*K*cos(theta))
Hm=hermiteH(m,(z-L)/k*e*K*sin(theta))
N=A*1i*k*a*ee*exp(B)*Hn*Hm
NNN=subs(N, 1i, -1i)
NN=subs(N, k, -k)
W1=N*NN/G^2
W2=N*NNN/f^2
etta=1*10^-3
delta=8.284*(K*etta)^(4/3)+12.978*(K*etta)^2
chiT=10^-10
w=-2
epsilon=10^-1
ff=chiT*(exp(-1.863*10^-2*delta)+w.^(-2)*exp(-1.9*10^-4*delta)-2*w.^(-1)*exp(-9.41*10^-3*delta))
Phi=0.388*10^(-8)*epsilon^(-1/3)*K^(-11/3)*(1+2.35*(K*etta)^(2/3))*ff
F=K*(W1+W2)*Phi
fun = matlabFunction(F)
q1 = int(F, theta, 0, 2*pi)
q2 = int(K*q1, K, 0, 1)
q3 = int(q2, z, 0, L)
m=4*pi*real(q3)
  댓글 수: 5
Torsten
Torsten 2022년 6월 13일
Please include your code.
Athira T Das
Athira T Das 2022년 6월 13일
clear all
clc
syms theta
syms z K
L=40
L = 40
lambda=0.532*10^(-6)
lambda = 5.3200e-07
a=10*(lambda*L)^0.5
a = 0.0461
A=1
A = 1
n=0
n = 0
m=0
m = 0
k=2*pi/(lambda)
k = 1.1810e+07
e =(k*a)/((k*a^2)-(1i*L))
e = 21.6777 + 0.0345i
ee =(k*a)/((k*a^2)+(1i*L))
ee = 21.6777 - 0.0345i
G=A*ee*a*hermiteH(n,0)*hermiteH(m,0)
G = 1.0000 - 0.0016i
f=norm(G)
f = 1.0000
B=(1i*ee*(L-z)/(2*k^2*a))*((k*a^2)+(1i*z))*K^2
B = 
Hn=hermiteH(n,(z-L)/k*e*K*cos(theta))
Hn = 
1
Hm=hermiteH(m,(z-L)/k*e*K*sin(theta))
Hm = 
1
N=A*1i*k*a*ee*exp(B)*Hn*Hm
N = 
NNN=subs(N, 1i, -1i)
NNN = 
NN=subs(N, k, -k)
NN = 
W1=N*NN/G^2
W1 = 
W2=N*NNN/f^2
W2 = 
etta=1*10^-3
etta = 1.0000e-03
delta=8.284*(K*etta)^(4/3)+12.978*(K*etta)^2
delta = 
chiT=10^-10
chiT = 1.0000e-10
w=-2
w = -2
epsilon=10^-1
epsilon = 0.1000
ff=chiT*(exp(-1.863*10^-2*delta)+w.^(-2)*exp(-1.9*10^-4*delta)-2*w.^(-1)*exp(-9.41*10^-3*delta))
ff = 
Phi=0.388*10^(-8)*epsilon^(-1/3)*K^(-11/3)*(1+2.35*(K*etta)^(2/3))*ff
Phi = 
F=K*(W1+W2)*Phi
F = 
fun = matlabFunction(F)
fun = function_handle with value:
@(K,z)1.0./K.^(8.0./3.0).*((K./1.0e+3).^(2.0./3.0).*(4.7e+1./2.0e+1)+1.0).*(exp(K.^2.*(z.*(2.680902008955589e-15+1.6844604112658e-12i)+-1.072360803582236e-13-6.737841645063198e-11i).*(z.*1i+pi.*8.0e+3).*-2.0).*(1.394878794885149e+14+1.145503965346681e-4i)-exp(K.^2.*(z.*(2.680902008955589e-15-1.6844604112658e-12i)+-1.072360803582236e-13+6.737841645063198e-11i).*(z.*1i-pi.*8.0e+3)).*exp(-K.^2.*(z.*(2.680902008955589e-15+1.6844604112658e-12i)+-1.072360803582236e-13-6.737841645063198e-11i).*(z.*1i+pi.*8.0e+3)).*1.394878794885149e+14).*(exp((K./1.0e+3).^(4.0./3.0).*(-1.57396e-3)-K.^2.*2.46582e-9)./4.0e+10+exp((K./1.0e+3).^(4.0./3.0).*(-7.795244e-2)-K.^2.*1.2212298e-7)./1.0e+10+exp((K./1.0e+3).^(4.0./3.0).*(-1.5433092e-1)-K.^2.*2.4178014e-7)./1.0e+10).*(-8.359206597323707e-9)
r=integral3(fun,0,2*pi,0,1,0,L)
Error using symengine>@(K,z)1.0./K.^(8.0./3.0).*((K./1.0e+3).^(2.0./3.0).*(4.7e+1./2.0e+1)+1.0).*(exp(K.^2.*(z.*(2.680902008955589e-15+1.6844604112658e-12i)+-1.072360803582236e-13-6.737841645063198e-11i).*(z.*1i+pi.*8.0e+3).*-2.0).*(1.394878794885149e+14+1.145503965346681e-4i)-exp(K.^2.*(z.*(2.680902008955589e-15-1.6844604112658e-12i)+-1.072360803582236e-13+6.737841645063198e-11i).*(z.*1i-pi.*8.0e+3)).*exp(-K.^2.*(z.*(2.680902008955589e-15+1.6844604112658e-12i)+-1.072360803582236e-13-6.737841645063198e-11i).*(z.*1i+pi.*8.0e+3)).*1.394878794885149e+14).*(exp((K./1.0e+3).^(4.0./3.0).*(-1.57396e-3)-K.^2.*2.46582e-9)./4.0e+10+exp((K./1.0e+3).^(4.0./3.0).*(-7.795244e-2)-K.^2.*1.2212298e-7)./1.0e+10+exp((K./1.0e+3).^(4.0./3.0).*(-1.5433092e-1)-K.^2.*2.4178014e-7)./1.0e+10).*(-8.359206597323707e-9)
Too many input arguments.

Error in integral3>@(y,z)fun(x(1)*ones(size(z)),y,z) (line 129)
@(y,z)fun(x(1)*ones(size(z)),y,z), ...

Error in integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral3>innerintegral (line 128)
Q1 = integral2Calc( ...

Error in integral3>@(x)innerintegral(x,fun,yminx,ymaxx,zminxy,zmaxxy,integral2options) (line 111)
f = @(x)innerintegral(x, fun, yminx, ymaxx, ...

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral3 (line 113)
Q = integralCalc(f,xmin,xmax,integralOptions);
m=4*pi*real(r)

댓글을 달려면 로그인하십시오.

채택된 답변

Torsten
Torsten 2022년 6월 13일
편집: Torsten 2022년 6월 13일
syms theta
syms z K
L=40;
lambda=0.532*10^(-6);
a=10*(lambda*L)^0.5;
A=1;
n=0;
m=0;
k=2*pi/(lambda);
e =(k*a)/((k*a^2)-(1i*L));
ee =(k*a)/((k*a^2)+(1i*L));
G=A*ee*a*hermiteH(n,0)*hermiteH(m,0);
f=norm(G);
B=(1i*ee*(L-z)/(2*k^2*a))*((k*a^2)+(1i*z))*K^2;
Hn=hermiteH(n,(z-L)/k*e*K*cos(theta));
Hm=hermiteH(m,(z-L)/k*e*K*sin(theta));
N=A*1i*k*a*ee*exp(B)*Hn*Hm;
NNN=subs(N, 1i, -1i);
NN=subs(N, k, -k);
W1=N*NN/G^2;
W2=N*NNN/f^2;
etta=1*10^-3;
delta=8.284*(K*etta)^(4/3)+12.978*(K*etta)^2;
chiT=10^-10;
w=-2;
epsilon=10^-1;
ff=chiT*(exp(-1.863*10^-2*delta)+w.^(-2)*exp(-1.9*10^-4*delta)-2*w.^(-1)*exp(-9.41*10^-3*delta));
Phi=0.388*10^(-8)*epsilon^(-1/3)*K^(-11/3)*(1+2.35*(K*etta)^(2/3))*ff;
F=K*(W1+W2)*Phi;
F = K*F; % Maybe superfluous and already done in the line before (F=K*(W1+W2)*Phi;) ; I refer here to the line q2 = int(K*q1, K, 0, 1) in your old code
fun = matlabFunction(F,'Vars',[theta,K,z]);
r=integral3(fun,0,2*pi,0,1,0,L)
r = 3.9950e-14 - 8.5054e-08i
m=4*pi*real(r)
m = 5.0203e-13
  댓글 수: 5
Torsten
Torsten 2022년 6월 13일
Removed the line
F = K*F;
But still the integration is unsuccessful.
I didn't want you to remove this line. I only wanted you to check whether the multiplication with K is correct.
Since you integrate over theta, maybe you have to take care of the functional determinant of a coordinate transformation.
Athira T Das
Athira T Das 2022년 6월 15일
@Torsten thanks for your help.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

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

Community Treasure Hunt

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

Start Hunting!

Translated by