Bessel Functions

The Problem Setup

I had to deal with some Bessel Functions in dealing with an analytical solution of steady state reaction-diffusion partial differential equation. The diffusion-reaction equation is a simplified, non-dimensionalized equation for a single species and single reaction, and is given by:

$\nabla ^{2}u=\phi ^{2}u$

with the following boundary condition on the boundary of the domain:

$v(1-u_{s})={\dfrac {\partial u}{\partial v}}$

For an infinite cylinder (analytical solutions on cylinders typically end up using Bessel Functions), the solution can be written in terms of the non-dimensionalized radial dimension $\rho$:

${\dfrac {1}{\rho }}{\dfrac {d}{d\rho }}\left(\rho {\dfrac {du}{d\rho }}\right)=\phi ^{2}u$

over the region $0\leq \rho \leq 1$, and the boundary conditions:

${\dfrac {du}{d\rho }}=0$

at $\rho =0$, and also

$v(1-u)={\dfrac {\partial u}{\partial \rho }}$

at $\rho =1$.

The Solution

The solution to the above equations is given on an infinite cylinder by the equation:

$u(\rho )={\dfrac {I_{0}(\phi \rho )}{I_{0}(\phi )+{\frac {\phi }{v}}I_{1}(\phi )}}$

where $I_{0}$ and $I_{1}$ are modified Bessel functions of the first kind.

Bessel Functions with Scipy

The Scipy documentation shows you how to call the $I_{v}$ modified Bessel functions:

• http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.iv.html

Implementing the solution given above in Python is dead simple:

import matplotlib.pylab as plt
from scipy.special import *

# Radial dimension
r = linspace(0,1,100)

# Thiele modulus (ratio of reaction versus diffusion driving forces)
phi = 0.1

# Sherwood number (rate of transfer of material across boundary)
v = 100.0

# Create parametric solution vector
for phi in range(10):
    y = [iv(0,phi*rho) / ( iv(0,phi) + (phi/v)*iv(1,phi) ) for rho in r]
    plt.plot(r,y)
    plt.xlabel('r')
    plt.ylabel('u')
    plt.hold(True)

which results in a plot of solutions for various Thiele modulii, for a given Sherwood number, shown below: