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 comentarios

Athira T Das
Athira T Das el 10 de Jun. de 2022
No error is there, but not getting a finite value after doing triple integration. Output is just the function with integral sign. I need an exact value of the integral.
Torsten
Torsten el 10 de Jun. de 2022
Try integrating numerically using "integral3".
Athira T Das
Athira T Das el 13 de Jun. de 2022
Showing error while using integral3, saying too many input arguments.
Torsten
Torsten el 13 de Jun. de 2022
Please include your code.
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)

Iniciar sesión para comentar.

 Respuesta aceptada

Torsten
Torsten el 13 de Jun. de 2022
Editada: Torsten el 13 de Jun. de 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 comentarios

If I change the limit from 0 to 1 to 0 to Inf, then the integration is unsuccessful due to the presence of Infinity.
But my function is defined at Infinity.
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=-0.08
w = -0.0800
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 = 
F=K*(W1+W2)*Phi;
F = K*F;
fun = matlabFunction(F,'Vars',[theta,K,z]);
r=integral3(fun,0,2*pi,0,Inf,0,L)
Warning: Inf or NaN value encountered.
Warning: The integration was unsuccessful.
r = NaN
m=4*pi*real(r)
m = NaN
Torsten
Torsten el 13 de Jun. de 2022
"fun" is the function you have defined, and integral3 is a reliable integrator.
What shall I say ? Time to debug your variable definitions.
Or maybe the line
F = K*F;
is superfluous because the multiplication by K is already done in
F=K*(W1+W2)*Phi;
And please enter a semicolon at the end of MATLAB lines in order not to output every temporary result.
Athira T Das
Athira T Das el 13 de Jun. de 2022
Removed the line
F = K*F;
But still the integration is unsuccessful
Change limit from Infinity to 10^3 then
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.
Torsten
Torsten el 13 de Jun. de 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 el 15 de Jun. de 2022
@Torsten thanks for your help.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Centro de ayuda y File Exchange.

Preguntada:

el 9 de Jun. de 2022

Comentada:

el 15 de Jun. de 2022

Community Treasure Hunt

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

Start Hunting!

Translated by