## Solving Bessel's Equation numerically

| categories: odes | View Comments

Solving Bessel's Equation numerically

# 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;