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: