From charlesreid1

(Created page with "Dealing with Bessel Functions I had to deal with some Bessel Functions in dealing with an analytical solution of steady state reaction-diffusion partial differential equation. T...")
 
No edit summary
Line 1: Line 1:
Dealing with 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:
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:
Line 32: Line 32:


at <math>\rho=1</math>.
at <math>\rho=1</math>.
=The Solution=
The solution to the above equations is given on an infinite cylinder by the equation:
<math>
u(\rho) = \dfrac{ I_0 (\phi \rho) }{ I_0(\phi) + \frac{\phi}{v} I_1(\phi) }
</math>
where <math>I_0</math> and <math>I_1</math> are modified Bessel functions of the first kind.
=Bessel Functions with Scipy=
The Scipy documentation shows you how to call the <math>I_v</math> 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:
<source lang="python">
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 solution vector
y1 = [iv(0,phi*rho) / ( iv(0,phi) + (phi/v)*iv(1,phi) ) for rho in r]
plot(r,y1)
# 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.set_xlabel('r')
    plt.set_ylabel('u')
    plt.hold(True)
</source>
which results in a plot of solutions for various Thiele modulii, for a given Sherwood number, shown below:
[[Category:Mathematics]]
[[Category:Python]]

Revision as of 19:53, 25 February 2014

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:

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 solution vector
y1 = [iv(0,phi*rho) / ( iv(0,phi) + (phi/v)*iv(1,phi) ) for rho in r]
plot(r,y1)

# 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.set_xlabel('r')
    plt.set_ylabel('u')
    plt.hold(True)

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