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:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \nabla^2 u = \phi^2 u }

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

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \rho} :

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \dfrac{1}{\rho} \dfrac{d}{d \rho} \left( \rho \dfrac{du}{d \rho} \right) = \phi^2 u }

over the region Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle 0 \leq \rho \leq 1} , and the boundary conditions:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \dfrac{du}{d \rho} = 0 }

at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \rho=0} , and also

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

at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \rho=1} .

The Solution

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

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle u(\rho) = \dfrac{ I_0 (\phi \rho) }{ I_0(\phi) + \frac{\phi}{v} I_1(\phi) } }

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle I_0} and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle I_1} are modified Bessel functions of the first kind.

Bessel Functions with Scipy

The Scipy documentation shows you how to call the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle 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: