This post is influenced by a homework problem from a numerical PDE class I took this past summer (2011). This problem examines the DMOLCH routine and Gear's method for a heat diffusion equation. We are to reproduce data found in a journal paper. The numerics to obtain can be found in the journal article "Temperature distribution in dental tissue after interaction with femtosecond laser pulses" by Pike et al in
Applied Optics, Vol. 46, Num. 34, Dec. 2007. I also include some discussion on IMSL, FORTRAN, MOL and many other topics associated with this problem. I will also have individual posts on these topics as well.
IMSL stands for the
International Mathematics and Statistics Library and consist of numerous commercial software library utilized in numerical analysis in various programming languages organized by
Rogue Wave Software [
1,
2]. According to the
2003 Visual Numerics IMSL Fortran Library User’s Guide Math Library Volume 2 of 2 - Mathematical Functions in Fortran [
3], the package routine
DMOLCH is the double precision version of the routine
MOLCH in the
Fortran 77 interface which of course is single precision. In the
Fortran 90 interface the specific interface names are
S MOLCH and
D MOLCH for double and single precision, respectively. The routine
MOLCH according to the 2003 guide,
“Solves a system of partial differential equations of the form \( u_t = f \left( x, t, u, u_x, u_{xx} \right) \) using the method of lines. The solution is represented with cubic Hermite polynomials.”
According to the 2010 guide,
Visual Numerics IMSL Fortran Numerical Library User’s Guide Math Library Version 7.0 [
4], the
MOLCH routine has been “deprecated” by the new routine
MMOLCH, where the definition remains the same as the quote above and the specific F90 interface names for single and double precision are
S MMOLCH and
D MMOLCH, respectivley.
The
method of lines (
MOL) is described by Wouwer et al [
5] in the second section of the first chapter in the book entitled
Adaptive Method of Lines [
6]. The
MOL attacks and aims to solve
partial differential equations (
PDEs) with the form
\begin{equation}
\textbf{u}_t = \textbf{f}(\textbf{u}), \quad \quad \textbf{x}_l < \textbf{x} < \textbf{x}_r, \quad \quad t > 0
\end{equation}
where
\begin{align*}
\textbf{u}_t &= \dfrac{\partial \textbf{u}}{\partial t} \\
\textbf{u} &= \text{vector of dependent variables} \\
t &= \text{initial value independent variable} \\
\textbf{x}& = \text{boundary value independent variables} \\
\textbf{f} &= \text{spatial differential operator} = \textbf{f} \left( \textbf{x}, t, \textbf{u}, \textbf{u}_x,\textbf{u}_{xx}, \ldots \right)
\end{align*}
Wouwer et al [5] divide the
MOL process into two sections. The first step discretizes the spatial derivatives, \( \textbf{u}_x, \textbf{u}_xx, \ldots \) using
finite difference (
FD) or
finite element (
FE) methods. The second step involves integrating the resulting system of semi-discrete
ordinary differential equations (
ODEs) in the initial value variable wrt to time, \( t \).
As an example of the
MOL, Wouwer et al [
5] utilize the linear advection equation to demonstrate the process.
\[
\dfrac{\partial u}{\partial t} = -v \dfrac{\partial u}{\partial x}
\]
To begin, the equation is discretized in the spatial realm of \( x \) into an \( N \) number of grid points which are spaced evenly by the distance \( \delta x \). The spatial derivative, \( \partial u/ \partial x \), is then approximated by its second order centered FD beginning at grid point \( i \) up to the maximum discretization point \( N \) so that
\[
\dfrac{\partial u}{\partial x} = -v \dfrac{u_{i+1} - u_{i-1}}{2 \Delta x} + O \left( \Delta x^2 \right), \quad i = 1, 2, \cdots, N
\]
Then assuming a value of \( v = 1 \) and substituting into \( \partial u/\partial t = -v \partial u/\partial x \) results in the equation
\[
\dfrac{\mathrm{d} u}{\mathrm{d} t} = -\dfrac{u_{i+1} - u_{i-1}}{2 \Delta x}, \quad i = 1, 2, \cdots, N
\]
Wouwer et al [
5] next mention that at this point a system of \( N \) ODEs are generated which can then be integrated by an already developed by a library ODE integrator. An initial condition must then be specified next at the \( N \) grid points and the false points \( u_0 \) and \( u_{N + 1} \) for \( i = 1 \) and \( i = N \). Thus, if an example of the equation is solved for two step functions separated by a time interval of \( t = 50 \) or simply a square pulse, with a spatial length of \( N =101 \) or 100 intervals with a length of \( \Delta x \), the solution exhibits oscillatory behavior as seen in the figure by Wouwer et al [5]. This oscillation is of course false and is caused by the numerics of the problem known as a numerical distortion.
Wouwer et al [5] further discuss some issues here about this problem that I feel like are good points to intake. First is a continuation of the numerical oscillation problem. Wouwer et al express the problem that the oscillations do not dissipate with an increase in grid resolution in the spatial direction \( x \) or, in other words, and increase in the number of intervals, even though the FD equation is of the second order so you would expect as \( \Delta x \) decreased the error would decrease \( O \left( \Delta x^2 \right) \). Thus, the observation is made that will a reasonable method of approximation of a linear PDE broken down into a system of ODEs is accomplished the solutions are inaccurate.
In fact, Wouwer et al [5] explain that the accuracy of the numerical solution is dependent upon the smoothness of the analytical solution. For instance, rather than the two discontinuous jumps encountered with the square wave, a “more smooth” function would be a triangular pulse where \( u\left(x,t\right) \) is continuous but \( \partial u/\partial x \) is not continuous. As shown in the figure the MOL numerical solution (dots) matches much closer with the actual solution (solid line). Finally, if the initial condition is even “smoother” such as a cosine function, then the oscillations diminish even more. Thus, if the initial condition contains non-smooth characteristics such as discontinuities, then the discontinuities spread through the time and space variables which is indicative of hyperbolic PDEs.
However, Gear’s method is also utilized by Pike et al [8]. Gear’s method as defined by Hoffman [7] was developed by Gear in 1971 [9]. Gear’s method contains much larger stability limits for a series of implicit FD equations (FDEs). The formula appears as
\[ y_{n+1} = \gamma \left[ \beta h f_{n+1} + \left( \alpha_0 y_n + \alpha_{-1} y_{n-1} + \alpha_{-2} y_{n-2} + \cdots \right)\right] \]
where k is the global order of the FDE and the coefficients \( \gamma \), \( \beta \), \( \alpha_i \left( i = 0, 1, 2, \ldots \right) \) are shown in the table from Hoffman [7].
In addition, an interesting comment from Hoffman [7] cites that the Gear FDEs are in a FORTRAN package named LSODE which came about from researchers at Lawrence Livermore National Laboratory. According to Hoffman, this package utilizes
"...an elaborate error control procedure based on using various-order Gear FDEs in conjunction with step size halving and doubling."
Additional notes from Hoffman explain that stiff ODEs contain extra issues for solutions. Usually explicit methods do not work for stiff ODEs. However, implicit methods such as the Gear method do work for stiff problems. Hoffman defines stiff ODEs as where
"...at least one eigenvalue has a large negative real part which causes the corresponding component of the solution to vary rapidly compared to the typical scale of variation displayed by the rest of the solution."
Velten [10] states that stiff ODEs pose issues since the numerical solutions may oscillate as demonstrated beforehand. Thus, small step sizes must be implemented to obtain a correct solution. Also of note is Velten’s mention of a few software packages such as Maxima and R. Velten gives another name for the Gear method as the backward differentiation formulas (BDF) method. A part of R’s odesolve package is lsoda which stands for Livermore Solver for Ordinary Differential Equations (automatic) which came about from Petzold and Hindmarsh at California’s Lawrence Livermore National Laboratory [11, 12]. The lsoda package has a unique feature in which it automatically switches between methods for stiff and nonstiff sets of ODEs which is one of the main differences from lsode. As a further side note, according to Soetaert et al [13] the package odesolve is planned to be deprecated and deSolve is the package successor. The nonstiff case utilizes the Adams method with a variable stepsize and a variable order of up to 12th order. For the stiff case of ODEs, the lsoda package utilizes the BDF method with a variable stepsize and a variable order of up to 5th order.
In order to answer the question of stiff ODEs, let's look at an example in A First Course in the Numerical Analysis of Differential Equations by Iserles [14].
Hoffman [7] states that stiffness appears for single linear and nonlinear ODEs, higher-order linear and nonlinear ODEs, and systems of linear and nonlinear ODEs. He lists several definitions for stiffness:
1. An ODE is stiff if the step size required for stability is much smaller than the step size required for accuracy.
2. An ODE is stiff if it contains some components of the solution that decay rapidly compared to other components of the solution.
3. A system of ODEs is stiff if at least one eigenvalue of the system is negative and large compared to the other eigen values of the system.
4. From a practical point of view, an ODE is stiff if the step sixe based on cost (i.e., computational time) is too large to obtain an accurate (i.e., stable) solution.
Hoffman [
7] also gives a couple of examples.
Also, since it is mentioned let’s talk about
Adams methods. According to Iserles [
14],
In progress...to be continued.
References:
[
1] IMSL Numerical Libraries. 07-23-2011, http://en.wikipedia.org/wiki/IMSL Numerical Libraries
[
2] IMSL R Numerical Libraries. imsl-numerical-libraries.aspx 07-23-2011, http://www.roguewave.com/products/
[
3] Visual Numerics IMSL Fortran Library User’s Guide Math Library Volume 2 of 2 - Mathematical Functions in Fortran. 07-23-2011, http://www.absoft.com/Support/Documentation/MathV2.pdf, 2003
[
4] Visual Numerics IMSL Fortran Numerical Library User’s Guide Math Library Version 7.0. 07-23-2011, http://www.roguewave.com/portals/0/products/imsl-numerical-libraries/fortran-library/docs/7.0/math/math.pdf, 2010
[
5] A. Vande Wouwer, P. Saucez, & W. E. Schiesser. Chapter 1 Introduction. A. Vande Wouwer, P. Saucez, & W. E. Schiesser (editors). Adaptive Method of Lines. Chapman & Hall/CRC, Boca Raton, FL. 2001
[
6] A. Vande Wouwer, P. Saucez, & W. E. Schiesser (editors). Adaptive Method of Lines. Chapman & Hall/CRC, Boca Raton, FL. 2001
[
7] J. D. Hoffman. Numerical Methods for Engineers and Scientists. 2 edition. Marcel Dekker, Inc., New York, NY. 2001
[
8] P. Pike, C. Parigger, R. Splinter, and P. Lockhart. Temperature distribution in dental tissue after interaction with femtosecond laser pulses. Applied Optics, Volume 46, Number 34, December 2007
[
9] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, New Jersey. 1971
[
10] K. Velten. Mathematical Modeling and Simulation: Introduction for Scientists and Engineers, Wiley-VCH, Weinheim, Germany. 2009
[
11] A. C. Hindmarsh. ODEPACK, a systematized collection of ODE solvers, Scientific Computing, North-Holland, Amsterdam, The Netherlands, pp. 55–64. 1983
[
12] L. R. Petzold.
Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM Journal on Scientific and Statistical Computing, 4, 136–48. 1983
[
13] K. Soetaert, T. Petzoldt, and R. W. Setzer.
Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33, 4, February 2010
[
14] A. Iserles. A First Course in the Numerical Analysis of Differential Equations. (2nd ed.).
Cambridge University Press, Cambridge, UK. 2009
Additional keywords:
finite difference equations (FDEs), Gear's method, LSODE, LSODA, Adams method, deSolve, R, Maxima, backward differentiation formulas (BDF) method, International Mathematics and Statistics Library (IMSL), Rogue Wave Software, DMOLCH, MOLCH, Fortran, Fortran 90, MMOLCH, Fortran 77, S MOLCH, D MOLCH, S MMOLCH, D MMOLCH, Method of Lines (MOL), Partial Differential Equations (PDEs), Finite Difference (FD), Finite Element (FE), ordinary differential equations (ODEs)
