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 $x^2 y'' + x y' + (x^2 - \nu^2)y=0$ 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 $z = y'$ and $z' = y''$ and performing the change of variables:

$$ y' = z$$

$$ z' = \frac{1}{x^2}(-x z - (x^2 - \nu^2) y$$

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

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 $1/x^2$ term, the ODEs are not defined at x=0.

ans =


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

ans =


xspan = [1e-15 10];
[x,Y] = ode45(@fbessel, xspan, Y0);

hold on
legend 'Numerical solution' 'analytical J_0(x)'

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

ans =


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;
blog comments powered by Disqus