How to prevent warning when running fplot of a function containing bessel functions?
1 view (last 30 days)
Show older comments
Hi everyone!
I am worried that this warning is omitting significant values in my code. I would be grateful if someone could tell me what is wrong with my code.
My code:
clearvars
%Sets respective variable values
t= 0.335.*10^-9;
ro=5;
ro1 = 0.26.*10^-6;
ro2 = 0.18.*10^-6;
z1=5;
ks = 10;
g = 3;
syms z r
qo=6;
%Sets definition for Tp(z) function
D(z) = (-besseli(0,z) .* besselk(1,z))-(besselk(0,z) .* besseli(1,z))
fplot(D(z),[0.1,2])
figure()
B(z) = (besselk(0,z) .* (qo/g) .* exp(-(z.^2)/(z1^2)) ) ./ (D(z))
pretty(B)
O(z) = int(B(z),z)
M(z) = real(O(z))
d1=digits(50);
fplot(vpa(M(z)),[0.01,15])
I recieve this regular output:
D(z) =
- besseli(0, z)*besselk(1, z) - besseli(1, z)*besselk(0, z)
B(z) =
-(2*exp(-z^2/25)*besselk(0, z))/(besseli(0, z)*besselk(1, z) + besseli(1, z)*besselk(0, z))
/ 2 \
| z |
exp| - -- | besselk(0, z) 2
\ 25 /
- ---------------------------------------------------------
besseli(0, z) besselk(1, z) + besseli(1, z) besselk(0, z)
O(z) =
-int((2*exp(-z^2/25)*besselk(0, z))/(besseli(0, z)*besselk(1, z) + besseli(1, z)*besselk(0, z)), z)
M(z) =
-int(2*real((exp(-z^2/25)*besselk(0, z))/(besseli(0, z)*besselk(1, z) + besseli(1, z)*besselk(0, z))), z)
Then this warning twice:
Warning: Inf or NaN value encountered.
> In integralCalc/midpArea (line 397)
In integralCalc (line 65)
In integral (line 87)
In symengine>@(z)integral(@(z)real((exp(z.^2.*-4.0e-2).*besselk(0.0,z))./(besseli(1.0,z).*besselk(0.0,z)+besselk(1.0,z).*besseli(0.0,z))).*2.0,0,z).*-1.0
In matlab.graphics.function.internal.sym2fn>tryBoth (line 73)
In matlab.graphics.function.internal.sym2fn>@(varargin)tryBoth(f1,f2,varargin{:}) (line 42)
In matlab.graphics.function.internal.sym2fn>@(varargin)arrayfun(fn,varargin{:}) (line 66)
In matlab.graphics.function.FunctionLine>getFunction
In matlab.graphics.function/FunctionLine/updateFunction
In matlab.graphics.function/FunctionLine/set.Function_I
In matlab.graphics.function/FunctionLine/set.Function
In matlab.graphics.function.FunctionLine
In fplot>singleFplot (line 245)
In fplot>@(f)singleFplot(cax,{f},limits,extraOpts,args) (line 200)
In fplot>vectorizeFplot (line 200)
In fplot (line 166)
In ResearchCode (line 27)
0 Comments
Accepted Answer
Torsten
on 19 Jul 2022
Edited: Torsten
on 19 Jul 2022
O is only determined up to a constant. You must fix this constant of integration to get reasonable results (I assumed O(1e-8) = 0).
Here is a numerical solution:
%Sets respective variable values
t= 0.335.*10^-9;
ro=5;
ro1 = 0.26.*10^-6;
ro2 = 0.18.*10^-6;
z1=5;
ks = 10;
g = 3;
%syms z r
qo=6;
%Sets definition for Tp(z) function
D = @(z)(-besseli(0,z) .* besselk(1,z))-(besselk(0,z) .* besseli(1,z));
z = 0.1:0.1:2;
figure(1)
plot(z,D(z))
B = @(z) (besselk(0,z) .* (qo/g) .* exp(-(z.^2)/(z1^2)) ) ./ D(z);
O = @(z) integral(@(u)B(u),1e-8,z);
z = 1e-8:0.01:15+1e-8;
plot(z,arrayfun(@(z)O(z),z))
More Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!