*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

*and consist of numerous commercial software library utilized in numerical analysis in various programming languages organized by*

**I**nternational**M**athematics and**S**tatistics**L**ibrary**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 speciﬁc interface names are

**S**

**MOLCH**and

**D MOLCH**for double and single precision, respectively. The routine

**MOLCH**according to the 2003 guide,

“According to the 2010 guide,Solves a system of partial diﬀerential 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.”

*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 deﬁnition remains the same as the quote above and the speciﬁc 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 ﬁrst chapter in the book entitled

*Adaptive Method of Lines*[

**6**]. The

**MOL**attacks and aims to solve

__(__

**p**artial**d**ifferential**e**quations**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 ﬁrst step discretizes the spatial derivatives, \( \textbf{u}_x, \textbf{u}_xx, \ldots \) using

**ﬁnite diﬀerence**(

**FD**) or

**ﬁnite element**(

**FE**) methods. The second step involves integrating the resulting system of semi-discrete

__(__

**o**rdinary**d**ifferential**e**quations**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 speciﬁed 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 deﬁned by Hoﬀman [**7**] was developed by Gear in 1971 [**9**]. Gear’s method contains much larger stability limits for a series of implicit**FD equations**(**FDE**s). 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 coeﬃcients \( \gamma \), \( \beta \), \( \alpha_i \left( i = 0, 1, 2, \ldots \right) \) are shown in the table from Hoﬀman [

**7**].
In addition, an interesting comment from Hoﬀman [

**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 Hoﬀman, 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 Hoﬀman explain that stiﬀ ODEs contain extra issues for solutions. Usually explicit methods do not work for stiﬀ ODEs. However, implicit methods such as the Gear method do work for stiﬀ problems. Hoﬀman deﬁnes stiﬀ 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 stiﬀ 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 diﬀerentiation formulas**(**BDF**)**method**. A part of*R*’s**odesolve**package is**lsoda**which stands for**L**ivermore**S**olver for**O**rdinary**D**iﬀerential Equations (**a**utomatic) 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 stiﬀ and nonstiﬀ sets of ODEs which is one of the main diﬀerences 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 nonstiﬀ case utilizes the**Adams method**with a variable stepsize and a variable order of up to 12th order. For the stiﬀ 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

Hoffman [

*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 deﬁnitions for stiﬀness:**1.**An ODE is stiﬀ if the step size required for stability is much smaller than the step size required for accuracy.

**2.**An ODE is stiﬀ if it contains some components of the solution that decay rapidly compared to other components of the solution.

**3.**A system of ODEs is stiﬀ 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 stiﬀ if the step sixe based on cost (i.e., computational time) is too large to obtain an accurate (i.e., stable) solution.

Hoﬀman [

**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. Hoﬀman. 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 Diﬀerential Equations, Prentice-Hall, Englewood Cliﬀs, 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 stiﬀ and nonstiﬀ systems of ordinary diﬀerential equations*. SIAM Journal on Scientiﬁc 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 Diﬀerential 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 Diﬀerence (FD), Finite Element (FE), ordinary differential equations (ODEs)

## No comments:

## Post a Comment