Bessel Functions: Difference between revisions
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: | ||
=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: