## Sunday, September 4, 2011

### Numerical Mathematics - Solving the heat diffusion equation

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 speciﬁc 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 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.
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 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 partial differential equations (PDEs) with the form $$\textbf{u}_t = \textbf{f}(\textbf{u}), \quad \quad \textbf{x}_l < \textbf{x} < \textbf{x}_r, \quad \quad t > 0$$ 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 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 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 (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 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 Livermore Solver for Ordinary Diﬀerential 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 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 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