Not getting finite result

9 Jun 2022
1 Answer
Answer Accepted
Updated 15 Jun 2022
12 Views (30 days)

0 votes

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 Comments
Show 3 older comments Hide 3 older comments

Athira T Das
Athira T Das on 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 on 10 Jun 2022
Try integrating numerically using "integral3".
Athira T Das
Athira T Das on 13 Jun 2022
Showing error while using integral3, saying too many input arguments.
Torsten
Torsten on 13 Jun 2022
Please include your code.
Athira T Das
Athira T Das on 13 Jun 2022
Ran in:
Ran in:
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)

Sign in to comment.

Sign in to answer this question.

 Accepted Answer

Torsten
Torsten on 13 Jun 2022
Edited: Torsten on 13 Jun 2022
Ran in:

0 votes

Ran in:
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 Comments
Show 3 older comments Hide 3 older comments

Athira T Das
Athira T Das on 13 Jun 2022
Ran in:
Ran in:
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 on 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.
Athira T Das
Athira T Das on 13 Jun 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 on 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 on 15 Jun 2022
@Torsten thanks for your help.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!