*** Solving Bessel's Equation numerically
:PROPERTIES:
:categories: ODE, math
:date: 2013/02/07 09:00:00
:updated: 2013/03/06 18:33:34
:END:
[[http://matlab.cheme.cmu.edu/2011/08/08/solving-bessels-equation-numerically/][Matlab post]]
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$.
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. If we start very close to zero instead, we avoid the problem.
#+BEGIN_SRC python
import numpy as np
from scipy.integrate import odeint
from scipy.special import jn # bessel function
import matplotlib.pyplot as plt
def fbessel(Y, x):
nu = 0.0
y = Y[0]
z = Y[1]
dydx = z
dzdx = 1.0 / x**2 * (-x * z - (x**2 - nu**2) * y)
return [dydx, dzdx]
x0 = 1e-15
y0 = 1
z0 = 0
Y0 = [y0, z0]
xspan = np.linspace(1e-15, 10)
sol = odeint(fbessel, Y0, xspan)
plt.plot(xspan, sol[:,0], label='numerical soln')
plt.plot(xspan, jn(0, xspan), 'r--', label='Bessel')
plt.legend()
plt.savefig('images/bessel.png')
#+END_SRC
#+RESULTS:
[[./images/bessel.png]]
You can see the numerical and analytical solutions overlap, indicating they are at least visually the same.