besseli - Modified Bessel function of the first kind (I sub alpha).
besselj - Bessel function of the first kind (J sub alpha).
besselk - Modified Bessel function of the second kind (K sub alpha).
bessely - Bessel function of the second kind (Y sub alpha).
besselh - Modified Bessel function of the third kind.
If alpha and x are arrays of the same size, the result y is also that size. If either input is a scalar, it is expanded to the other input's size. If one input is a row vector and the other is a column vector, the result y is a two-dimensional table of function values.
I_alpha and K_alpha are 2 independant solutions of the modified Bessel's differential equation :
2 2 2
x y(x)" + x y(x)' - (x + alpha ) y(x) = 0 ,
J_alpha and Y_alpha are 2 independant solutions of the Bessel's differential equation :
2 2 2
x y(x)" + x y(x)' + (x - alpha ) y(x) = 0 ,
// besselI functions
// ==================
x = linspace(0.01,10,5000)';
xbasc()
subplot(2,1,1)
plot2d(x,besseli(0:4,x), style=2:6)
legend('I'+string(0:4),2);
xtitle("Some modified bessel functions of the first kind")
subplot(2,1,2)
plot2d(x,besseli(0:4,x,1), style=2:6)
legend('I'+string(0:4),1);
xtitle("Some modified scaled bessel functions of the first kind")
// besselJ functions
// =================
x = linspace(0,40,5000)';
xbasc()
plot2d(x,besselj(0:4,x), style=2:6, leg="J0@J1@J2@J3@J4")
legend('I'+string(0:4),1);
xtitle("Some bessel functions of the first kind")
// use the fact that J_(1/2)(x) = sqrt(2/(x pi)) sin(x)
// to compare the algorithm of besselj(0.5,x) with a more direct formula
x = linspace(0.1,40,5000)';
y1 = besselj(0.5, x);
y2 = sqrt(2 ./(%pi*x)).*sin(x);
er = abs((y1-y2)./y2);
ind = find(er > 0 & y2 ~= 0);
xbasc()
subplot(2,1,1)
plot2d(x,y1,style=2)
xtitle("besselj(0.5,x)")
subplot(2,1,2)
plot2d(x(ind), er(ind), style=2, logflag="nl")
xtitle("relative error between 2 formulae for besselj(0.5,x)")
// besselK functions
// =================
x = linspace(0.01,10,5000)';
xbasc()
subplot(2,1,1)
plot2d(x,besselk(0:4,x), style=2:6, rect=[0,0,6,10])
legend('K'+string(0:4),1);
xtitle("Some modified bessel functions of the second kind")
subplot(2,1,2)
plot2d(x,besselk(0:4,x,1), style=2:6, rect=[0,0,6,10])
legend('K'+string(0:4),1);
xtitle("Some modified scaled bessel functions of the second kind")
// besselY functions
// =================
x = linspace(0.1,40,5000)'; // Y Bessel functions are unbounded for x -> 0+
xbasc()
plot2d(x,bessely(0:4,x), style=2:6, rect=[0,-1.5,40,0.6])
legend('Y'+string(0:4),4);
xtitle("Some bessel functions of the second kind")
// besselH functions
// =================
x=-4:0.025:2; y=-1.5:0.025:1.5;
X=x(:)*ones(y);Y=ones(x')*y;
H = besselh(0,1,X+%i*Y);
clf();f=gcf();
xset("fpf"," ")
f.color_map=jetcolormap(17);
contour2d(x,y,abs(H),0:0.2:3.2,rect=[-4,-1.5,2,1.5])
Slatec : dbesi.f, zbesi.f, dbesj.f, zbesj.f, dbesk.f, zbesk.f, dbesy.f, zbesy.f, zbesh.f
Drivers to extend definition area (Serge Steer INRIA): dbesig.f, zbesig.f, dbesjg.f, zbesjg.f, dbeskg.f, zbeskg.f, dbesyg.f, zbesyg.f, zbeshg.f