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 Kommentare

Athira T Das
Athira T Das am 10 Jun. 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 am 10 Jun. 2022
Try integrating numerically using "integral3".
Athira T Das
Athira T Das am 13 Jun. 2022
Showing error while using integral3, saying too many input arguments.
Torsten
Torsten am 13 Jun. 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)

Melden Sie sich an, um zu kommentieren.

 Akzeptierte Antwort

Torsten
Torsten am 13 Jun. 2022
Bearbeitet: Torsten am 13 Jun. 2022

0 Stimmen

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 Kommentare

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 am 13 Jun. 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.
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 am 13 Jun. 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 am 15 Jun. 2022
@Torsten thanks for your help.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Hilfe-Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by