## Solving Bessel's Equation numerically

August 08, 2011 at 01:53 PM | categories: odes | View Comments

# Solving Bessel's Equation numerically

Reference Ch 5.5 Kreysig, Advanced Engineering Mathematics, 9th ed.

Bessel's equation comes up often in engineering problems such as heat transfer. The solutions to this equation are the Bessel functions. To solve this equation numerically, we must convert it to a system of first order ODEs. This can be done by letting and and performing the change of variables:

if we take the case where , the solution is known to be the Bessel function , which is represented in Matlab as besselj(0,x). The initial conditions for this problem are: and .

```
function bessel_ode
```

clear all; close all; clc; % define the initial conditions y0 = 1; z0 = 0; Y0 = [y0 z0];

There is a problem with our system of ODEs at x=0. Because of the term, the ODEs are not defined at x=0.

fbessel(0,Y0)

ans = 0 NaN

if we start *very* close to zero instead, we avoid the problem

fbessel(1e-15,Y0)

ans = 0 -1

xspan = [1e-15 10]; [x,Y] = ode45(@fbessel, xspan, Y0); figure hold on plot(x,Y(:,1),'b-'); plot(x,besselj(0,x),'r--') legend 'Numerical solution' 'analytical J_0(x)'

You can see the numerical and analytical solutions overlap, indicating they are the same.

```
'done'
```

ans = done

function dYdx = fbessel(x,Y) % definition of the Bessel equation as a system of first order ODEs nu = 0; y = Y(1); z = Y(2); dydx = z; dzdx = 1/x^2*(-x*z - (x^2 - nu^2)*y); dYdx = [dydx; dzdx]; % categories: ODEs % tags: math % post_id = 763; %delete this line to force new post;