Working with Special Functions in .NET

Painting QR Code JIS X 0510 in .NET Working with Special Functions

Figure 11-3

A plot of J1(x)

These characteristics apply in general to Bessel functions of the first kind Let s generate a plot that compares J0(x), J1(x), and J2(x):

>> x = [0:01:20]; u = besselj(0,x); v = besselj(1,x); w = besselj(2,x); >> plot(x,u,x,v,'--',x,w,'-'),xlabel('x'),ylabel('BesselJ(n,x)'), grid on, legend('bessel0(x)','bessel1(x)','bessel2(x)')

Looking at the plot, shown in Figure 11-4, it appears that the oscillations become smaller as n gets larger Bessel functions of the first kind can also be defined for negative integers They are related to Bessel functions of the first kind for positive integers via: J n ( x ) = ( 1)n J n ( x ) Let s do some symbolic computation to find the derivative of a Bessel function:

>> syms n x y >> diff(besselj(n,x)) ans = besselj(n+1,x)+n/x*besselj(n,x)

Figure 11-4

A plot of J0(x), J1(x), and J2(x)

When integrals involving Bessel functions are encountered, MATLAB can also be used For example:

>> int(x^n*besselj(n 1,x)) ans = x^n*besselj(n,x)

That is, modulo a constant of integration:

x n J n 1 ( x )dx = x n J n ( x )

Let s compute the first five terms of the Taylor series expansions of J0(x) and J1(x):

MATLAB also has the other Bessel functions built in Bessel functions of the second kind are implemented with bessely(n, x) In Figure 11-5, we show a plot of the first three Bessel functions of the second kind generated with MATLAB We can also implement another type of Bessel function in MATLAB, the Hankel function These can be utilized by calling besselh(nu, k, z) There are two types of Hankel functions (first kind and second kind), where the kind is denoted by k in MATLAB If we leave the k out of the argument and call besselh(nu, z) MATLAB defaults to the Hankel function of the first kind

Figure 11-5

EXAMPLE 11-2 A cylindrical traveling wave can be shown to be described by:

where we take the real part Implement the real part of this function in MATLAB and plot as a function of r (in meters) for t = 0, p /2w, 3p /4w Assume that k = 500 cm 1 over 0 r 10 m SOLUTION 11-2 First we convert k into the proper units: 100 cm = 50000 m 1 k = 500 cm 1 m Let s generate some plots to see how k impacts the behavior of the real part of the Hankel function First we generate an array for the range of values over which we want to plot:

>> r = linspace(0,10);

Let s suppose that k were unity In this case we could define the Hankel function as:

>> u = besselh(0,r);

You can evaluate the Hankel function at different points to see that it s complex For example:

This indicates that for the function in question to be physical, we have to take the real part, a common practice when studying traveling electromagnetic waves, for example So we plot the real part of this function:

>> plot(r,real(u)),xlabel('r'),ylabel('Hankel(0,r)')

Figure 11-6

The result, shown in Figure 11-6, indicates that this function is a decaying oscillatory function Not surprising since this is a type of Bessel function and we have already seen that Bessel functions are decaying oscillatory functions! Now let s add in the wave number k The code above becomes:

>> >> >> >> r = linspace(0,10); k = 50000; w = besselh(0,k*r); plot(r,real(w)),xlabel('r'),ylabel('Hankel(0,kr)')

Two items are immediately apparent by looking at the graph, shown in Figure 11-7 The first is that the function appears to oscillate quite a bit more rapidly over the distance r Secondly, we notice that the amplitude is much smaller in this case Now let s take a look at the function at different times First considering t = 0, we have:

( (r , 0 ) = H 01) ( kr )

This is what we have in Figure 11-7 Now let s let t = p/2w:

( ( ( (r , /2 ) = H 01) ( kr )e i ( /2 ) = H 01) ( kr )e i /2 = iH 01) ( kr )