Not getting finite result

2 visualizzazioni (ultimi 30 giorni)
Athira T Das
Athira T Das il 9 Giu 2022
Commentato: Athira T Das il 15 Giu 2022
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 Commenti
Torsten
Torsten il 13 Giu 2022
Please include your code.
Athira T Das
Athira T Das il 13 Giu 2022
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)

Accedi per commentare.

Risposta accettata

Torsten
Torsten il 13 Giu 2022
Modificato: Torsten il 13 Giu 2022
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 Commenti
Torsten
Torsten il 13 Giu 2022
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 il 15 Giu 2022
@Torsten thanks for your help.

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by