How to prevent warning when running fplot of a function containing bessel functions?

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)

Answer 1

[Torsten]

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))

[Matthew Niermann]

Thank you this worked perfectly!