[Scilab-users] Bessel functions

[Claus Futtrup]

Hi

I've made a small script to play around with Bessel functions of the 
first kind... but this is very basic and I'm wondering if there's a 
smarter way.

Below is the script I made (the function + a small test which plots the 
result).

Best regards,
Claus

// bessel_test.sce

function  z=Jn(n,x)
     // The following power series approximates the nth Bessel
     // function of the first kind for each input x
     // In acoustics x = 2ka, defines the iput frequency k = omega / c and size
     // of the piston radiator
     powerseries  =  0;
     for  m=0:19  // actually it should be infinity, but 10 approximates OK ...
         powerseries_m  =  ((-1)^m  /  (factorial(m)  *  factorial(m  +  n)))  *  (x/2)^(2*m);
         powerseries  =  powerseries  +  powerseries_m;  // sum the powerseries
     end
     z  =  ((x/2)^n)  .*  powerseries;
endfunction

x_array  =  0:0.1:9.9;  // define 100 points on the x-axis
                      // with 2ka (x) from 0 to 10, and with m = 0-19,
                      // this approximation is reasonably good for 2ka<  10

z_0_array  =  Jn(0,x_array);  // Calculate J0

z_1_array  =  Jn(1,x_array);  // Calculate J1

scf();
a  =  gca();
plot(x_array,z_0_array,'-b');
plot(x_array,z_1_array,'-r');
xtitle("Bessel functions","x (2ka)","output (z)");
legend("J0","J1");


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20150124/1566f2ce/attachment.htm>