deal.II Homepage

deal.II Homepage: A Finite Element Differential Equations Analysis Library - October 9th, 2011: deal.II 7.1 released

What is deal.II? 

deal.II is a C++ program library targeted at the computational solution of partial differential equations using adaptive finite elements. It uses state-of-the-art programming techniques to offer you a modern interface to the complex data structures and algorithms required. 

The main aim of deal.II is to enable rapid development of modern finite element codes, using among other aspects adaptive meshes and a wide array of tools classes often used in finite element program. Writing such programs is a non-trivial task, and successful programs tend to become very large and complex. We believe that this is best done using a program library that takes care of the details of grid handling and refinement, handling of degrees of freedom, input of meshes and output of results in graphics formats, and the like. Likewise, support for several space dimensions at once is included in a way such that programs can be written independent of the space dimension without unreasonable penalties on run-time and memory consumption. 

deal.II is widely used in many academic and commercial projects. For its creation, its principal authors have received the 2007 J. H. Wilkinson Prize for Numerical Software. It is also part of the industry standard SPEC CPU 2006 benchmark suite used to determine the speed of computers and compilers, and comes pre-installed on the machines offered by the commercial Sun Grid program. 

deal.II emerged from work at the Numerical Methods Group at Universität Heidelberg, Germany, which is at the forefront of adaptive finite element methods and error estimators. Today, it is maintained by two of its original authors at Texas A&M University, and dozens of contributors and several hundred users are scattered around the world (see our credits page for a detailed list of people contributing to deal.II).

What deal.II can offer you? 

If you are active in the field of adaptive finite element methods, deal.II might be the right library for your projects. Among other features, it offers: 

Support for one, two, and three space dimensions, using a unified interface that allows to write programs almost dimension independent. 

Handling of locally refined grids, including different adaptive refinement strategies based on local error indicators and error estimators. Both h, p, and hp refinement is fully supported for continuous and discontinuous elements. 

Support for a variety of finite elements: Lagrange elements of any order, continuous and discontinuous; Nedelec and Raviart-Thomas elements of any order; elements composed of other elements. 

Parallelization on single machine through the Threading Build Blocks and across nodes via MPI. deal.II has been shown to scale to at least 16k processors. 

Extensive documentation: all documentation is available online in a logical tree structure to allow fast access to the information you need. If printed it comprises more than 500 pages of tutorials, several reports, and presently some 5,000 pages of programming interface documentation with explanations of all classes, functions, and variables. All documentation comes with the library and is available online locally on your computer after installation. 

Modern software techniques that make access to the complex data structures and algorithms as transparent as possible. The use of object oriented programming allows for program structures similar to the structures in mathematical analysis. 

A complete stand-alone linear algebra library including sparse matrices, vectors, Krylov subspace solvers, support for blocked systems, and interface to other packages such as Trilinos, PETSc and METIS.

Support for several output formats, including many common formats for visualization of scientific data.

Portable support for a variety of computer platforms and compilers.
Free source code under an Open Source license, and the invitation to contribute to further development of the library.

Numerical Math - The Finite Element Method (FEM)

From Landau [1], the finite element method (FEM) is explained as solving a PDE where the whole region or domain is subsectioned into smaller areas known as elements. Next an initial solution is formulated for the PDE in each of the elements.  The next step comes in the form of modifying the parameters of the initial formalization known as a best fit.

Pepper and Heinrich [2] open their FEM book by explaining that the method is a numerical formulation that solves physics and engineering models described by differential equations  (DEs). Similar to finite difference methods (FDM), FEM models a geometric region which is then broken up into a number of smaller subregions which results in a network known as a mesh. One difference between FDM and FEM includes mesh type (geometrically). The FDM requires that the network be made of orthogonal rows and columns (squares and rectangles) while the FEM does not require this limitation and can be, in fact, any shape such as triangles and/or quadrilaterals in two dimensions and tetrahedrons and/or hexahedrons in three dimensions.  Next, each FE is initialized with an approximate functions of the unknown variables for which to be solved. Expansions determine the variables approximations and appear as a linear or higher-order polynomials functions. These expansions, in turn, depend upon the geometric shape of the elements and location known as nodes. A second difference between FDM and FEM noted by Pepper and Heinrich [3] entails the solution method. The FEM integrates over each subregion then adds or connects them to make up the whole. This integration and summation engenders a set of finite linear equations per section which can then be solved utilizing linear algebra methods. Jiang [6] also discusses that FEM does not...

The network of elements and nodes or where these elements connect make up discrete systems such as a trusses, circuits, and fluid transport pipes [5]. In order to solve for the system variables such as displacements, electric potentials, and pressures, one can begin with a known and simple parameter such as

Hooke’s Law:

\begin{equation} F = \dfrac{\Delta L}{R} = E A \dfrac{\Delta L}{L} \end{equation}

Ohm’s Law:

\begin{equation} i = \dfrac{\Delta V}{R} = \dfrac{A}{\rho} \dfrac{\Delta V}{L} \end{equation}

Poiseuille’s Law:

\begin{equation} \dot{m} = \dfrac{\Delta p}{R} = \dfrac{\rho \pi D^4}{128 \mu} \dfrac{\Delta p}{L} \end{equation}

Jiang [6] notes that FEM has been utilized as one of the most general numerical techniques to solve DEs and done so with great success. Jiang even quotes Oden as

“Perhaps no other family of approximation methods has had a greater impact on the theory and application of numerical methods during the twentieth century”

In progress...to be continued.

---

References:

[1] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu. The Finite Element Method - It’s Basis and Fundamentals, 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

[2] D. W. Pepper and J. C. Heinrich. The Finite Element Method: Basic Concepts and Applications, Taylor & Francis Hemisphere Publishing Corporation,. 1992

[4] J. C. Heinrich and D. W. Pepper. The Intermediate Finite Element Method: Fluid Flow and Heat Transfer Applications, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1999

[5] G. Comini, S. D. Giudice and C. Nonino. Finite Element Analysis in Heat Transfer: Basic Formulation & Linear Problems, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1994

[6] B.-N. Jiang. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics, Springer-Verlag, Berlin, Germany. 1998

[7] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. The Finite Element Method for Fluid Dynamics, (Volume 3) 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

Numerical Math - Square wire problem - finite difference method

The square-wire problem is described by the 2-D Laplace and Poisson equations [1].

\begin{equation} \dfrac{\partial^2 U \left( x, y \right)}{\partial x^2} + \dfrac{\partial^2 U \left( x, y \right)}{\partial y^2} = \left\{ \begin{array}{cr} 0 & \text{Laplace's Equation}\\ -4 \pi \rho \left( \mathbf{x} \right) & \text{Poisson's Equation} \end{array} \right. \end{equation}

Next, in order to obtain the second order derivatives a Taylor series expansion with steps of positive and negative values in both \( x \) and \( y \) are executed for \( U \).

\begin{equation} U \left( x + \Delta x, y \right) = U \left( x, y \right) + \dfrac{\partial U}{\partial x} \Delta x + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 + \cdots, \end{equation}

\begin{equation} U \left( x - \Delta x, y \right) = U \left( x, y \right) - \dfrac{\partial U}{\partial x} \Delta x + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 - \cdots, \end{equation}

\begin{equation} U \left( x, y + \Delta y \right) = U \left( x, y \right) + \dfrac{\partial U}{\partial y} \Delta y + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 + \cdots, \end{equation}

\begin{equation} U \left( x, y - \Delta y \right) = U \left( x, y \right) - \dfrac{\partial U}{\partial y} \Delta y + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 - \cdots, \end{equation}

Now, adding both equations each for \( x \) and \( y \) the odd terms cancel, and truncating at \( O \left( \Delta^4 \right) \) obtains

\begin{equation} U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) = 2 U \left( x, y \right) + \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 + O \left( \Delta x^4 \right) \end{equation}

\begin{equation} U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) = 2 U \left( x, y \right) + \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 + O \left( \Delta y^4 \right) \end{equation}

Rearranging to segregate the second derivative engenders

\begin{equation}  \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 = U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right) \end{equation}

\begin{equation} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 = U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right) \end{equation}

Finally, dividing produces

\begin{equation} \dfrac{\partial^2 U}{\partial x^2} = \dfrac{U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right)}{\Delta x^2} \end{equation}

\begin{equation} \dfrac{\partial^2 U}{\partial y^2} = \dfrac{U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right)}{\Delta y^2} \end{equation}

Next, substituting the second order derivatives into the Poisson equation produces

\begin{equation} \begin{array}{l} \dfrac{U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right)}{\Delta x^2} \\ \hspace{1.15 in} + \dfrac{U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right)}{\Delta y^2} = -4 \pi \rho \left( x, y \right) \end{array} \end{equation}

If the grid spacing is the same in both the \( x \) and \( y \) coordinates, that is \( \Delta x = \Delta y \), then the equation reduces to

\begin{equation} \begin{array}{l} U \left( x + \Delta, y \right) + U \left( x - \Delta, y \right) + U \left( x, y + \Delta \right) \\ \hspace{2 in} + \hspace{1 mm} U \left( x, y - \Delta \right) - 4 U \left( x, y \right) = -4 \pi \rho \left( x, y \right) \Delta^2 \end{array} \end{equation}

Algebraically, the solution \( U \left( x, y \right) \) can be shown as

\begin{equation} \begin{array}{l} U \left( x, y \right) \simeq \dfrac{1}{4} \left[ U \left( x + \Delta, y \right) + U \left( x - \Delta, y \right) \right. \\ \hspace{2 in} \left. + \hspace{1 mm} U \left( x, y + \Delta \right) + U \left( x, y - \Delta \right) \right] + \pi \rho \left( x, y \right) \Delta^2 \end{array} \end{equation}

Lattice spacings can be identified as \( x = x_0 \pm \Delta \) and \( y = y_0 \pm \Delta \) which correspond to \(i \pm 1 \) and \( j \pm 1 \) where \( i, j = 0, 1, \ldots, N_{max - 1} \).  Thus, the finite difference algorithm appears as

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) + \pi \rho_{i, j}\Delta^2 \end{equation}

However, for simplicity Landau et al just ask to solve for the homogeneous equation

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \end{equation}

A few notes here from Landau et al.  For one a direct solution cannot be obtain from this algorithm.  That is, a process must be repeated several times in order for a solution to converge!  Thus, this is where efficiency and proficiency becomes very important in numerical methods.

First, a guess is imposed then swept over the domain which finds values for all variables at each point or node.  The solution has "converged" when the current generation does not change compared to the previous generation of values within a "tolerance" or error/precision (absolute value of the difference between old values and new) specified by the user or technique.  Another outcome could be that the values emerge as not expected which at that point the user must go back and figure out from where the error is coming.  Landau et al give the definition of relaxed as when the first guess has converged into a solution.

Landau et al pose "Does it always converge, and if so, does it converge fast enough to be useful?"

Another important note is the coordinate system used.  This problem was done in rectangular Cartesian coordinates.  Other problems may work better in other coordinate systems such as cylindrical or spherical depending upon the geometry of the problem in which case the Laplacian expansion \( \Delta^2 \) will be different.

One way to solve the algorithm is to utilize the Jacobi method.  This most basic method simply evaluates each point, sweeping across the domain.  The initialization and BCs symmetry are preserved in this way.  However, in order to speed up this method, a simple modification can be implemented where values that are known from the current run compared to other nodes can be updated immediately.  Thus, reducing the amount of needed storage by half since two nodes will be know for one in this case and the previous sweep will not be necessary to store since information is updated instantaneously when available.  This technique is called the Gauss-Seidel method.  The GS method leads to quicker convergence, less round-off error, and less needed storage.  However, symmetry is not conserved which the user hopes does not affected the solution outcome.

If the sweep begins at the top left corner then the GS algorithm looks like

\begin{equation} U_{i, j} = \dfrac{1}{4} \left[ U^{(old)}_{i + 1, j} + U^{(new)}_{i - 1, j} + U^{(old)}_{i, j + 1} + U^{(new)}_{i, j - 1} \right] \end{equation}

Octave code:

%This code was written by Tim Barber for the numerical PDE class taught by Dr. Parigger in the Summer of 2011

%This file corresponds specifically to HW set #1

%Problem 1.3 Square Wire Problem: Chapter 17 of the Landau-course book

discusses the square-wire problem, viz. elliptic PDE. Section 17.4.2

elaborates on the solution using the finite difference method. Your task

here is to (i) derive the finite difference algorithm Eq. 17.27, (ii)

implement the code using source programing (see the listed source

code), and (iii) display the results using gnuplot (along the lines of

Section 17.4.2).

clf

clc

clear

%input data manually

N_max = 100;

iter = 1000;

%initialize matrix

for i = 1:N_max

for j = 1:N_max

if j = 100

U(i, 100) = 99;

else

U(i, j) = 0;

j = j + 1;

end

end

i = i +1;

end

%iterate matrix using the Gauss-Seidel technique

for n = 1:iter + 1

for i = 2:N_max - 1

for j = 2:N_max - 1

U(i, j) = (1/4)*(U(i + 1, j) + U(i - 1, j)

+ U(i, j + 1) + U(i, j - 1));

j = j + 1;

end

i = i +1;

end

n = n + 1;

end

%plot of U(x, y)

mesh (U)

%label title and axes

title('HW 1.3', "fontsize", 20)

ylabel('y', "fontsize", 16)

xlabel('x', "fontsize", 16)

zlabel('U(x, y)', "fontsize", 16)

%commands to save plot as a .png file

print("HW_1_3_mesh.png")

replot

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics -

a Introductory Computational Science, Princeton University Press, Princeton, New Jersey.

2008

Partial Differential Equations - Types

One very important step for solving, working with, and understanding partial differential equations (PDEs) is knowing the classifications of types.

A general second order PDE can be written as (Özişik 1993, Tannehill et al 1997, Hoffman 2001, Chung 2002, Arfken & Weber 2005, Chung 2010):

\[ A \dfrac{\partial^2 \phi}{\partial x^2} + B \dfrac{\partial^2 \phi}{\partial x \partial y} + C \dfrac{\partial^2 \phi}{\partial y^2} + D \dfrac{\partial \phi}{\partial x}+ E \dfrac{\partial \phi}{\partial y} + F \phi + G \left( x, y \right) = 0 \]

The types of PDEs is based on the coefficients \( A \),  \( B \), and  \( C \) which produces the characteristic equations also known as discriminants:

\[ \begin{array}{lll} B^2 - 4AC < 0 & \text{Elliptic} & \text{Complex characteristic curves} \\ B^2 - 4AC = 0 & \text{Parabolic} & \text{Real and repeated characteristic curves}\\ B^2 - 4AC > 0 & \text{Hyperbolic}  & \text{Real and distinct characteristic curves} \end{array} \]

Some examples of each type of equation are:

Elliptic:

Laplace's Differential Equation:

\[ \nabla^2 \phi = 0 \]

in 2-d Cartesian:

\[ \dfrac{\partial^2 \phi}{\partial x^2} + \dfrac{\partial^2 \phi}{\partial y^2} = 0 \]

Steady Heat Diffusion Equation:


\[ \dfrac{\partial^2 T}{\partial x^2} + \dfrac{\partial^2 T}{\partial y^2} = 0 \]

Poisson's Differential Equation (non-homogeneous Laplace Equation):


\[ \nabla^2 \phi = f \]



in 2-d Cartesian:

\[ \dfrac{\partial^2 \phi}{\partial x^2} + \dfrac{\partial^2 \phi}{\partial y^2} = f \left( x, y \right) \]


Parabolic:

Heat or Unsteady Diffusion Equation:

\[ \dfrac{\partial \phi}{\partial t} = \alpha \nabla^2 \phi = 0 \]

in 1-d Cartesian:

\[ \dfrac{\partial \phi}{\partial t} = \alpha \dfrac{\partial^2 \phi}{\partial x^2} = 0 \]


Hyperbolic:

Wave Equation:

\[ \dfrac{\partial^2 \phi}{\partial t^2} = c \nabla^2 \phi = 0 \]

in 1-d Cartesian:

\[ \dfrac{\partial^2 \phi}{\partial t^2} = c \dfrac{\partial^2 \phi}{\partial x^2} = 0 \]

Physical classification breaks up PDE's for further understanding (Tannehill et al 1997, Hoffman 2001).

Equilibrium Problems:


The first is called an equilibrium problem.  The solution of the PDE must be satisfied in a closed domain where conditions are set on the boundary.  Equilibrium problems are also known as boundary value problems and include examples such as steady-state temperature distributions, incompressible inviscid flows, and equilibrium stress distributions in solids.  Dominated by elliptical equations, equilibrium problems are also called jury problems as the boundary conditions determine the solution of every point in the domain.


Marching Problems:


Other equations of note:


First-Order Linear Wave Equation or the Advection Equation:


\[  \dfrac{\partial u}{\partial t} + c \dfrac{\partial u}{\partial x} = 0 \]

Inviscid Burgers Equation or Nonlinear First-Order Wave Equation:


\[  \dfrac{\partial u}{\partial t} + u \dfrac{\partial u}{\partial x} = 0 \]


Burgers' Equation or Nonlinear Wave Equation with Diffusion:


\[  \dfrac{\partial u}{\partial t} + u \dfrac{\partial u}{\partial x} = v \dfrac{\partial^2 u}{\partial x^2} \]


Tricomi Equation with Diffusion:


\[  y \dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2} = 0 \]

---

References:


Arfken, G. B. & Weber, H.-J. Mathematical Methods for Physicists (6th ed.). Elsevier Academic Press. Burlington, MA. 2005

Chung, T. J. Computational Fluid Dynamics, Cambridge University Press, Cambridge, UK. 2002

Chung, T. J. Computational Fluid Dynamics, 2nd ed. Cambridge University Press, Cambridge, UK. 2010

Hoffman, J. D. Numerical Methods for Engineers and Scientists. CRC Press, Boca Raton, FL, 2nd edition, 2001.

Özişik, M. N. Heat Conduction. John Wiley & Sons, Inc., New York, NY, 2nd edition, 1993.

Smith, G. D. Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford Applied Mathematics & Computing Science Series. Oxford University Press, New York, NY, 3rd edition, 1985.

Tannehill, J. C., Anderson, D. A., and Pletcher, R. H. Computational Fluid Mechanics and Heat Transfer, 2nd ed. Taylor & Francis  Hemisphere, New York, NY. 1997

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 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)