%% This document created by Scientific Word (R) Version 3.0

\documentclass[12pt,thmsa]{article}
\usepackage{sw20jart}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{TCIstyle=article/art4.lat,jart,sw20jart}
%TCIDATA{Created=Sun Feb 13 21:16:24 2000}
%TCIDATA{LastRevised=Mon May 01 09:07:47 2000}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}



\begin{document}
\parindent0cm

%\Large
%

\begin{center}
\textbf{{\Huge Partial Differential Equations}}
\end{center}

\vspace{0.5cm}

\begin{center}
W. E. Schiesser\\[0pt]Iacocca D307\\[0pt]111 Research Drive\\[0pt]Lehigh
University\\[0pt]Bethlehem, PA 18015\\[0pt]

\vspace{0.5cm}

(610) 758-4264 (office)\\[0pt](610) 758-5057 (fax)\\[0pt]

\vspace{0.5cm}

wes1@lehigh.edu \\[0pt]\vspace{0.2cm} http://www.lehigh.edu/\~{\thinspace}wes1/wes1.html
\end{center}

\newpage

\begin{center}
\textbf{{\Huge \textbf{\ Partial Differential Equations}}}
\end{center}

\vspace{0.2cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace}wes1/apci/28apr00.tex\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/28apr00.ps\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/28apr00.pdf\\[0pt]files.txt
\end{center}

\vspace{1cm}

\begin{enumerate}
\item  Geometric and related classifications

\item  Time and spatial discretization

\item  Method of Lines

\item  Finite differences

\item  Finite elements

\item  Finite volumes

\item  Weighted residuals

\item  Adaptive grids

\item  References and software
\end{enumerate}

\newpage

\textbf{\ Why differential equations?}

\vspace{0.5cm}

Differential equations are possibly the most widely used form of mathematics
in science and engineering.

\vspace{0.5cm}

\textbf{\ What are differential equations?}

\vspace{0.5cm}

Differential equations are equations with one or more derivatives.

\vspace{0.5cm}

\textbf{\ What are partial differential equations (PDEs)?}

\vspace{0.5cm}

PDEs are differential equations with two or more independent variables,
typically time and one or more spatial dimensions.

\vspace{0.5cm}

\textbf{\ Why solve PDEs using computers?}

\vspace{0.5cm}

Essentially all realistic problems in differential equations are:

\begin{itemize}
\item  High order (of order $n\times n$)

\begin{itemize}
\item $10^{0}\times10^{0}$ systems are now routine

\item $10^{1}\times10^{1}$ to $10^{2}\times10^{2}$ systems are now rather commonplace

\item $10^{3}\times10^{3}$ systems are at the forefront of applications
\end{itemize}

\item  Nonlinear
\end{itemize}

\vspace{0.5cm}

We can solve differential equations (PDEs) analytically only if $n\leq2$, and
they are linear.

\vspace{0.5cm}

Representative analytical methods include separation of variables, Laplace transforms

\vspace{1cm}

\textbf{\ What is the origin of such high order, nonlinear problems?}

\vspace{0.5cm}

Examples: Solid state device simulation; dynamics of chemical reactors,
separation systems

\vspace{0.5cm}

PDEs - (more than one independent variable, e.g., $x,t)$%

\[
\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial x^{2}};
\text{solution: }u(x,t)
\]

The solution is the dependent variable(s) as a function of the independent
variables (that satisfies the PDE(s) and all of its auxiliary conditions).

\vspace{0.5cm}

PDEs are frequently written with subscript notation, e.g.,%

\[
\frac{\partial u}{\partial t}\Leftrightarrow u_{t},\frac{\partial^{2}%
u}{\partial x^{2}}\Leftrightarrow u_{xx}
\]

so the PDE can be written in subscript notation as%

\[
u_{t}=u_{xx}
\]

This notation makes the statement of PDEs more compact, and is also useful in
their coding, e.g., for an arrays, we can use%

\[
u\Leftrightarrow u(i),u_{t},\Leftrightarrow ut(i),u_{xx}\Leftrightarrow
uxx(i)
\]

\vspace{0.5cm}

In general, we will obtain the solution in numerical form, e.g., $u(x,t)$ will
be tabulated (or plotted) values of $u(x,t)$ vs $x$ and $t$..

\vspace{0.5cm}

\textbf{Geometric and Related Classifications}

\vspace{0.5cm}

The geometric classification of PDEs as \textit{elliptic, hyperbolic and
parabolic} has a rigorous definition for a single linear, second order PDE.
Here we introduce a geometric classification that is less rigorous, but more general.

\vspace{0.5cm}%

\[%
\begin{array}
[c]{cccc}%
\text{Order in }z\text{ (BV)} & \text{Order in }t\text{ (IV)} &
\text{Classification} & \text{Example}\\
&  &  & \\
1 & 1 & \text{First order hyperbolic} & \dfrac{\partial u}{\partial
t}=-v\dfrac{\partial u}{\partial z}\\
&  &  & \text{ (advection equation)}\\
&  &  & \\
2 & 2 & \text{Second order hyperbolic} & \dfrac{\partial^{2}u}{\partial t^{2}%
}=c^{2}\dfrac{\partial^{2}u}{\partial z^{2}}\\
&  &  & \text{ (wave equation)}\\
&  &  & \\
2 & 1 & \text{Parabolic} & \dfrac{\partial u}{\partial t}=\alpha
\dfrac{\partial^{2}u}{\partial z^{2}}\\
&  &  & \text{ (Fourier's or Fick's second law)}\\
&  &  & \\
2\text{ }(\text{in }y,z\text{)} & 0 & \text{Elliptic} & \dfrac{\partial^{2}%
u}{\partial y^{2}}+\dfrac{\partial^{2}u}{\partial z^{2}}=0\\
&  &  & \text{ (Laplace's equation)}%
\end{array}
\]

Combinations of these classes of PDEs are also possible. For example%

\begin{equation}
\dfrac{\partial u}{\partial t}=-v\dfrac{\partial u}{\partial z}+D\dfrac
{\partial^{2}u}{\partial z^{2}}\text{ }\tag{1}%
\end{equation}

is a \textit{hyperbolic-parabolic (or convective diffusion)} equation which
could have the auxiliary conditions%

\begin{equation}
u(z,0)=f(z)\tag{2}%
\end{equation}%

\begin{equation}
-D\dfrac{\partial u(0,t)}{\partial z}=v(u_{0}-u(0,t))\tag{3}%
\end{equation}%

\begin{equation}
\dfrac{\partial u(L,t)}{\partial z}=0\tag{4}%
\end{equation}

BCs (3) and (4) (note there are two since eq. (1) is second order in $z$) are
the \textbf{Wilhelm Danckwerts} BCs (named after two well known ChEs);
generally, eq. (4) is not appropriate for unsteady state problems, and can be
replaced by%

\begin{equation}
\dfrac{\partial u(L,t)}{\partial t}=-v\dfrac{\partial u(L,t)}{\partial
z}\tag{5}%
\end{equation}

Note that eq. (5) is an acceptable BC for eq. (1) since it is only first order
in $z$, while eq. (1) is second order in $z$ (generally, a BC is at least one
order lower than the highest order derivative in the PDE).

\vspace{0.5cm}

To repeat, we have the basic question, what do we mean by the solution to a
PDE? The answer in general is \textit{the dependent variable as a function of
the independent variables}. Thus, for equations (1) to (4), the solution is a
function, $u(z,t),$ in analytical or numerical form.

\vspace{0.5cm}

Eq. (1) can be extended in \textbf{cylindrical coordinates} to include a
radial coordinate, $r$%

\begin{equation}
\dfrac{\partial u}{\partial t}=-v\dfrac{\partial u}{\partial z}+D\left(
\dfrac{\partial^{2}u}{\partial z^{2}}\text{ }+\dfrac{\partial^{2}u}{\partial
r^{2}}\text{ }+\frac{1}{r}\dfrac{\partial u}{\partial r}\right) \tag{6}%
\end{equation}

Since eq. (6) is second order in $r$, it requires two BCs, e.g.,%

\begin{equation}
\dfrac{\partial u(0,z,t)}{\partial r}=0\tag{7}%
\end{equation}

Eq. (7) is a \textit{Neumann BC} (reflecting symmetry around the centerline
$r=0$) since it defines the first derivative of the solution, $u(r,z,t)$, at
$r=0$.

\vspace{0.5cm}

A second BC might be
\begin{equation}
u(r_{0},z,t)=u_{w}\tag{8}%
\end{equation}

Eq. (8) is a \textit{Dirichlet BC} since it defines the dependent variable,
$u(r,z,t),$ at the wall, $r=r_{0}.$

\vspace{0.5cm}

Eq. (3) is a combination of a Dirichlet BC and a Neumann BC (it contains both
$u(0,t)$ and $\dfrac{\partial u(0,t)}{\partial z}$), and is generally termed a
\textit{BC of the third type or a Robin BC}.

\vspace{0.5cm}

Thus, the possible BCs are \textit{Dirichlet, Neumann, third type, linear or
nonlinear}.

\vspace{0.5cm}

The coefficient $\dfrac{1}{r}$ in eq. (6) is termed a \textit{variable
coefficient} (it is a function of the independent variable $r$). Note,
however, that eq. (6) is linear (PDEs with variable coefficients may be linear
or nonlinear).

\vspace{0.5cm}

If the third angular coordinate is added to eq. (6), it becomes%

\begin{equation}
\dfrac{\partial u}{\partial t}=-v\dfrac{\partial u}{\partial z}+D\left(
\dfrac{\partial^{2}u}{\partial z^{2}}\text{ }+\dfrac{\partial^{2}u}{\partial
r^{2}}\text{ }+\frac{1}{r}\dfrac{\partial u}{\partial r}+\frac{1}{r^{2}}%
\dfrac{\partial^{2}u}{\partial\theta^{2}}\right) \tag{9}%
\end{equation}

Eq. (9) is a \textit{linear, initial value, hyperbolic-parabolic PDE (with
variable coefficients) in cylindrical coordinates}; it requires two BCs in
$\theta$ (generally unnecessary in most applications because of no variation
of $u(r,\theta,z,t)$ with $\theta$ (note the ``units'' of the various
derivatives in eq. (9)). Eq. (9), with $v=0$, is the full \textit{Fourier's
second law in cylindrical coordinates} (a parabolic PDE). Note that eq. (9) is
first order in $t,$ and second order in $r,\theta$ and $z$, and therefore
requires the requisite number of auxiliary conditions in each independent variable

\vspace{0.5cm}%

\[%
\begin{array}
[c]{cccc}%
\text{Independent} & \text{Highest order} & \text{Type of auxiliary} &
\text{Example}\\
\text{variable} & \text{derivative} & \text{condition} & \\
t & 1 & \text{Initial condition} & u(r,\theta,z,0)=u_{0}\\
&  &  & \\
r & 2 & \text{Neumann, linear} & \partial u(r_{0},\theta,z,t)/\partial r=0\\
&  & \text{Third type, nonlinear} & -k\partial u(r_{1},\theta,z,t)/\partial
r=\\
&  &  & \epsilon(u_{a}^{4}-u^{4}(r_{1},\theta,z,t))\\
&  &  & \\
\theta & 2 & \text{Neumann, linear} & \partial u(r,0,z,t)/\partial\theta=0\\
&  & \text{Neumann, linear} & \partial u(r,\pi,z,t)/\partial\theta=0\\
&  &  & \\
z & 2 & \text{Neumann, linear} & \partial u(r,\theta,0,t)/\partial z=0\\
&  & \text{Neumann, linear} & \partial u(r,\theta,L,t)/\partial z=0
\end{array}
\]

\vspace{0.5cm}

These ICs and BCs imply ranges for the independent variables that are not
specified above. For example, the second BC%

\[
\partial u(r_{0},\theta,z,t)/\partial r=0
\]

means%

\[
\partial u(r_{0},\theta,z,t)/\partial r=0,0\leq\theta\leq\pi,0\leq z\leq
L,t\geq0
\]

The equation in spherical coordinates is:%

\begin{equation}
\dfrac{\partial u}{\partial t}=D\left\{  \dfrac{\partial^{2}u}{\partial r^{2}%
}\text{ }+\frac{2}{r}\dfrac{\partial u}{\partial r}+\frac{1}{r^{2}\sin\theta
}\frac{\partial}{\partial\theta}\left(  \sin\theta\frac{\partial u}%
{\partial\theta}\right)  +\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}%
u}{\partial\varphi^{2}}\right\} \tag{10}%
\end{equation}

Eq. (10) is a \textit{linear, initial value, parabolic PDE (with variable
coefficients) in spherical coordinates} (Fourier's second law in spherical coordinates).

\vspace{0.5cm}

If we have a reaction in the sphere, with a significant heat effect%

\begin{equation}
\dfrac{\partial u}{\partial t}=D\left\{  \dfrac{\partial^{2}u}{\partial r^{2}%
}\text{ }+\frac{2}{r}\dfrac{\partial u}{\partial r}+\frac{1}{r^{2}\sin\theta
}\frac{\partial}{\partial\theta}\left(  \sin\theta\frac{\partial u}%
{\partial\theta}\right)  +\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}%
u}{\partial\varphi^{2}}\right\}  +\Delta Hk_{0}e^{-E/(Ru)}\tag{11}%
\end{equation}

Eq. (10) is a \textit{nonlinear, initial value, parabolic PDE (with variable
coefficients) in spherical coordinates.}

\vspace{0.5cm}

All of these PDEs can be stated in a coordinate-independent format, e.g.,%

\begin{equation}
\dfrac{\partial u}{\partial t}=D\nabla\cdot\nabla u-\nabla\cdot(\mathbf{v}%
u)-R\tag{12}%
\end{equation}

where:

\newpage

$\nabla\cdot$ (\textit{divergence of a vector}):

\vspace{0.5cm}%

\[%
\begin{array}
[c]{cc}%
\text{Coordinate system} & \text{Components}\\
\text{Cartesian} & \left[
\begin{array}
[c]{c}%
\left[  \nabla\right]  _{x}=\dfrac{\partial}{\partial x}\\
\\
\left[  \nabla\right]  _{y}=\dfrac{\partial}{\partial y}\\
\\
\left[  \nabla\right]  _{z}=\dfrac{\partial}{\partial z}%
\end{array}
\right] \\
& \\
\text{cylindrical} & \left[
\begin{array}
[c]{c}%
\left[  \nabla\right]  _{r}=\dfrac{1}{r}\dfrac{\partial}{\partial r}(r\text{
})\\
\\
\left[  \nabla\right]  _{\theta}=\dfrac{1}{r}\dfrac{\partial}{\partial\theta
}\\
\\
\left[  \nabla\right]  _{z}=\dfrac{\partial}{\partial z}%
\end{array}
\right] \\
& \\
\text{spherical} & \left[
\begin{array}
[c]{c}%
\left[  \nabla\right]  _{r}=\dfrac{1}{r^{2}}\dfrac{\partial}{\partial r}%
(r^{2}\text{ })\\
\\
\left[  \nabla\right]  _{\theta}=\dfrac{1}{r\sin\theta}\dfrac{\partial
}{\partial\theta}(\sin\theta\text{ })\\
\\
\left[  \nabla\right]  _{\phi}=\dfrac{1}{r\sin\theta}\dfrac{\partial}%
{\partial\phi}%
\end{array}
\right]
\end{array}
\]

\newpage

$\nabla$ (\textit{gradient of a scalar)} :%

\[%
\begin{array}
[c]{cc}%
\text{Coordinate system} & \text{Components}\\
\text{Cartesian } & \left[
\begin{array}
[c]{c}%
\left[  \nabla\right]  _{x}=\dfrac{\partial}{\partial x}\\
\\
\left[  \nabla\right]  _{y}=\dfrac{\partial}{\partial y}\\
\\
\left[  \nabla\right]  _{z}=\dfrac{\partial}{\partial z}%
\end{array}
\right] \\
& \\
\text{cylindrical} & \left[
\begin{array}
[c]{c}%
\left[  \nabla\right]  _{r}=\dfrac{\partial}{\partial r}\\
\\
\left[  \nabla\right]  _{\theta}=\dfrac{1}{r}\dfrac{\partial}{\partial\theta
}\\
\\
\left[  \nabla\right]  _{z}=\dfrac{\partial}{\partial z}%
\end{array}
\right] \\
& \\
\text{spherical} & \left[
\begin{array}
[c]{c}%
\left[  \nabla\right]  _{r}=\dfrac{\partial}{\partial r}\\
\\
\left[  \nabla\right]  _{\theta}=\dfrac{1}{r}\dfrac{\partial}{\partial\theta
}\\
\\
\left[  \nabla\right]  _{\phi}=\dfrac{1}{r\sin\theta}\dfrac{\partial}%
{\partial\phi}%
\end{array}
\right]
\end{array}
\]

\newpage

$\nabla\cdot\nabla$ (\textit{divergence of the gradient of a scalar}):

\vspace{0.5cm}%

\[%
\begin{array}
[c]{cc}%
\text{Coordinate system} & \text{Component}\\
\text{Cartesian} & \dfrac{\partial^{2}}{\partial x^{2}}+\dfrac{\partial^{2}%
}{\partial y^{2}}+\dfrac{\partial^{2}}{\partial z^{2}}\\
& \\
\text{cylindrical} & (\dfrac{\partial^{2}}{\partial r^{2}}+\dfrac{1}{r}%
\dfrac{\partial}{\partial r})+\dfrac{1}{r^{2}}\dfrac{\partial^{2}}%
{\partial\theta^{2}}+\dfrac{\partial^{2}}{\partial z^{2}}\\
& \\
\text{spherical} & \dfrac{1}{r^{2}}\dfrac{\partial}{\partial r}(r^{2}\text{
}\dfrac{\partial}{\partial r})+\dfrac{1}{r^{2}\sin\theta}\dfrac{\partial
}{\partial\theta}(\sin\theta\text{ }\dfrac{\partial}{\partial\theta}%
)+\dfrac{1}{r^{2}\sin^{2}\theta}\dfrac{\partial^{2}}{\partial\phi^{2}}%
\end{array}
\]

The derivation of $\nabla\cdot\nabla$ (the Laplacian) follows directly from
the preceding components of $\nabla\cdot$ (divergence of a vector) and
$\nabla$ (gradient of a scalar).

\vspace{0.5cm}

Cartesian coordinates:%

\[
\nabla\cdot\nabla=(\mathbf{i}\dfrac{\partial}{\partial x}+\mathbf{j}%
\dfrac{\partial}{\partial y}+\mathbf{k}\dfrac{\partial}{\partial z}%
)\cdot(\mathbf{i}\dfrac{\partial}{\partial x}+\mathbf{j}\dfrac{\partial
}{\partial y}+\mathbf{k}\dfrac{\partial}{\partial z})=\dfrac{\partial^{2}%
}{\partial x^{2}}+\dfrac{\partial^{2}}{\partial y^{2}}+\dfrac{\partial^{2}%
}{\partial z^{2}}
\]

Cylindrical coordinates:%

\[
\nabla\cdot\nabla=(\mathbf{i}_{r}\frac{1}{r}\dfrac{\partial}{\partial
r}(r)+\mathbf{j}_{\theta}\dfrac{1}{r}\dfrac{\partial}{\partial\theta
}+\mathbf{k}_{z}\dfrac{\partial}{\partial z})\cdot(\mathbf{i}_{r}%
\dfrac{\partial}{\partial r}+\mathbf{j}_{\theta}\dfrac{1}{r}\dfrac{\partial
}{\partial\theta}+\mathbf{k}_{z}\dfrac{\partial}{\partial z})
\]%

\[
=\frac{1}{r}\dfrac{\partial}{\partial r}(r\dfrac{\partial}{\partial r}%
)+\dfrac{1}{r}\dfrac{\partial}{\partial\theta}(\dfrac{1}{r}\dfrac{\partial
}{\partial\theta})+\dfrac{\partial}{\partial z}(\dfrac{\partial}{\partial z})
\]%

\[
=\frac{1}{r}(\dfrac{\partial}{\partial r}+r\dfrac{\partial^{2}}{\partial
r^{2}})+\dfrac{1}{r^{2}}\dfrac{\partial}{\partial\theta}\dfrac{\partial
}{\partial\theta}+\dfrac{\partial}{\partial z}\dfrac{\partial}{\partial z}
\]%

\[
=(\dfrac{\partial^{2}}{\partial r^{2}}+\frac{1}{r}\dfrac{\partial}{\partial
r})+\dfrac{1}{r^{2}}\dfrac{\partial^{2}}{\partial\theta^{2}}+\dfrac
{\partial^{2}}{\partial z^{2}}
\]

Spherical coordinates:%

\[
\nabla\cdot\nabla=(\mathbf{i}_{r}\dfrac{1}{r^{2}}\dfrac{\partial}{\partial
r}(r^{2}\text{ })+\mathbf{j}_{\theta}\dfrac{1}{r\sin\theta}\dfrac{\partial
}{\partial\theta}(\sin\theta\text{ })+\mathbf{k}_{\phi}\dfrac{1}{r\sin\theta
}\dfrac{\partial}{\partial\phi})\cdot(\mathbf{i}_{r}\dfrac{\partial}{\partial
r}+\mathbf{j}_{\theta}\dfrac{1}{r}\dfrac{\partial}{\partial\theta}%
+\mathbf{k}_{\phi}\dfrac{1}{r\sin\theta}\dfrac{\partial}{\partial\phi})
\]%

\[
=\dfrac{1}{r^{2}}\dfrac{\partial}{\partial r}(r^{2}\text{ }\dfrac{\partial
}{\partial r})+\dfrac{1}{r\sin\theta}\dfrac{\partial}{\partial\theta}%
(\sin\theta\text{ }\dfrac{1}{r}\dfrac{\partial}{\partial\theta})+\dfrac
{1}{r\sin\theta}\dfrac{\partial}{\partial\phi}(\dfrac{1}{r\sin\theta}%
\dfrac{\partial}{\partial\phi})
\]%

\[
=\dfrac{1}{r^{2}}\dfrac{\partial}{\partial r}(r^{2}\text{ }\dfrac{\partial
}{\partial r})+\dfrac{1}{r^{2}\sin\theta}\dfrac{\partial}{\partial\theta}%
(\sin\theta\text{ }\dfrac{\partial}{\partial\theta})+\dfrac{1}{r^{2}\sin
^{2}\theta}\dfrac{\partial^{2}}{\partial\phi^{2}}
\]

To summarize, the classification of a PDE and its auxiliary conditions is
specified in terms of

\vspace{0.5cm}

\begin{itemize}
\item  Geometric classification: elliptic, hyperbolic, parabolic

\item  Coordinate system

\item  Number of dimensions

\item  Initial value or boundary value

\item  Linearity

\item  Coefficients: constant or variable

\item  Boundary conditions: Dirichlet, Neumann, third type; linear or nonlinear
\end{itemize}

\vspace{0.5cm}

This classification is useful for:

\begin{itemize}
\item  Defining the type of PDE(s)

\item  Suggesting numerical methods for the solution of the PDE(s)
\end{itemize}

\vspace{0.5cm}

Of the three classes of PDEs, hyperbolic PDEs are the most difficult to
integrate numerically. This can be illustrated with some famous examples of
first-order, hyperbolic PDEs:

\vspace{0.5cm}

(1) \textbf{\ }\textit{One dimensional, isothermal Euler equations written in
conservation form}:%

\[
\mathbf{u}_{t}+\left[  \mathbf{f}(\mathbf{u})\right]  _{x}=\mathbf{0}
\]

where%

\[
\mathbf{u}=\left[
\begin{array}
[c]{c}%
\rho\\
\rho u
\end{array}
\right]  ,\mathbf{f}(\mathbf{u})=\left[
\begin{array}
[c]{c}%
\rho u\\
\rho u^{2}+p
\end{array}
\right]
\]

Substitution of $\mathbf{u}$ and $\mathbf{f}(\mathbf{u})$ into the general
equation gives the one-dimensional Euler equations%

\[%
\begin{array}
[c]{cc}%
\rho_{t}+(\rho u)_{x}=0 & \text{continuity}\\
& \\
(\rho u)_{t}+(\rho u^{2}+p)_{x}=0 & \text{momentum}%
\end{array}
\]

Also, an equation of state (EOS) is required (to give a $2\times2$ system for
$\rho(x,t),$ $u(x,t)$)%

\[
p=p(\rho,T)
\]

These PDEs are:

\vspace{0.5cm}

\begin{itemize}
\item  Hyperbolic (first order in $x$ and $t$)

\item  Source of sharp fronts and discontinuities

\item  Nonlinear

\item  Difficult to solve numerically (and impossible analytically)

\item  Limited in their application, e.g., they have no heat transfer or
reaction terms
\end{itemize}

\vspace{0.5cm}

CFD codes are available for the solution of these Euler equations.

\vspace{0.5cm}

(2) \textit{One dimensional, nonisothermal Euler equations with chemical
reaction written in conservation form}:%

\[
\mathbf{u}_{t}+\left[  \mathbf{f}(\mathbf{u})\right]  _{x}=\mathbf{g(u)}
\]

where%

\[
\mathbf{u}=\left[
\begin{array}
[c]{c}%
\rho\\
\rho u\\
E\\
\rho Y_{1}\\
\vdots\\
\rho Y_{N}%
\end{array}
\right]  ,\mathbf{f}(\mathbf{u})=\left[
\begin{array}
[c]{c}%
\rho u\\
\rho u^{2}+p\\
(E+p)u\\
\rho uY_{1}\\
\vdots\\
\rho uY_{N}%
\end{array}
\right]  ,\mathbf{g}(\mathbf{u})=\left[
\begin{array}
[c]{c}%
0\\
0\\
0\\
r_{1}(T,\rho,Y_{1},\cdots,Y_{N})\\
\vdots\\
r_{N}(T,\rho,Y_{1},\cdots,Y_{N})
\end{array}
\right]
\]

where $Y_{i}$ is the mass fraction of component $i$ and $N$ is the number of
reacting species; the energy is given by%

\[
E=-p+\frac{\rho u^{2}}{2}+\rho h
\]

where $h$ is the enthalpy per unit mass.

\vspace{0.5cm}

Substitution of $\mathbf{u}$, $\mathbf{f}(\mathbf{u})$ and $\mathbf{g}%
(\mathbf{u})$ into the general equation gives the one-dimensional Euler
equations with chemical reactions%

\[%
\begin{array}
[c]{cc}%
\rho_{t}+(\rho u)_{x}=0 & \text{total continuity}\\
& \\
(\rho u)_{t}+(\rho u^{2}+p)_{x}=0 & \text{momentum}\\
& \\
E_{t}+\left[  (E+p)u\right]  _{x}=0 & \text{energy}\\
& \\
\left(  \rho Y_{1}\right)  _{t}+\left[  \rho uY_{1}\right]  _{x}=r_{1}%
(T,\rho,Y_{1},\cdots,Y_{N}) & \text{component 1 continuity}\\
\vdots & \vdots\\
\left(  \rho Y_{N}\right)  _{t}+\left[  \rho uY_{N}\right]  _{x}=r_{N}%
(T,\rho,Y_{1},\cdots,Y_{N}) & \text{component N continuity}%
\end{array}
\]

Note that an explicit heat transfer term and a heat of reaction term do not
appear in the energy balance (these can be added as in the model tubular
reactor PDEs derived in the first lecture). Again, an equation of state (EOS)
is required%

\[
p=p(\rho,E,Y_{1},\cdots Y_{N}),T=T(\rho,E,Y_{1},\cdots Y_{N})
\]

This is a $(3+N)\times(3+N)$ for $\rho(x,t),u(x,t),E(x,t),Y_{1}(x,t),\cdots
Y_{N}(x,t)$.

\vspace{0.5cm}

CFD codes are available for the solution of these Euler equations.

\vspace{0.5cm}

(3) \textit{Multidimensional, isothermal Euler equations (without the energy
equation or chemical reaction) written in conservation form:}%

\[
\mathbf{u}_{t}+\nabla\cdot\left[  \mathbf{f}(\mathbf{u})\right]  =\mathbf{0}
\]

where%

\[
\mathbf{u}=\left[
\begin{array}
[c]{c}%
\rho\\
\rho\mathbf{v}%
\end{array}
\right]  ,\mathbf{f}(\mathbf{u})=\left[
\begin{array}
[c]{c}%
\rho v\\
\rho\mathbf{vv}+p
\end{array}
\right]
\]

$\mathbf{v}$ is a vector, e.g., Substitution of $\mathbf{u}$ and
$\mathbf{f}(\mathbf{u})$ into the general equation gives the one-dimensional
Euler equations%

\[%
\begin{array}
[c]{cc}%
\rho_{t}+\nabla\cdot(\rho\mathbf{v})=0 & \text{continuity}\\
& \\
(\rho\mathbf{v})_{t}+\nabla\cdot(\rho\mathbf{vv)}+\nabla p=0 & \text{momentum}%
\end{array}
\]

where $\nabla$ is an operator applied to a vector in the continuity equation,
and to a scalar and a tensor in the momentum equation. Considering $\nabla$ a
bit further, in Cartesian coordinates (with $\mathbf{v}=\mathbf{i}%
v_{x}+\mathbf{j}v_{y}+\mathbf{k}v_{z})$%

\[
\nabla=\mathbf{i}\frac{\partial}{\partial x}+\mathbf{j}\frac{\partial
}{\partial y}+\mathbf{k}\frac{\partial}{\partial z}\text{ }
\]

For example%

\[
\nabla\cdot(\rho\mathbf{v})=\left(  \mathbf{i}\frac{\partial}{\partial
x}+\mathbf{j}\frac{\partial}{\partial y}+\mathbf{k}\frac{\partial}{\partial
z}\right)  \cdot\left(  \mathbf{i}\rho v_{x}+\mathbf{j}\rho v_{y}%
+\mathbf{k}\rho v_{z}\right)
\]%

\[
=\frac{\partial\left(  \rho v_{x}\right)  }{\partial x}+\frac{\partial(\rho
v_{y})}{\partial y}+\frac{\partial\left(  \rho v_{z}\right)  }{\partial z}
\]

and the (scalar) continuity equation is%

\[
\frac{\partial\rho}{\partial t}+\frac{\partial\left(  \rho v_{x}\right)
}{\partial x}+\frac{\partial(\rho v_{y})}{\partial y}+\frac{\partial\left(
\rho v_{z}\right)  }{\partial z}=0
\]

Note that for a constant density (incompressible) fluid,%

\[
\rho_{t}=0
\]%

\[
\frac{\partial\left(  \rho v_{x}\right)  }{\partial x}+\frac{\partial(\rho
v_{y})}{\partial y}+\frac{\partial\left(  \rho v_{z}\right)  }{\partial
z}=\rho(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial
y}+\frac{\partial v_{z}}{\partial z})+v_{x}\frac{\partial\rho}{\partial
x}+v_{y}\frac{\partial\rho}{\partial x}+v_{z}\frac{\partial\rho}{\partial
x}=0
\]

or%

\[
\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}%
+\frac{\partial v_{z}}{\partial z}=0
\]%

\[
\nabla\cdot\mathbf{v}=0
\]

Thus, an incompressible fluid is \textit{divergence free}.

\vspace{0.5cm}

Consider the term%

\[
\nabla\cdot(\rho\mathbf{vv)}
\]

The product $\mathbf{vv}$ is actually a second order (nine-component) tensor.
Thus
\[
\rho\mathbf{vv=}\left[
\begin{array}
[c]{ccc}%
\rho v_{x}v_{x} & \rho v_{x}v_{y} & \rho v_{x}v_{z}\\
\rho v_{y}v_{x} & \rho v_{y}v_{y} & \rho v_{y}v_{z}\\
\rho v_{z}v_{x} & \rho v_{z}v_{y} & \rho v_{z}v_{z}%
\end{array}
\right]
\]%

\[
\nabla\cdot(\rho\mathbf{vv)=}\left[
\begin{array}
[c]{c}%
\mathbf{i}\left\{  \dfrac{\partial}{\partial x}v_{x}(\rho v_{x})+\dfrac
{\partial}{\partial y}v_{y}(\rho v_{x})+\dfrac{\partial}{\partial z}v_{z}(\rho
v_{x})\right\} \\
\\
\mathbf{j}\left\{  \dfrac{\partial}{\partial x}v_{x}(\rho v_{y})+\dfrac
{\partial}{\partial y}v_{y}(\rho v_{y})+\dfrac{\partial}{\partial z}v_{z}(\rho
v_{y})\right\} \\
\\
\mathbf{k}\left\{  \dfrac{\partial}{\partial x}v_{x}(\rho v_{z})+\dfrac
{\partial}{\partial y}v_{y}(\rho v_{z})+\dfrac{\partial}{\partial z}v_{z}(\rho
v_{z})\right\}
\end{array}
\right]
\]

Then, using this result (a vector) in the momentum equation%

\[
(\rho\mathbf{v})_{t}+\nabla\cdot(\rho\mathbf{vv)}+\nabla p=0
\]

and equating like components%

\[
\dfrac{\partial}{\partial t}\left(  \rho v_{x}\right)  +\dfrac{\partial
}{\partial x}(\rho v_{x}v_{x})+\dfrac{\partial}{\partial y}(\rho v_{x}%
v_{y})+\dfrac{\partial}{\partial z}(\rho v_{x}v_{z})+\frac{\partial
P}{\partial x}=0\text{ (}\mathbf{i}\text{ component)}
\]%

\[
\dfrac{\partial}{\partial t}(\rho v_{y})+\dfrac{\partial}{\partial x}(\rho
v_{y}v_{x})+\dfrac{\partial}{\partial y}(\rho v_{y}v_{y})+\dfrac{\partial
}{\partial z}(\rho v_{y}v_{z})+\frac{\partial P}{\partial y}=0\text{
(}\mathbf{j}\text{ component)}
\]%

\[
\dfrac{\partial}{\partial t}\left(  \rho v_{z}\right)  +\dfrac{\partial
}{\partial x}(\rho v_{z}v_{x})+\dfrac{\partial}{\partial y}(\rho v_{z}%
v_{y})+\dfrac{\partial}{\partial z}(\rho v_{z}v_{z})+\frac{\partial
P}{\partial z}=0\text{ (}\mathbf{k}\text{ component)}
\]

Again, these PDEs are:

\begin{itemize}
\item  Hyperbolic (first order in $x,y,z$ and $t$)

\item  Source of sharp fronts and discontinuities (but now in three dimensions)

\item  Nonlinear

\item  Difficult to solve numerically (and impossible analytically)

\item  Limited in their application, e.g., they have no heat transfer,
reaction terms
\end{itemize}

\vspace{0.5cm}

CFD codes are available for the solution of these Euler equations.

\vspace{0.5cm}

These equations can be simplified through the continuity equation, e.g., for
the $\mathbf{i}$ component:%

\[
v_{x}\dfrac{\partial}{\partial t}\left(  \rho\right)  +v_{x}\dfrac{\partial
}{\partial x}(\rho v_{x})+v_{x}\dfrac{\partial}{\partial y}(\rho v_{y}%
)+v_{x}\dfrac{\partial}{\partial z}(\rho v_{z})
\]%

\[
\rho\dfrac{\partial}{\partial t}\left(  v_{x}\right)  +\rho v_{x}%
\dfrac{\partial}{\partial x}(\rho v_{x})+\rho v_{y}\dfrac{\partial}{\partial
y}(\rho v_{x})+\rho v_{z}\dfrac{\partial}{\partial z}(\rho v_{x}%
)+\frac{\partial P}{\partial x}=0
\]

or%

\[
v_{x}\left\{  \dfrac{\partial\rho}{\partial t}+\dfrac{\partial(\rho v_{x}%
)}{\partial x}+\dfrac{\partial(\rho v_{y})}{\partial y}+\dfrac{\partial(\rho
v_{z})}{\partial z}\right\}
\]
\[
\rho\dfrac{\partial v_{x}}{\partial t}+\rho v_{x}\dfrac{\partial v_{x}%
}{\partial x}+\rho v_{y}\dfrac{\partial v_{x}}{\partial y}+\rho v_{z}%
\dfrac{\partial v_{x}}{\partial z}+\frac{\partial P}{\partial x}=0
\]

and the first row is zero through the continuity equation. In general,%

\[
\rho\mathbf{v}_{t}+\rho\mathbf{v\cdot}\nabla\mathbf{v}+\nabla p=0
\]

which is an alternate form of the Euler equations.

\vspace{0.5cm}

Finally, if a viscosity term is added to the momentum equation%

\[
(\rho\mathbf{v})_{t}+\nabla\cdot(\rho\mathbf{vv)}+\nabla p-\mu\nabla
^{2}\mathbf{v}=0
\]

which are the \textit{Navier Stokes equations.}

\vspace{0.5cm}

The preceding PDEs are written in Cartesian coordinates. As another general
approach, we can start with PDEs written in terms of vector operators that are
independent of the coordinate system, e.g.,%

\[%
\begin{array}
[c]{cc}%
\rho_{t}+\nabla\cdot(\rho\mathbf{v})=0 & \text{continuity}\\
& \\
(\rho\mathbf{v})_{t}+\nabla\cdot(\rho\mathbf{vv)}+\nabla p=0 & \text{momentum
(Euler)}\\
& \\
\rho\mathbf{v}_{t}+\rho\mathbf{v\cdot}\nabla\mathbf{v}+\nabla p=0 &
\text{momentum (Euler)}\\
& \\
(\rho\mathbf{v})_{t}+\nabla\cdot(\rho\mathbf{vv)}+\nabla p-\mu\nabla
^{2}\mathbf{v}=0 & \text{momentum (Navier Stokes)}\\
& \\
\rho\mathbf{v}_{t}+\rho\mathbf{v\cdot}\nabla\mathbf{v}+\nabla p-\mu\nabla
^{2}\mathbf{v}=0 & \text{momentum (Navier Stokes)}%
\end{array}
\]

With these coordinate-independent equations as the starting point, we can go
to the three preceding tables to write equations in Cartesian, cylindrical or
spherical coordinates. This approach can also be applied to other coordinate
systems, perhaps fitted to an body that is not easily described in one of the
well-established coordinate systems (body-fitted coordinates are an active
area of CFD).

\vspace{0.5cm}

As an example of this simplification process, consider the x-component of the
Navier Stokes equations

\vspace{0.5cm}%

\[
\rho\dfrac{\partial v_{x}}{\partial t}+\rho v_{x}\dfrac{\partial v_{x}%
}{\partial x}+\rho v_{y}\dfrac{\partial v_{x}}{\partial y}+\rho v_{z}%
\dfrac{\partial v_{x}}{\partial z}+\frac{\partial P}{\partial x}-\mu
(\dfrac{\partial^{2}v_{x}}{\partial x^{2}}+\dfrac{\partial^{2}v_{x}}{\partial
y^{2}}+\dfrac{\partial^{2}v_{x}}{\partial z^{2}})=0
\]

If we consider a one-dimensional problem with no pressure gradient, this
equation reduces to%

\[
\rho\dfrac{\partial v_{x}}{\partial t}+\rho v_{x}\dfrac{\partial v_{x}%
}{\partial x}-\mu\dfrac{\partial^{2}v_{x}}{\partial x^{2}}=0
\]

which is \textit{Burgers' equation}. This is an unusual (and widely studied)
PDE because

\vspace{0.5cm}

\begin{itemize}
\item  It is nonlinear, yet has a known, exact (analytical solution)

\item  The solution exhibits moving fronts that can be made arbitarily sharp
by decreasing the kinematic viscosity $\upsilon=\mu/\rho.$
\end{itemize}

\vspace{0.5cm}

Burgers' equation for $\upsilon=\mu/\rho=0$ reduces to the \textit{inviscid
Burgers' equation} or \textit{nonlinear advection equation}.%

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

This equation can propagate discontinuties, and therefore provides a stringent
test problem for PDEs.

\vspace{0.5cm}

As general as all of the preceding equations may appear, they are actually
rather limited, e.g., they do not contain heat transfer or heat of reaction
terms. As an added thought, general purpose codes that solve these equations
may not be as general as we would like. We should at least be sure that we
understand the underlying PDEs of these codes.

\vspace{0.5cm}

Therefore, a second approach to developing a PDE model is to start with a
basic incremental volume and write the appropriate conservation equations for
it (including all of the effects that are considered important). In other
words, derive the mathematical model first, then look for, or write, a code
that solves the model equations. The process should not be reversed, i.e., we
should not start with a general code, then modify the model so that it fits
into the code (and thereby possibly drop important terms that should be included).

\vspace{0.5cm}

This might seem like an inefficient way to proceed, i.e., not using an
(expensive), readily available code, but rather writing a new code for the
particular mathematical model. But the latter ensures that we are solving the
appropriate model, and not one that has been modified to fit into the software
(or worse, using an available code incorrectly by not solving the equations we
think we are solving). In other words, the equations are the fundamental
starting point, and the calculations (software) to solve the equations are secondary.

\vspace{0.5cm}

Additionally, we can examine a code written for the particular problem system
in whatever detail is necessary to fully understand the calculations; a
library code, however, may not provide the details of the calculations in the
detail required to fully understand the calculations. As a minimum, we should
fully understand what the library code is doing.

\vspace{0.5cm}

As a specific example of this second approach, consider the PDE for a heated
pipe that originates from an energy balance written on a section of length
$\Delta z$. This obviously is a simple physical system, but it leads to a PDE
that is surprisingly difficult to solve numerically, depending on the
conditions of the problem.%

\[
A\Delta z\rho C_{p}\frac{\partial T}{\partial t}=Av\rho C_{p}T_{z-\Delta
z}-Av\rho C_{p}T_{z}+\pi D\Delta zh(T_{w}-T)
\]

Division by $A\Delta z\rho C_{p}$ with $A=\pi D^{2}/4$ and minor rearrangement gives%

\[
\frac{\partial T}{\partial t}=-v\frac{\left(  T_{z}-T_{z-\Delta z}\right)
}{\Delta z}+\frac{\pi Dh}{(\pi D^{2}/4)\rho C_{p}}(T_{w}-T)
\]

In the limit as $\Delta z\rightarrow0$, we arrive at the PDE%

\begin{equation}
\frac{\partial T}{\partial t}=-v\frac{\partial T}{\partial z}+\frac{4h}{D\rho
C_{p}}(T_{w}-T)\tag{1}%
\end{equation}

The initial condition (IC) is%

\begin{equation}
T(z,0)=T_{0}\tag{2}%
\end{equation}

and the boundary condition (BC) is%

\begin{equation}
T(0,t)=T_{e}\tag{3}%
\end{equation}

Eq. (1) is a first order, linear hyperbolic PDE. Eqs. (1) to (3) provide a
stringent test problem for numerical methods and software. We will use the
analytical solution to evaluate the numerical methods.

\vspace{0.5cm}

If%

\begin{equation}
u=T-T_{w},\alpha=\frac{4h}{D\rho C_{p}}\tag{4}%
\end{equation}

eqs. (1) to (3) become%

\begin{equation}
\frac{\partial u}{\partial t}=-v\frac{\partial u}{\partial z}-\alpha u\tag{5}%
\end{equation}%

\begin{equation}
u(z,0)=T_{0}-T_{w}\tag{6}%
\end{equation}%

\begin{equation}
u(0,t)=T_{e}-T_{w}\tag{7}%
\end{equation}

If the Laplace transform of $u(z,t)$ with respect to $t$ is defined as%

\[
L_{t}\left\{  u(z,t)\right\}  =\int_{0}^{\infty}u(z,t)e^{-st}dt=U(z,s)
\]

eqs. (5) and (6) transform to%

\[
sU-(T_{0}-T_{w})=-v\frac{dU}{dz}-\alpha U
\]

or%

\begin{equation}
\frac{dU}{dz}=-(1/v)(s+\alpha)U+(1/v)(T_{0}-T_{w})\tag{8}%
\end{equation}

Eq. (8) is a linear, constant coefficient, nonohomogeneous ODE with
homogeneous and particular solutions%

\[
U_{h}=c_{1}e^{-(1/v)(s+\alpha)z}
\]%

\[
U_{p}=\frac{T_{0}-T_{w}}{s+\alpha}
\]

Thus%

\[
U(z,s)=c_{1}e^{-(1/v)(s+\alpha)z}+\frac{T_{0}-T_{w}}{s+\alpha}
\]

Application of BC (7) to evaluate constant $c_{1}$ gives%

\[
\frac{T_{e}-T_{w}}{s}=c_{1}+\frac{T_{0}-T_{w}}{s+\alpha}
\]

or%

\[
c_{1}=\frac{T_{e}-T_{w}}{s}-\frac{T_{0}-T_{w}}{s+\alpha}
\]

Thus%

\[
U(z,s)=\left(  \frac{T_{e}-T_{w}}{s}-\frac{T_{0}-T_{w}}{s+\alpha}\right)
e^{-(1/v)(s+\alpha)z}+\frac{T_{0}-T_{w}}{s+\alpha}
\]%

\begin{equation}
=e^{-(\alpha/v)z}\left(  \frac{T_{e}-T_{w}}{s}-\frac{T_{0}-T_{w}}{s+\alpha
}\right)  e^{-(z/v)s}+\frac{T_{0}-T_{w}}{s+\alpha}\tag{9}%
\end{equation}

Inversion of eq. (9) gives%

\[
u(z,t)=L_{t}^{-1}\left\{  U(z,s)\right\}
\]%

\[
=e^{-(\alpha/v)z}\left\{  (T_{e}-T_{w})-\left(  T_{0}-T_{w}\right)
e^{-\alpha\{t-(z/v)\}}\right\}  h(t-(z/v))+(T_{0}-T_{w})e^{-\alpha t}
\]

where%

\[%
\begin{array}
[c]{cc}%
h(t)=0 & t<0\\
h(t)=1 & t>0
\end{array}
\]

Thus%

\[
T(z,t)=T_{w}+u(z,t)
\]%

\begin{equation}
=T_{w}+e^{-(\alpha/v)z}\left\{  (T_{e}-T_{w})-\left(  T_{0}-T_{w}\right)
e^{-\alpha\{t-(z/v)\}}\right\}  h(t-(z/v))+(T_{0}-T_{w})e^{-\alpha t}\tag{10}%
\end{equation}

Check (being a little ``loose'' with $h(z,t)$ and it's derivative at $t=z/v$):%

\[
T(z,0)=T_{w}+0+(T_{0}-T_{w})e^{-0}=T_{0}
\]%

\[
T(0,t)=T_{w}+e^{0}\left\{  (T_{e}-T_{w})-\left(  T_{0}-T_{w}\right)
e^{-\alpha t}\right\}  +(T_{0}-T_{w})e^{-\alpha t}=T_{e}
\]
\[%
\begin{array}
[c]{cc}%
PDE\text{, eq. (1)} & \text{From eq. (10)}\\
&\\
\partial T/\partial t & e^{-(\alpha/v)z}\left\{  -\left(  T_{0}-T_{w}\right)
e^{-\alpha\{t-(z/v)\}}(-\alpha)\right\}  +(T_{0}-T_{w})(-\alpha)e^{-\alpha
t}\\
&\\
v\partial T/\partial z & v\left[  -(\alpha/v)e^{-(\alpha/v)z}(T_{e}%
-T_{w})\right] \\
&\\
-\alpha(T_{w}-T) & -\alpha T_{w}+\alpha[T_{w}+e^{-(\alpha/v)z}\left\{
(T_{e}-T_{w})-\left(  T_{0}-T_{w}\right)  e^{-\alpha\{t-(z/v)\}}\right\}
h(t-(z/v))\\
&\\
& +(T_{0}-T_{w})e^{-\alpha t}]\\
&\\
0 & 0
\end{array}
\]

Eq. (10) consists of two parts (as defined by $h(t-(z/v))$):%

\begin{equation}
t<z/v,T(z,t)=T_{w}+(T_{0}-T_{w})e^{-\alpha t}\tag{10a}%
\end{equation}%
\begin{equation}
t>z/v,T(z,t)=T_{w}+e^{-(\alpha/v)z}\left\{  (T_{e}-T_{w})-\left(  T_{0}%
-T_{w}\right)  e^{-\alpha\{t-(z/v)\}}\right\}  +(T_{0}-T_{w})e^{-\alpha
t}\tag{10b}%
\end{equation}

Eqs. (10a) and (10b) give a discontinuity at $t=z/v$. To see this, if we
consider eq. (10) at $t=z/v$,

\vspace{0.5cm}

(a) with $u(t-z/v)=0,$
\[
T(z,t)=T_{w}+e^{-(\alpha/v)z}\left\{  (T_{e}-T_{w})-\left(  T_{0}%
-T_{w}\right)  e^{-\alpha\{t-(z/v)\}}\right\}  0+(T_{0}-T_{w})e^{-\alpha t}
\]
\begin{equation}
=T_{w}+(T_{0}-T_{w})e^{-\alpha t}\tag{10c}%
\end{equation}

(b) with $u(t-z/v)=1$%

\[
T(z,t)=T_{w}+e^{-\alpha t}\left\{  (T_{e}-T_{w})-\left(  T_{0}-T_{w}\right)
e^{-\alpha\{0\}}\right\}  (1)+(T_{0}-T_{w})e^{-\alpha t}
\]
\begin{equation}
=T_{w}+(T_{e}-T_{0})e^{-\alpha t}+(T_{0}-T_{w})e^{-\alpha t}=T_{w}%
+(T_{e}-T_{w})e^{-\alpha t}\tag{10d}%
\end{equation}

Only when $T_{0}=T_{e}$ are eqs. (10c) and (10d) the same, i.e., the initial
and entering conditions are the same. Otherwise, with $T_{e}\neq T_{0}$, a
discontinuity is introduced that propagates through the system. \textit{This
is a basic property of hyperbolic PDEs, i.e., they propagate discontinuities,
which is the reason they are the most difficult of the three classes of PDEs
(elliptic, parabolic, hyperbolic) to solve numerically}; in particular, the
discontinuity often occurs through an incompatibility between the initial and
boundary conditions, e.g., $T_{e}\neq T_{0}$.

\vspace{0.5cm}

There are some other interesting special cases of eqs. (10a) and (10b)%

\begin{equation}%
\begin{array}
[c]{ccc}%
t=0 &  & \\
& t<z/v & T(z,t)=T_{0}\\
& t>z/v & \text{NA}%
\end{array}
\tag{11a}%
\end{equation}%
\begin{equation}%
\begin{array}
[c]{ccc}%
t\rightarrow\infty &  & \\
&  t<z/v & \text{NA}\\
& t>z/v & T(z,t)=T_{w}+e^{-(\alpha/v)z}(T_{e}-T_{w})
\end{array}
\tag{11b}%
\end{equation}

Within eqs. (11b),%

\begin{equation}%
\begin{array}
[c]{cccc}%
t\rightarrow\infty &  &  & \\
&  t>z/v & z=0 & T(z,t)=T_{e}\\
& t>z/v & z\rightarrow\infty &  T(z,t)=T_{w}%
\end{array}
\tag{11c}%
\end{equation}

Also,%

\begin{equation}%
\begin{array}
[c]{ccc}%
\alpha=0 &  & \\
& t<z/v & T(z,t)=T_{0}\\
& t>z/v & T(z,t)=T_{e}%
\end{array}
\tag{11d}%
\end{equation}

Eqs. (11d) in particular illustrates the difficulty in computing numerical
solutions to first order, hyperbolic PDEs, i.e., the discontinuity generated
at $z=0$ where $T(z,t)$ switches from $T_{0}$ to $T_{e}$, ) (again, first
order, hyperbolic PDEs propagate discontinuities).

\vspace{0.5cm}

To summarize the case $\alpha=0$, eqs. (1) to (3) reduce to%

\begin{equation}
\frac{\partial T}{\partial t}+v\frac{\partial T}{\partial z}=0\tag{12a}%
\end{equation}%

\begin{equation}
T(z,0)=T_{0}\tag{12b}%
\end{equation}%

\begin{equation}
T(0,t)=T_{e}\tag{12c}%
\end{equation}

for which the solution is%

\[
T(z,t)=T_{0},t<z/v
\]%

\[
T(z,t)=T_{e},t>z/v
\]

or%

\begin{equation}
T(z,t)=T_{0}+(T_{e}-T_{0})h(t-z/v)\tag{13}%
\end{equation}

If $T_{0}=0,T_{e}=1$, this reduces to%

\begin{equation}
T(z,t)=h(t-z/v)\tag{14}%
\end{equation}

where again $h(t)$ is the Heaviside unit step function. Note that eq. (14)
indicates the Heaviside unit step function propagates along the system with
velocity $v$.

\vspace{0.5cm}

Another approach to eq. (11b) is to start with eq. (1) at steady state%

\begin{equation}
\frac{\partial T}{\partial t}=-v\frac{\partial T}{\partial z}+\frac{4h}{D\rho
C_{p}}(T_{w}-T)=0\tag{15}%
\end{equation}

Since eq. (15) has only one independent variable, $z$, it is actually an ODE.
The corresponding auxiliary condition is eq. (3) at steady state%

\begin{equation}
T(0)=T_{e}\tag{16}%
\end{equation}

The solution to eq. (15) will be a homogeneous part plus a particular part%

\[
T_{h}(z)=c_{1}e^{-(\alpha/v)z},\alpha=\frac{4h}{D\rho C_{p}}
\]
\[
T_{p}(z)=T_{w}
\]

Thus
\[
T(z)=T_{w}+c_{1}e^{-(\alpha/v)z}
\]

Then from eq. (16),%

\[
T_{e}=T_{w}+c_{1}e^{-0}
\]

and the solution is%

\[
T(z)=T_{w}+(T_{e}-T_{w})e^{-(\alpha/v)z}
\]

which is eq. (11b).

\vspace{0.5cm}

We now consider:

\vspace{0.5cm}

\begin{itemize}
\item  How eqs. (1) to (3) can be programmed

\item  The comparison of this numerical solution with the analytical solution,
eq. (10).
\end{itemize}

\vspace{0.5cm}

In order to program a PDE, we generally start by defining a spatial grid, in
this case, a grid in $z$. An index, $i$, is used to denote a particular
position along the grid. The spatial derivative in eq. (1), $\dfrac{\partial
T}{\partial z},$ is then approximated algebraically, in this case by a finite
difference, at each of the grid points. For example, at the general grid point $i$%

\[
\dfrac{\partial T_{i}}{\partial z}\cong\frac{T_{i}-T_{i-1}}{\Delta z}
\]
This leads to a system of ODEs, one at each grid point, that approximate the
PDE of eq. (1)%

\begin{equation}
\dfrac{dT_{i}}{dz}=-v(\frac{T_{i}-T_{i-1}}{\Delta z})+\frac{\pi Dh}{(\pi
D^{2}/4)\rho C_{p}}(T_{w}-T_{i}),i=1,2,\cdots,n\tag{17}%
\end{equation}

(note that in deriving eq. (17), we have essentiually reversed the process of
deriving eq. (1) starting with an incremental heat balance). This system of
simultanenous ODEs can then be integrated by any of the ODE\ integration
methods we discussed previously, e.g., Euler's method, the modified Euler
method, RKF45, BDF, etc. In the following program we use Euler's method to
provide coding that is easy to understand (but, of course, Euler's method is
only first order correct, so we might eventually want to switch to a more
acccurate ODE integrator, perhaps with automatic step size adjustment).

\vspace{0.5cm}

The details of the procedure are most easily understood by considering the
actual computer programming, as illustrated by the following program.

\vspace{0.5cm}

\textbf{Time and Spatial Discretization}

\vspace{0.5cm}

We now have to discretize the time derivative and the spatial derivative in
eq. (1). A Matlab program that does this is listed below (this program
illustrates all of the essential features for a complete PDE solution):
\begin{verbatim}

%   Fluid flow through a heated pipe
% 
%   The PDE that models the heated pipe is
%
%      T = -v*T  + (4*hc)/(D*rho*Cp)*(Tw - T)                (1)
%       t      z
%
%   with the initial and boundary conditions
%
%      T(z,0) = T0, T(0,t) = Te                          (2)(3)
%
%   where
%
%      T     fluid temperature
%
%      t     time
%
%      z     position along the pipe
%
%      v     fluid velocity
%
%      hc    fluid to pipe wall heat transfer coefficient
%
%      D     pipe diameter
%
%      rho   fluid density
%
%      Cp    fluid specific heat
%
%      T0    initial fluid temperature
%
%      Te    entering fluid temperature
%
%      Tw    pipe wall temperature
%
%      L     pipe length
%
%   Of particular interest is the fluid exiting temperature, T(L,t).
%
%   The method of lines (MOL) solution of eq. (1) is coded below.
%   Specifically, the spatial derivative in eq. (1) is replaced
%   by a three point, finite difference (FD) approximation, both
%   centered and upwind.  The resulting system of ODEs in t,
%   defined on a 21-point grid in z, is then integrated by the
%   Euler method.
%
%   The analytical solution to eqs. (1) to (3) is also programmed
%   for comparison with the numerical solution.
%
%   The following code is run for the four special cases:
%
%      Case 1: hc = 2, T0 = 25, Te = 25, Tw = 100
%
%      Case 2: hc = 0.2, T0 = 25, Te = 25, Tw = 100
%
%      Case 3: hc = 0, T0 = 25, Te = 25, Tw (NA)
%
%      Case 4: hc = 0, T0 = 25, Te = 100, TW (NA)
%
% Open output files
  fid1=fopen('htex1.out','w');
  fid2=fopen('htex1.plt','w');
%
% Model parameters (defined in the comments above)
  cp=1.0;
  rho=1.0;
  v=50.0;
  d=1.0;
  zl=100.0;
  n=21;
  dz=zl/(n-1);
%
% Select three point, finite differencing (FD) of spatial
% derivative
%
%   ifd = 1: Centered approximations
%
%   ifd = 2: Two point upwind approximation
%
%   ifd = 3: Five point, biased upwind approximation
%
%   ifd = 4: van Leer flux limiter
%
  ifd=1;
%
% Cases listed above
  for ncase=1:4
%
% Initial, final times, integration interval, number of Euler
% steps for each output
  t=0.0; tf=5.0; h=0.001; nout=250;
  if ncase==1
     hc=2.0;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==2
     hc=0.2;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==3
     hc=0.0;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==4
     hc=0.0;
     T0=25.0;
     Te=100.0;
     Tw=100.0;
  end
%
% Initial conditions
  for i=1:n
    T(i)=T0;
  end
%
% Write h
  fprintf(fid1,'\n\n ncase = %5d    h = %10.3e\n\n',ncase,h);
%
% Write heading
  fprintf(fid1,'   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff\n');
%
% Integrate until t = tf
  while t < tf*1.0001
%
% Monitor solution by displaying t
  t
%
% Write selected output, including analytical solution and difference
% between numerical and analytical solutions
  alpha=4.0*hc/(d*rho*cp);
  Ta=Tw+(T0-Tw)*exp(-alpha*t);
  if(t > zl/v)Ta=Ta+exp(-alpha/v*zl)*((Te-Tw)-(T0-Tw)*...
                      exp(-alpha*(t-zl/v))); end
  diff=T(n)-Ta;
  fprintf(fid1,'%5.2f%10.2f%10.2f%10.2f%10.3f\n',...
          t,t-zl/v,T(n),Ta,diff);
  fprintf(fid2,'%5.2f%10.2f%10.2f%10.2f%10.3f\n',...
          t,t-zl/v,T(n),Ta,diff);
%
% Take nout Euler steps
  for iout=1:nout
%
%   Temporal derivatives
%
%     Centered approximations
      if(ifd==1)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss002(0.0,zl,n,T);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of three point centered approximation
      end
%
%     Two point upwind approximation
      if(ifd==2)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss012(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of two point upwind approximation
      end
%
%     Five point, biased upwind approximation
      if(ifd==3)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss020(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of five point, biased upwind approximation
      end
%
%     van Leer flux limiTer
      if(ifd==4)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=vanl2(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of van Leer flux limiter
      end
%
%   Take Euler step
    for i=1:n
      T(i)=T(i)+Tt(i)*h;
    end
    t=t+h;
%
% Next Euler step
  end
%
% Next output
  end
%
% Next case
  end

\end{verbatim}

We can note the following points about this program:

\vspace{0.5cm}

(1) The model parameters are first set
\begin{verbatim}

% Model parameters (defined in the comments above)
  cp=1.0;
  rho=1.0;
  v=50.0;
  d=1.0;
  zl=100.0;
  n=21;
  dz=zl/(n-1);

\end{verbatim}

(2) A spatial differentiator for the derivative $T_{z}$ is selected
\begin{verbatim}

% Select finite differencing (FD) of spatial derivative
%
%   ifd = 1: Centered approximations
%
%   ifd = 2: Two point upwind approximation
%
%   ifd = 3: Five point, biased upwind approximation
%
%   ifd = 4: van Leer flux limiter
%
  ifd=1;

\end{verbatim}

The details of these differentiators will be considered below.

\vspace{0.5cm}

(3) The integration parameters and the specific values of the model parameters
for one of four cases are selected
\begin{verbatim}

%   The following code is run for the four special cases:
%
%      Case 1: hc = 2, T0 = 25, Te = 25, Tw = 100
%
%      Case 2: hc = 0.2, T0 = 25, Te = 25, Tw = 100
%
%      Case 3: hc = 0, T0 = 25, Te = 25, Tw (NA)
%
%      Case 4: hc = 0, T0 = 25, Te = 100, TW (NA)

% Cases listed above
  for ncase=1:4
%
% Initial, final times, integration interval, number of Euler
% steps for each output
  t=0.0; tf=5.0; h=0.001; nout=250;
  if ncase==1
     hc=2.0;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==2
     hc=0.2;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==3
     hc=0.0;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==4
     hc=0.0;
     T0=25.0;
     Te=100.0;
     Tw=100.0;
  end

\end{verbatim}

If the integration step $h$ in the Euler integration is increased above $h =
0.001$, the stability limit of the Euler method is exceeded, i.e.,
$|h\lambda|< 2$, and the calculation goes unstable.

\vspace{0.5cm}

(4) The 21 initial conditions are set (for the 21-point grid in z)
\begin{verbatim}

% Initial conditions
  for i=1:n
    T(i)=T0;
  end

\end{verbatim}

(5) A heading for the output is next
\begin{verbatim}

%
% Write h
  fprintf(fid1,'\n\n ncase = %5d    h = %10.3e\n\n',ncase,h);
%
% Write heading
  fprintf(fid1,'   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff\n');

\end{verbatim}

Note that both the numerical and analytical solutions, and their difference
are displayed in fid1.

\vspace{0.5cm}

(6) The time stepping begins at t = 0 and ends at t = tf
\begin{verbatim}

% Integrate until t = tf
  while t < tf*1.0001
%
% Monitor solution by displaying t
  t

\end{verbatim}

(7) The analytical solution. eq. (10), is evaluated, and the numerical and
analytical solutions and their differenerence are then written to files, fid1
and fid2
\begin{verbatim}

% Write selected output, including analytical solution and difference
% between numerical and analytical solutions
  alpha=4.0*hc/(d*rho*cp);
  Ta=Tw+(T0-Tw)*exp(-alpha*t);
  if(t > zl/v)Ta=Ta+exp(-alpha/v*zl)*((Te-Tw)-(T0-Tw)*...
                      exp(-alpha*(t-zl/v))); end
  diff=T(n)-Ta;
  fprintf(fid1,'%5.2f%10.2f%10.2f%10.2f%10.3f\n',...
          t,t-zl/v,T(n),Ta,diff);
  fprintf(fid2,'%5.2f%10.2f%10.2f%10.2f%10.3f\n',...
          t,t-zl/v,T(n),Ta,diff);

\end{verbatim}

(8) nout (= 250) Euler steps are taken for the n (= 21) ODEs in $t$
\begin{verbatim}

% Take nout Euler steps
  for iout=1:nout

\end{verbatim}

(9) The calculation of the time derivatives of eq. (1), $T_{t}$, requires
first the evaluation of the spatial derivative, $T_{z}$. This is done by one
of the four spatial differentiators selected by ifd (see (2) above). The RHS
of eq (1) is then computed in each case
\begin{verbatim}

%   Temporal derivatives
%
%     Centered approximations
      if(ifd==1)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss002(0.0,zl,n,T);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of three point centered approximation
      end
%
%     Two point upwind approximation
      if(ifd==2)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss012(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of two point upwind approximation
      end
%
%     Five point, biased upwind approximation
      if(ifd==3)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss020(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of five point, biased upwind approximation
      end
%
%     van Leer flux limiter
      if(ifd==4)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=vanl2(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of van Leer flux limiter
      end

\end{verbatim}

(10) Once the n (= 21) time derivatives are computed, the n ODEs are
integrated by the Euler method to move the solution along (as parallel lines
in $z$)
\begin{verbatim}

%   Take Euler step
    for i=1:n
      T(i)=T(i)+Tt(i)*h;
    end
    t=t+h;

\end{verbatim}

(11) The loop for the nout (= 250) Euler steps is then completed and the
numerical and analytical solutions are printed before the next set of nout
Euler steps is taken (see (7) above). First, the test for the final time, tf,
is checked (see (6) above).
\begin{verbatim}

% Next Euler step
  end
%
% Next output
  end
%
% Next case
  end

\end{verbatim}

The preceding Matlab code is in

\vspace{0.5cm}

\begin{center}
http://www.lehigh.edu/~wes1/apci/htex1.m
\end{center}

\vspace{0.5cm}

The extension of this code to systems of ODEs/PDEs is straightforward; their
time derivatives would be added, and the corresponding ODEs would be
integrated, either by the Euler method, or by a more accurate method, e.g.,

\begin{itemize}
\item  The modified Euler method, the classical Runge Kutta method, the Runge
Kutta Fehlberg method, either fixed step or variable step

\item  A library integrator called where the Euler integration is programmed,
e.g., a call to RKF45, LSODE, LSODES, DASSL, a NAG ODE integrator, a Matlab
ODE integrator (ode23, ode45)

\vspace{0.5cm}
\end{itemize}

This modular and flexible approach to programming PDEs is termed the
\textit{method of lines} (MOL). Additionally, the algorithm for calculating
the spatial derivative(s) can easily be changed (by calling another routine in
place of one of the library differentiators).

\vspace{0.5cm}

The output of the preceding program for ifd = 1 (three point centered FD) is
listed below (and is plotted in htex1d1.ps)
\begin{verbatim}

 ncase =     1    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     89.93     89.85     0.081
 0.50     -1.50     98.65     98.63     0.022
 0.75     -1.25     99.82     99.81     0.004
 1.00     -1.00     99.98     99.97     0.001
 1.25     -0.75    100.00    100.00     0.000
 1.50     -0.50    100.00    100.00     0.000
 1.75     -0.25    100.00    100.00     0.000
 2.00      0.00    100.00    100.00     0.000
 2.25      0.25    100.00    100.00     0.000
 2.50      0.50    100.00    100.00     0.000
 2.75      0.75    100.00    100.00     0.000
 3.00      1.00    100.00    100.00     0.000
 3.25      1.25    100.00    100.00     0.000
 3.50      1.50    100.00    100.00     0.000
 3.75      1.75    100.00    100.00     0.000
 4.00      2.00    100.00    100.00     0.000
 4.25      2.25    100.00    100.00     0.000
 4.50      2.50    100.00    100.00     0.000
 4.75      2.75    100.00    100.00     0.000
 5.00      3.00    100.00    100.00     0.000


 ncase =     2    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     38.60     38.60     0.005
 0.50     -1.50     49.73     49.73     0.008
 0.75     -1.25     58.85     58.84     0.010
 1.00     -1.00     66.31     66.30     0.011
 1.25     -0.75     72.42     72.41     0.011
 1.50     -0.50     77.41     77.41    -0.003
 1.75     -0.25     81.37     81.51    -0.134
 2.00      0.00     84.10     84.86    -0.756
 2.25      0.25     85.22     84.86     0.363
 2.50      0.50     85.03     84.86     0.168
 2.75      0.75     84.69     84.86    -0.164
 3.00      1.00     84.81     84.86    -0.047
 3.25      1.25     84.90     84.86     0.039
 3.50      1.50     84.81     84.86    -0.052
 3.75      1.75     84.82     84.86    -0.035
 4.00      2.00     84.85     84.86    -0.005
 4.25      2.25     84.82     84.86    -0.036
 4.50      2.50     84.83     84.86    -0.025
 4.75      2.75     84.84     84.86    -0.019
 5.00      3.00     84.83     84.86    -0.030


 ncase =     3    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.00     25.00     0.000
 1.00     -1.00     25.00     25.00     0.000
 1.25     -0.75     25.00     25.00     0.000
 1.50     -0.50     25.00     25.00     0.000
 1.75     -0.25     25.00     25.00     0.000
 2.00      0.00     25.00     25.00     0.000
 2.25      0.25     25.00     25.00     0.000
 2.50      0.50     25.00     25.00     0.000
 2.75      0.75     25.00     25.00     0.000
 3.00      1.00     25.00     25.00     0.000
 3.25      1.25     25.00     25.00     0.000
 3.50      1.50     25.00     25.00     0.000
 3.75      1.75     25.00     25.00     0.000
 4.00      2.00     25.00     25.00     0.000
 4.25      2.25     25.00     25.00     0.000
 4.50      2.50     25.00     25.00     0.000
 4.75      2.75     25.00     25.00     0.000
 5.00      3.00     25.00     25.00     0.000


 ncase =     4    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.00     25.00     0.000
 1.00     -1.00     25.00     25.00     0.000
 1.25     -0.75     25.03     25.00     0.027
 1.50     -0.50     25.63     25.00     0.627
 1.75     -0.25     30.73     25.00     5.733
 2.00      0.00     51.17     25.00    26.166
 2.25      0.25     90.33    100.00    -9.665
 2.50      0.50    117.36    100.00    17.363
 2.75      0.75    103.95    100.00     3.951
 3.00      1.00     88.96    100.00   -11.039
 3.25      1.25    103.31    100.00     3.309
 3.50      1.50    104.99    100.00     4.991
 3.75      1.75     93.70    100.00    -6.300
 4.00      2.00    102.14    100.00     2.136
 4.25      2.25    102.54    100.00     2.540
 4.50      2.50     95.71    100.00    -4.288
 4.75      2.75    102.78    100.00     2.784
 5.00      3.00    100.18    100.00     0.176

\end{verbatim}

We can note the following points about this output:

\vspace{0.5cm}

(1) The three point centered FD approximations give good accuracy for Cases 1
to 3 which are for relatively smooth problems. In particular, the FD solution
is exact for Case 3 (the constant solution is differentiated exactly).

\vspace{0.5cm}

(2) The centered FD approximations do not give an accurate solution for the
discontinuous solution of Case 4. This result indicates that \textit{centered
approximations should not be used for nonsmooth solutions of first order
hyperbolic PDEs}.

\vspace{0.5cm}

The library differentiator for the three point centered FD approximations is
listed below
\begin{verbatim}

      function [ux]=dss002(xl,xu,n,u)
%...
%...  FUNCTION DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A
%...                                                  X
%...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU
%...
%...  ARGUMENT LIST
%...
%...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
%...
%...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
%...
%...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
%...             BOUNDARY POINTS (INPUT)
%...
%...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
%...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
%...             TO BE COMPUTED (INPUT)
%...
%...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
%...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
%...             (OUTPUT)
%...
%...  THE WEIGHTING COEFFICIENTS CAN BE SUMMARIZED AS
%...
%...          -3   4  -1
%...
%...     1/2  -1   0   1
%...
%...           1  -4   3
%...
%...  WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 2, M =
%...  1, P = 0, 1, 2 (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFER-
%...  ENTIATION, MATH. GAZ., VOL. 25, 1941).
%...
%...  COMPUTE THE SPATIAL INCREMENT
      dx=(xu-xl)/(n-1);
      r2fdx=1./(2.*dx);
      nm1=n-1;
%...
%...  LEFT END POINT.  THE FOLLOWING CODING HAS BEEN FORMATTED
%     SO THAT THE NUMERICAL WEIGHTING COEFFICIENTS CAN BE MORE
%...  EASILY ASSOCIATED WITH THE BICKLEY MATRIX LISTED ABOVE)
      ux(1)=r2fdx*...
      (     -3.   *u(  1)     +4.   *u(  2)     -1.   *u(  3));
%...
%...  INTERIOR POINTS
      for i=2:nm1
      ux(i)=r2fdx*...
      (     -1.   *u(i-1)     +0.   *u(  i)     +1.   *u(i+1));
      end
%...
%...  RIGHT END POINT
      ux(n)=r2fdx*...
      (      1.   *u(n-2)     -4.   *u(n-1)     +3.   *u(  n));

\end{verbatim}

The derivation of the FD approximations will be considered in the next section.

\vspace{0.5cm}

The preceding output for Case = 4 illustrates the first type of numerical
distortion that is common in the numerical integration of first order
hyperbolic PDEs, that is, \textit{numerical oscillation}. In order to
eliminate this oscillation, we consider two point upwind (noncentered)
approximations, as implemented with ifd = 2. The output of the preceding
program for ifd = 2 is listed below (and is plotted in htex1d2.ps)
\begin{verbatim}

 ncase =     1    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     89.93     89.85     0.081
 0.50     -1.50     98.65     98.63     0.022
 0.75     -1.25     99.82     99.81     0.004
 1.00     -1.00     99.98     99.97     0.001
 1.25     -0.75    100.00    100.00     0.000
 1.50     -0.50    100.00    100.00     0.000
 1.75     -0.25    100.00    100.00    -0.001
 2.00      0.00    100.00    100.00    -0.001
 2.25      0.25    100.00    100.00    -0.001
 2.50      0.50    100.00    100.00    -0.001
 2.75      0.75    100.00    100.00    -0.001
 3.00      1.00    100.00    100.00    -0.001
 3.25      1.25    100.00    100.00    -0.001
 3.50      1.50    100.00    100.00    -0.001
 3.75      1.75    100.00    100.00    -0.001
 4.00      2.00    100.00    100.00    -0.001
 4.25      2.25    100.00    100.00    -0.001
 4.50      2.50    100.00    100.00    -0.001
 4.75      2.75    100.00    100.00    -0.001
 5.00      3.00    100.00    100.00    -0.001


 ncase =     2    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     38.60     38.60     0.005
 0.50     -1.50     49.73     49.73     0.008
 0.75     -1.25     58.85     58.84     0.010
 1.00     -1.00     66.30     66.30     0.003
 1.25     -0.75     72.34     72.41    -0.073
 1.50     -0.50     77.00     77.41    -0.414
 1.75     -0.25     80.25     81.51    -1.256
 2.00      0.00     82.21     84.86    -2.644
 2.25      0.25     83.22     84.86    -1.635
 2.50      0.50     83.66     84.86    -1.194
 2.75      0.75     83.83     84.86    -1.026
 3.00      1.00     83.89     84.86    -0.971
 3.25      1.25     83.90     84.86    -0.955
 3.50      1.50     83.91     84.86    -0.950
 3.75      1.75     83.91     84.86    -0.949
 4.00      2.00     83.91     84.86    -0.949
 4.25      2.25     83.91     84.86    -0.949
 4.50      2.50     83.91     84.86    -0.949
 4.75      2.75     83.91     84.86    -0.949
 5.00      3.00     83.91     84.86    -0.949


 ncase =     3    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.00     25.00     0.000
 1.00     -1.00     25.00     25.00     0.000
 1.25     -0.75     25.00     25.00     0.000
 1.50     -0.50     25.00     25.00     0.000
 1.75     -0.25     25.00     25.00     0.000
 2.00      0.00     25.00     25.00     0.000
 2.25      0.25     25.00     25.00     0.000
 2.50      0.50     25.00     25.00     0.000
 2.75      0.75     25.00     25.00     0.000
 3.00      1.00     25.00     25.00     0.000
 3.25      1.25     25.00     25.00     0.000
 3.50      1.50     25.00     25.00     0.000
 3.75      1.75     25.00     25.00     0.000
 4.00      2.00     25.00     25.00     0.000
 4.25      2.25     25.00     25.00     0.000
 4.50      2.50     25.00     25.00     0.000
 4.75      2.75     25.00     25.00     0.000
 5.00      3.00     25.00     25.00     0.000


 ncase =     4    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.01     25.00     0.008
 1.00     -1.00     25.25     25.00     0.247
 1.25     -0.75     27.24     25.00     2.243
 1.50     -0.50     34.27     25.00     9.275
 1.75     -0.25     47.86     25.00    22.861
 2.00      0.00     64.76     25.00    39.764
 2.25      0.25     79.80    100.00   -20.201
 2.50      0.50     90.08    100.00    -9.924
 2.75      0.75     95.76    100.00    -4.243
 3.00      1.00     98.40    100.00    -1.604
 3.25      1.25     99.46    100.00    -0.544
 3.50      1.50     99.83    100.00    -0.168
 3.75      1.75     99.95    100.00    -0.047
 4.00      2.00     99.99    100.00    -0.012
 4.25      2.25    100.00    100.00    -0.003
 4.50      2.50    100.00    100.00    -0.001
 4.75      2.75    100.00    100.00     0.000
 5.00      3.00    100.00    100.00     0.000

\end{verbatim}

We can note the following points about this output:

\vspace{0.5cm}

(1) The two point point FD approximations give good accuracy for Cases 1 to 3
which again are for relatively smooth problems. In particular, the FD solution
is exact for Case 3 (the constant solution is differentiated exactly).

\vspace{0.5cm}

(2) The upwind FD approximations do not give an accurate solution for the
discontinuous solution of Case 4. Although the previous numerical oscillation
has been eliminated, there is now significant \textit{numerical diffusion},
which is the second common form of numerical distortion in the solution of
first order hyperbolic PDEs.

\vspace{0.5cm}

The library differentiator for the two point upwind FD approximations is
listed below. Note that the direction of the flow is required to use upwinding
(it is the sixth input argument of the routine, v). This is a consequence of
the ``flow of information'' in the solution that is taken into account by the
approximation, i.e., upwind points have a greater influence on the solution
than downwind points. This characteristic indicates an important feature of
hyperbolic PDEs, that is, they have a preferred direction; physically, this
preferred direction in convective systems is the direction of flow.

\vspace{0.5cm}

Experience has indicated that some form of upwinding (to take into account the
preferred direction) is essential for strongly hyperbolic systems. In other
words, centered approximations for hyperbolic systems generally don't work
(consider the ifd = 1 results discussed previously). In order to use
upwinding, the direction of flow must be known (which may not be easy to
determine in complex, multidimensional flow systems). Upwinding in the wrong
direction (``downwinding'') usually leads to unstable solutions.
\begin{verbatim}

      function [ux]=dss012(xl,xu,n,u,v)
%...
%...  FUNCTION DSS012 IS AN APPLICATION OF FIRST-ORDER DIRECTIONAL
%...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
%...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
%...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS WITH THE
%...  SIMPLEST FORM
%...
%...                            U  + V*U  = 0                        (1)
%...                             T      X
%...
%...  THE FIRST FOUR PARAMETERS, XL, XU, N AND U, ARE THE SAME AS
%...  FOR FUNCTION DSS002.  THE FIFTH PARAMETER, V, MUST BE PROVIDED
%...  TO DSS012 SO THAT THE DIRECTION OF FLOW IN EQUATION (1) CAN BE
%...  USED TO SELECT THE APPROPRIATE FINITE DIFFERENCE APPROXIMATION
%...  FOR THE FIRST-ORDER SPATIAL DERIVATIVE IN EQUATION (1), U .
%...  THE CONVENTION FOR THE SIGN OF V IS                      X
%...
%...     FLOW LEFT TO RIGHT             V GT 0
%...     (I.E., IN THE DIRECTION        (I.E., THE SIXTH ARGUMENT IS
%...     OF INCREASING X)               POSITIVE IN CALLING DSS012)
%...
%...     FLOW RIGHT TO LEFT             V LT 0
%...     (I.E., IN THE DIRECTION        (I.E., THE SIXTH ARGUMENT IS
%...     OF DECREASING X)               NEGATIVE IN CALLING DSS012)
%...
%...  COMPUTE THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFERENCE
%...  APPROXIMATION DEPENDING ON THE SIGN OF V IN EQUATION (1).
      dx=(xu-xl)/(n-1);
      if v > 0
%...
%...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
              ux(1)=(u(2)-u(1))/dx;
              for i=2:n
                 ux(i)=(u(i)-u(i-1))/dx;
              end
      end
%...
%...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
      if v < 0
              nm1=n-1;
              for i=1:nm1
                 ux(i)=(u(i+1)-u(i))/dx;
              end
              ux(n)=(u(n)-u(n-1))/dx;
      end

\end{verbatim}

The origin of the FD approximation is evident (it can also be considered a
``stirred-tanks-in-series'' approximation).

\vspace{0.5cm}

The preceding two outputs for Case = 4 illustrates the two types of numerical
error that are common in the numerical integration of first order hyperbolic
PDEs, that is, \textit{numerical oscillation} and \textit{numerical
diffusion}. In order to reduce this error, we can consider combining the two
approaches, i.e., using a biased upwind approximation (which is not fully
centered or fully upwind, but a combination of the two). This is done by
executing the preceding program with ifd = 3 for a five point biased upwind
(5pbu) approximation coded in function dss020 listed subsequently (note that
in calling the 5pbu routine, the direction of flow is required since the
approximations use upwinding); the output is listed below (and is plotted in htex1d3.ps)
\begin{verbatim}

 ncase =     1    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     89.93     89.85     0.081
 0.50     -1.50     98.65     98.63     0.022
 0.75     -1.25     99.82     99.81     0.004
 1.00     -1.00     99.98     99.97     0.001
 1.25     -0.75    100.00    100.00     0.000
 1.50     -0.50    100.00    100.00     0.000
 1.75     -0.25    100.00    100.00     0.000
 2.00      0.00    100.00    100.00     0.000
 2.25      0.25    100.00    100.00     0.000
 2.50      0.50    100.00    100.00     0.000
 2.75      0.75    100.00    100.00     0.000
 3.00      1.00    100.00    100.00     0.000
 3.25      1.25    100.00    100.00     0.000
 3.50      1.50    100.00    100.00     0.000
 3.75      1.75    100.00    100.00     0.000
 4.00      2.00    100.00    100.00     0.000
 4.25      2.25    100.00    100.00     0.000
 4.50      2.50    100.00    100.00     0.000
 4.75      2.75    100.00    100.00     0.000
 5.00      3.00    100.00    100.00     0.000


 ncase =     2    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     38.60     38.60     0.005
 0.50     -1.50     49.73     49.73     0.008
 0.75     -1.25     58.85     58.84     0.010
 1.00     -1.00     66.31     66.30     0.007
 1.25     -0.75     72.44     72.41     0.036
 1.50     -0.50     77.34     77.41    -0.075
 1.75     -0.25     81.65     81.51     0.142
 2.00      0.00     84.58     84.86    -0.279
 2.25      0.25     84.87     84.86     0.008
 2.50      0.50     84.86     84.86     0.004
 2.75      0.75     84.86     84.86     0.003
 3.00      1.00     84.86     84.86    -0.003
 3.25      1.25     84.86     84.86     0.002
 3.50      1.50     84.86     84.86    -0.001
 3.75      1.75     84.86     84.86     0.001
 4.00      2.00     84.86     84.86     0.000
 4.25      2.25     84.86     84.86     0.000
 4.50      2.50     84.86     84.86     0.000
 4.75      2.75     84.86     84.86     0.000
 5.00      3.00     84.86     84.86     0.000


 ncase =     3    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.00     25.00     0.000
 1.00     -1.00     25.00     25.00     0.000
 1.25     -0.75     25.00     25.00     0.000
 1.50     -0.50     25.00     25.00     0.000
 1.75     -0.25     25.00     25.00     0.000
 2.00      0.00     25.00     25.00     0.000
 2.25      0.25     25.00     25.00     0.000
 2.50      0.50     25.00     25.00     0.000
 2.75      0.75     25.00     25.00     0.000
 3.00      1.00     25.00     25.00     0.000
 3.25      1.25     25.00     25.00     0.000
 3.50      1.50     25.00     25.00     0.000
 3.75      1.75     25.00     25.00     0.000
 4.00      2.00     25.00     25.00     0.000
 4.25      2.25     25.00     25.00     0.000
 4.50      2.50     25.00     25.00     0.000
 4.75      2.75     25.00     25.00     0.000
 5.00      3.00     25.00     25.00     0.000


 ncase =     4    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00    -0.003
 0.75     -1.25     24.97     25.00    -0.031
 1.00     -1.00     25.26     25.00     0.256
 1.25     -0.75     24.19     25.00    -0.809
 1.50     -0.50     28.03     25.00     3.026
 1.75     -0.25     15.27     25.00    -9.732
 2.00      0.00     67.21     25.00    42.208
 2.25      0.25    103.99    100.00     3.991
 2.50      0.50     97.68    100.00    -2.315
 2.75      0.75    101.67    100.00     1.673
 3.00      1.00     99.05    100.00    -0.946
 3.25      1.25    100.51    100.00     0.514
 3.50      1.50     99.74    100.00    -0.260
 3.75      1.75    100.12    100.00     0.119
 4.00      2.00     99.95    100.00    -0.045
 4.25      2.25    100.01    100.00     0.010
 4.50      2.50    100.00    100.00     0.004
 4.75      2.75     99.99    100.00    -0.008
 5.00      3.00    100.01    100.00     0.008

\end{verbatim}

We can note the following points about this output:

\vspace{0.5cm}

(1) The five point biased upwind (5pbu) FD approximations give good accuracy
for Cases 1 to 3 which again are relatively smooth problems. In particular,
the FD solution is exact for Case 3 (the constant solution is differentiated exactly).

\vspace{0.5cm}

(2) The 5pbu FD approximations give an improvement over centered
approximations (ifd = 1) and low order upwinding (ifd = 2) by reducing both
the numerical diffusion and oscillation, but these numerical errors still
persist. This is a consequence of Godunov's theorem to be discussed
subsequently. Briefly, we have probably done about as well as we can do with
linear (FD) approximations, i.e., computing approximations to derivatives
using linear combinations of the dependent variable, for example
\begin{verbatim}

      ux(  i)=r4fdx*...
      (  -1.     *u(i-3)...
         +6.     *u(i-2)...
        -18.     *u(i-1)...
        +10.     *u(i  )...
         +3.     *u(i+1));

\end{verbatim}

We will therefore next turn to a nonlinear FD approximation as the final case
(ifd = 4).

\vspace{0.5cm}

The library differentiator for the 5pbu FD approximations (ifd = 3) is listed below
\begin{verbatim}

      function [ux]=dss020(xl,xu,n,u,v)
%...
%...  SUBROUTINE DSS020 IS AN APPLICATION OF FOURTH-ORDER DIRECTIONAL
%...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
%...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
%...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS-
%...  CUSSED IN SUBROUTINE DSS012.  THE COEFFICIENTS OF THE FINITE
%...  DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W.
%...  G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL
%...  GAZETTE, PP. 19-27, 1941, N = 4, M = 1, P = 0, 1, 2, 3, 4.  THE
%...  IMPLEMENTATION IS THE **FIVE-POINT BIASED UPWIND FORMULA** OF
%...  M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE
%...  ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69,
%...  AUGUST, 1978
%...
%...  COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION
%...  CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER-
%...  ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT).
      dx=(xu-xl)/(n-1);
      r4fdx=1./(12.*dx);
%...
%...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
      if v > 0.
      ux(  1)=r4fdx*...
      ( -25.     *u(  1)...
        +48.     *u(  2)...
        -36.     *u(  3)...
        +16.     *u(  4)...
         -3.     *u(  5));
      ux(  2)=r4fdx*...
      (  -3.     *u(  1)...
        -10.     *u(  2)...
        +18.     *u(  3)...
         -6.     *u(  4)...
         +1.     *u(  5));
      ux(  3)=r4fdx*...
      (  +1.     *u(  1)...
         -8.     *u(  2)...
         +0.     *u(  3)...
         +8.     *u(  4)...
         -1.     *u(  5));
      nm1=n-1;
      for i=4:nm1
      ux(  i)=r4fdx*...
      (  -1.     *u(i-3)...
         +6.     *u(i-2)...
        -18.     *u(i-1)...
        +10.     *u(i  )...
         +3.     *u(i+1));
      end
      ux(  n)=r4fdx*...
      (  +3.     *u(n-4)...
        -16.     *u(n-3)...
        +36.     *u(n-2)...
        -48.     *u(n-1)...
        +25.     *u(n  ));
      end
%...
%...  (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
      if v < 0.
      ux(  1)=r4fdx*...
      ( -25.     *u(  1)...
        +48.     *u(  2)...
        -36.     *u(  3)...
        +16.     *u(  4)...
         -3.     *u(  5));
      nm3=n-3;
      for i=2:nm3
      ux(  i)=r4fdx*...
      (  -3.     *u(i-1)...
        -10.     *u(i  )...
        +18.     *u(i+1)...
         -6.     *u(i+2)...
         +1.     *u(i+3));
      end
      ux(n-2)=r4fdx*...
      (  +1.     *u(n-4)...
         -8.     *u(n-3)...
         +0.     *u(n-2)...
         +8.     *u(n-1)...
         -1.     *u(n  ));
      ux(n-1)=r4fdx*...
      (  -1.     *u(n-4)...
         +6.     *u(n-3)...
        -18.     *u(n-2)...
        +10.     *u(n-1)...
         +3.     *u(n  ));
      ux(  n)=r4fdx*...
      (  +3.     *u(n-4)...
        -16.     *u(n-3)...
        +36.     *u(n-2)...
        -48.     *u(n-1)...
        +25.     *u(n  ));
      end

\end{verbatim}

The 5pbu FD approximations are derived from the Taylor series as discussed in
the next section.

\vspace{0.5cm}

To conclude this discussion of time and spatial discretization, we consider
the use of the van Leer flux limiter to compute spatial derivatives (by
setting ifd = 4 in the preceding program). The objective is to eliminate
numerical diffusion and oscillation as much as possible (in the solution of
first order, hyperbolic PDEs). The output is first listed (and plotted in htex1d4.ps)
\begin{verbatim}

 ncase =     1    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     89.93     89.85     0.081
 0.50     -1.50     98.65     98.63     0.022
 0.75     -1.25     99.82     99.81     0.004
 1.00     -1.00     99.98     99.97     0.001
 1.25     -0.75    100.00    100.00     0.000
 1.50     -0.50    100.00    100.00     0.000
 1.75     -0.25    100.00    100.00     0.000
 2.00      0.00    100.00    100.00     0.000
 2.25      0.25    100.00    100.00     0.000
 2.50      0.50    100.00    100.00     0.000
 2.75      0.75    100.00    100.00     0.000
 3.00      1.00    100.00    100.00     0.000
 3.25      1.25    100.00    100.00     0.000
 3.50      1.50    100.00    100.00     0.000
 3.75      1.75    100.00    100.00     0.000
 4.00      2.00    100.00    100.00     0.000
 4.25      2.25    100.00    100.00     0.000
 4.50      2.50    100.00    100.00     0.000
 4.75      2.75    100.00    100.00     0.000
 5.00      3.00    100.00    100.00     0.000


 ncase =     2    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     38.60     38.60     0.005
 0.50     -1.50     49.73     49.73     0.008
 0.75     -1.25     58.85     58.84     0.010
 1.00     -1.00     66.31     66.30     0.011
 1.25     -0.75     72.42     72.41     0.011
 1.50     -0.50     77.42     77.41     0.011
 1.75     -0.25     81.51     81.51     0.003
 2.00      0.00     84.23     84.86    -0.628
 2.25      0.25     84.79     84.86    -0.073
 2.50      0.50     84.77     84.86    -0.087
 2.75      0.75     84.77     84.86    -0.088
 3.00      1.00     84.77     84.86    -0.087
 3.25      1.25     84.77     84.86    -0.087
 3.50      1.50     84.77     84.86    -0.087
 3.75      1.75     84.77     84.86    -0.087
 4.00      2.00     84.77     84.86    -0.087
 4.25      2.25     84.77     84.86    -0.087
 4.50      2.50     84.77     84.86    -0.087
 4.75      2.75     84.77     84.86    -0.087
 5.00      3.00     84.77     84.86    -0.087


 ncase =     3    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.00     25.00     0.000
 1.00     -1.00     25.00     25.00     0.000
 1.25     -0.75     25.00     25.00     0.000
 1.50     -0.50     25.00     25.00     0.000
 1.75     -0.25     25.00     25.00     0.000
 2.00      0.00     25.00     25.00     0.000
 2.25      0.25     25.00     25.00     0.000
 2.50      0.50     25.00     25.00     0.000
 2.75      0.75     25.00     25.00     0.000
 3.00      1.00     25.00     25.00     0.000
 3.25      1.25     25.00     25.00     0.000
 3.50      1.50     25.00     25.00     0.000
 3.75      1.75     25.00     25.00     0.000
 4.00      2.00     25.00     25.00     0.000
 4.25      2.25     25.00     25.00     0.000
 4.50      2.50     25.00     25.00     0.000
 4.75      2.75     25.00     25.00     0.000
 5.00      3.00     25.00     25.00     0.000


 ncase =     4    h =  1.000e-03

   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff
 0.00     -2.00     25.00     25.00     0.000
 0.25     -1.75     25.00     25.00     0.000
 0.50     -1.50     25.00     25.00     0.000
 0.75     -1.25     25.00     25.00     0.000
 1.00     -1.00     25.00     25.00     0.000
 1.25     -0.75     25.00     25.00     0.000
 1.50     -0.50     25.01     25.00     0.008
 1.75     -0.25     31.29     25.00     6.286
 2.00      0.00     66.42     25.00    41.423
 2.25      0.25     91.65    100.00    -8.351
 2.50      0.50     98.61    100.00    -1.390
 2.75      0.75     99.82    100.00    -0.182
 3.00      1.00     99.98    100.00    -0.021
 3.25      1.25    100.00    100.00    -0.002
 3.50      1.50    100.00    100.00     0.000
 3.75      1.75    100.00    100.00     0.000
 4.00      2.00    100.00    100.00     0.000
 4.25      2.25    100.00    100.00     0.000
 4.50      2.50    100.00    100.00     0.000
 4.75      2.75    100.00    100.00     0.000
 5.00      3.00    100.00    100.00     0.000

\end{verbatim}

We can note the following points about this output:

\vspace{0.5cm}

(1) The van Leer limited gives good accuracy for Cases 1 to 3 which again are
relatively smooth problems. In particular, the van Leer limiter is exact for
Case 3 (the constant solution is differentiated exactly).

\vspace{0.5cm}

(2) The van Leer limiter gives an improvement over centered approximations
(ifd = 1), low order upwinding (ifd = 2) and higher order biased upwinding
(ifd = 3) by reducing further the numerical diffusion and eliminating the
oscillation. Clearly, flux limiters are effective, and their performance can
be improved by going to higher order limiters (e.g., van Leer is second order;
Superbee is third order).

\vspace{0.5cm}

The library differentiator for the van Leer limiter (ifd = 4) is listed below
(this routine is based on a Fortran routine first coded by WES and improved by
Alain vande Wouwer and Philippe Saucez)
\begin{verbatim}

      function [ux]=vanl2(xl,xu,n,u,v)
      dx=(xu-xl)/(n-1);
      delta=1.0e-05;
      if v >= 0.0
        for i=3:n-1
          if(abs(u(i)-u(i-1))<delta)
            phi(i)=0.0;
          else
            r(i)=(u(i+1)-u(i))/(u(i)-u(i-1));
            phi(i)=(r(i)+abs(r(i)))/(1.0+abs(r(i)));
          end
          if(abs(u(i-1)-u(i-2))<delta)
            phi(i-1)=0.0;
          else
            r(i-1)=(u(i)-u(i-1))/(u(i-1)-u(i-2));
            phi(i-1)=(r(i-1)+abs(r(i-1)))/(1.0+abs(r(i-1)));
          end
          flux2=u(i  )+(u(i  )-u(i-1))*phi(i  )/2.0;
          flux1=u(i-1)+(u(i-1)-u(i-2))*phi(i-1)/2.0;
          ux(i)=(flux2-flux1)/dx;
        end
        ux(1)=(-u(1)+u(2))/dx;
        ux(2)=(-u(1)+u(2))/dx;
        ux(n)=(u(n)-u(n-1))/dx;
      end
      if v < 0.0
        for i=2:n-2
          if(abs(u(i)-u(i+1))<delta)
            phi(i)=0.0;
          else
            r(i)=(u(i-1)-u(i))/(u(i)-u(i+1));
            phi(i)=(r(i)+abs(r(i)))/(1.0+abs(r(i)));
          end
          if(abs(u(i+1)-u(i+2))<delta)
            phi(i+1)=0.0;
          else
            r(i+1)=(u(i)-u(i+1))/(u(i+1)-u(i+2));
            phi(i+1)=(r(i+1)+abs(r(i+1)))/(1.0+abs(r(i+1)));
          end
          flux2=u(i  )+(u(i  )-u(i+1))*phi(i  )/2.0;
          flux1=u(i+1)+(u(i+1)-u(i+2))*phi(i+1)/2.0;
          ux(i)=(flux2-flux1)/dx;
        end
      ux(1)=(u(2)-u(1))/dx;
      ux(n-1)=(-u(n-1)+u(n))/dx;
      ux(n  )=(-u(n-1)+u(n))/dx;
      end

\end{verbatim}

The van Leer limiter has two basic sections for positive and negative
velocities and therefore the direction of flow is required as an input. Within
each case (positive or negative velocity),

\vspace{0.5cm}

\begin{itemize}
\item  van Leer computes a ratio at grid point i, r(i), that is a measure of
the change in the slope of the solution

\item  A nonlinear limiting function, phi(i), is then applied to r(i)

\item  From that limited slope, the ``flux'' is computed at two neighboring
grid points

\item  The difference in the two fluxes is then the numerical approximation of
the spatial derivative, ux(i), that is returned from the function to be used
in the PDE

\item  For phi(i) = 0, van Leer reduces to two point upwinding

\item  More generally, van Leer is a form of nonlinear, two point upwinding as
evident when expressed in the following way (from M. Berzins). \vspace
{0.5cm}For the basic hyperbolic equation written in conservative form%

\[
u_{t}+F(u)_{x}=0
\]

the convective term can be approximated as:%

\[
F(u)_{x}=\left[  1+\frac{B(r_{j})}{2}-\frac{B(r_{j-1})}{2r_{j-1}}\right]
\frac{u_{j}-u_{j-1}}{\Delta x}
\]

where%

\[
B(r_{j})=\frac{r_{j}+\left|  r_{j}\right|  }{1+r_{j}},r_{j}=\frac
{u_{j+1}-u_{j}}{u_{j}-u_{j-1}}
\]

Note that when $r_{j}=0,B(r_{j})=0,$ and we have just two point upwinding (ifd
= 2). The nonlinearity of $B(r_{j})$ is clear.
\end{itemize}

\vspace{0.5cm}

To conclude this section, we have considered strongly convective systems that
have as the basic differential group%

\[
u_{t}+F(u)_{x}
\]

or for a constant velocity system%

\[
u_{t}+vu_{x}
\]

If we solve the linear advection equation%

\begin{equation}
u_{t}+vu_{x}=0\tag{1}%
\end{equation}

subject to the initial and boundary conditions%

\[
u(x,0)=f(x),u(0,t)=g(t)
\]

and if the functions $f(x)$ and $g(t)$ introduce a discontinuity, e.g.,%

\begin{equation}
u(x,0)=0,u(0,t)=1\tag{2)(3}%
\end{equation}

we have a form of the \textit{Riemann problem}.

\vspace{0.5cm}

The solution to eqs. (1), (2) and (3) is%

\begin{equation}
u(x,t)=h(x-vt)\tag{4}%
\end{equation}

where $h(t)$ is the Heaviside function%

\[%
\begin{array}
[c]{c}%
h(t)=0,t<0\\
h(t)=1,t>0
\end{array}
\]

In a sense, eq (4) indicates that this Riemann problem is impossible to solve
numerically since at $x-vt=0,$ the derivative $u_{x}$ is infinite. Thus, about
the best we can do is to use an \textit{approximate Riemann solver}$,$ and the
preceding discussion indicates that a flux limiter is a good choice (at least
as good as we have right now). In other words, we will generally not be able
to solve problems exactly that have moving discontinuites, e.g., eq. (4).

\vspace{0.5cm}

On the other hand, we might argue that physically, there are few problems that
have this characteristic; some form of smoothing (e.g., diffusion) is
generally going to occur that will reduce the numerical severity of the
problem. However, problems with steep moving fronts occur in many important
applications such as adsorption and chromatography. Therefore, the continuing
development of effective numerical methods for strongly convective
(hyperbolic) PDEs is an active area of research, particularly the testing of
algorithms, e.g., van Leer, Superbee. For example:

\vspace{0.5cm}

\begin{itemize}
\item  Toro, E. F., {\textit{Riemann Solvers and Numerical Methods for Fluid
Dynamics}}, Springer, Berlin, 1st ed., (1997) has 387 references to the CFD
literature (2nd ed. appeared in 1999)

\item  Hirsch, C., {\textit{Numerical Computation of Internal and External
Flows}}, Vol I: {\textit{Fundamentals of Numerical Discretizations}}, Wiley,
1988; Vol II: {\textit{Computational Methods for Inviscid and Viscous Flows}},
Wiley, 1990
\end{itemize}

\vspace{0.5cm}

Finally, linear approximations can be effective if the problem is sufficiently
smooth, but they are limited by Godunov's theorem which states (loosely):
There are no linear approximations to the Riemann problem, higher than first
order, that are nonoscillatory. The preceding heat exchanger example
demonstrated this conclusion.

\vspace{0.5cm}

The files pertaining to this heat exchanger example are in:

\vspace{0.5cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace}wes1/apci/htex1.ppt\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/htex1.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/dss002.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/dss012.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/dss020.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/vanl2.m\\[0pt]
\end{center}

\vspace{0.5cm}

\textbf{Method of Lines}

\vspace{0.5cm}

The preceding heat exchanger example illustrates a basic approach to the
numerical integration of PDEs generally termed the {\textit{method of lines}}
(MOL). This is essentially a two-step process as the heat exchanger example illustrated:

\vspace{0.5cm}

\begin{itemize}
\item  The spatial (boundary value) derivatives are approximated
algebraically, for example, by finite differences or a flux limiter, on a
spatial grid

\item  The resulting system of ODEs in the initial value variable (typically
time) is then integrated by an initial value ODE integrator (e.g., Euler's
method, the modified Euler method, RKF45, BDF)
\end{itemize}

\vspace{0.5cm}

The advantages of this approach are:

\vspace{0.5cm}

\begin{itemize}
\item  The temporal and spatial integrations are separated in the sense that
they can be treated separately in the coding; this adds a very attractive
degree of flexibility, including the use of library routines for the temporal
and spatial integrations

\item  We can build on the new developments in ODE integrators to do PDEs

\item  All of the major classes of PDEs can be accommodated (elliptic,
hyperbolic, parabolic). Elliptic problems will generally require the addition
of a pseudo initial value derivative that is then integrated to effectively
zero (a steady state or equilibrium condition); this is offer easier than
attacking the elliptic problem directly

\item  Systems of ODEs, DAEs and PDEs (in one, two and three spatial
dimensions plus time), linear and nonlinear, can be naturally accommodated
within the MOL; for example, ODEs or DAEs can serve as boundary conditions for PDEs
\end{itemize}

\vspace{0.5cm}

The MOL is not a single method of solution of ODE/DAE/PDE systems. Rather, it
is a flexible, open-ended approach that is limited only by the ingenuity of
the analyst.

\vspace{0.5cm}

To illustrate the MOL analysis of a PDE system, we return to the problem
considered in the first lecture (28jan00.tex, 28jan00.ps, 28jan00.pdf): A
tubular reactor, as depicted in the accompanying diagram (see file
reactor.ppt), has a significant heat effect so that cooling at the wall is
being considered to remove the heat of reaction.

\vspace{0.5cm}

Consequently, significant radial heat and concentration profiles develop. The
reactor is therefore modeled in terms of radial position, axial position and
time; since there are three independent variables, we will use PDEs.

\vspace{0.5cm}

Our objective in this analysis is to determine what comes out of the reactor,
i.e., the average concentrations and temperature, so that we can determine if
it achieves the required operating performance. To do this, we construct a
mathematical model.

\vspace{0.5cm}

The material balance for an incremental element of length $\Delta z$ is (see
accompanying diagram)%

\[
2\pi r\Delta r\Delta z\frac{\partial c_{a}}{\partial t}=2\pi r\Delta
zq_{m}|_{r}-2\pi(r+\Delta r)\Delta zq_{m}|_{r+\Delta r}
\]
\[
+2\pi r\Delta rvc_{a}|_{z-\Delta z}-2\pi r\Delta rvc_{a}|_{z}
\]
\[
-2\pi r\Delta r\Delta zk_{r}c_{a}^{2}
\]

Division by $2\pi r\Delta r\Delta z$ and minor rearrangement gives%

\[
\frac{\partial c_{a}}{\partial t}=-\frac{(r+\Delta r)q_{m}|_{r+\Delta
r}-rq_{m}|_{r}}{r\Delta r}-v(\frac{c_{a}|_{z}-c_{a}|_{z-\Delta z}}{\Delta
z})-k_{r}c_{a}^{2}
\]

or in the limit $r\rightarrow0,\Delta z\rightarrow0$,%

\[
\frac{\partial c_{a}}{\partial t}=-\frac{1}{r}\frac{\partial(rq_{m})}{\partial
r}-v\frac{\partial c_{a}}{\partial z}-k_{r}c_{a}^{2}
\]

If we now assume Fick's first law for the flux with a constant diffusivity, $D$%

\[
q_{m}=-D\frac{\partial c_{a}}{\partial r}
\]

we obtain%

\[
\frac{\partial c_{a}}{\partial t}=\frac{D}{r}\frac{\partial(r\dfrac{\partial
c_{a}}{\partial r})}{\partial r}-v\frac{\partial c_{a}}{\partial z}-k_{r}%
c_{a}^{2}
\]

or%

\begin{equation}
\frac{\partial c_{a}}{\partial t}=D(\dfrac{\partial^{2}c_{a}}{\partial r^{2}%
}+\frac{1}{r}\dfrac{\partial c_{a}}{\partial r})-v\frac{\partial c_{a}%
}{\partial z}-k_{r}c_{a}^{2}\tag{1}%
\end{equation}

Equation (1) is the required material balance for $c_{a}(r,z,t)$.

\vspace{0.5cm}

The energy balance for the incremental section is%

\[
2\pi r\Delta r\Delta z\rho C_{p}\frac{\partial T}{\partial t}=2\pi r\Delta
zq_{h}|_{r}-2\pi(r+\Delta r)\Delta zq_{h}|_{r+\Delta r}
\]
\[
+2\pi r\Delta rv\rho C_{p}T|_{z-\Delta z}-2\pi r\Delta rv\rho C_{p}T|_{z}
\]
\[
-\Delta H2\pi r\Delta r\Delta zk_{r}c_{a}^{2}
\]

Division by $2\pi r\Delta r\Delta z\rho C_{p}$ gives%

\[
\frac{\partial T}{\partial t}=-\frac{(r+\Delta r)q_{h}|_{r+\Delta r}%
-rq_{h}|_{r}}{r\Delta r\rho C_{p}}-v\frac{(T|_{z}-T|_{z-\Delta z})}{\Delta
z}-\frac{\Delta Hk_{r}}{\rho C_{p}}c_{a}^{2}
\]

or in the limit $\Delta r\rightarrow0,\Delta z\rightarrow0$,%

\[
\frac{\partial T}{\partial t}=-\frac{1}{\rho C_{p}r}\frac{\partial(rq_{h}%
)}{\partial r}-v\frac{\partial T}{\partial z}-\frac{\Delta Hk_{r}}{\rho C_{p}%
}c_{a}^{2}
\]

If we now assume Fourier's first law for the flux with a constant
conductivity, $k$%

\[
q_{h}=-k\frac{\partial T}{\partial r}
\]

we obtain%

\[
\frac{\partial T}{\partial t}=\frac{k}{\rho C_{p}r}\frac{\partial
(r\dfrac{\partial T}{\partial r})}{\partial r}-v\frac{\partial T}{\partial
z}-\frac{\Delta Hk_{r}}{\rho C_{p}}c_{a}^{2}
\]

or%

\begin{equation}
\frac{\partial T}{\partial t}=\frac{k}{\rho C_{p}}(\dfrac{\partial^{2}%
T}{\partial r^{2}}+\frac{1}{r}\dfrac{\partial T}{\partial r})-v\frac{\partial
T}{\partial z}-\frac{\Delta Hk_{r}}{\rho C_{p}}c_{a}^{2}\tag{2}%
\end{equation}

Eq. (2) is the required energy balance for $T(r,z,t).$

\vspace{0.5cm}

The reaction rate constant, $k_{r}$, is given by%

\begin{equation}
k_{r}=k_{0}e^{-E/(RT)}\tag{3}%
\end{equation}

\vspace{0.5cm}

The variables and parameters of eqs. (1), (2) and (3) are summarized in the
following table (in cgs units)

\vspace{0.5cm}
\[%
\begin{array}
[c]{ccc}%
\text{Reactant concentration} & c_{a} & \text{Solution to eq. (1)}\\
\text{Temperature} & T & \text{Solution to eq. (2)}\\
\text{Reaction rate constant} & k_{r} & \text{From eq. (3)}\\
\text{Time} & t & \\
\text{Radial position} & r & \\
\text{Axial position} & z & \\
\text{Entering concentration} & c_{a0} & 0.01\\
\text{Entering temperature} & T_{0} & 305\\
\text{Wall temperature} & T_{w} & 305,355\\
\text{Reactor radius} & r_{0} & 2\\
\text{Reactor length} & z_{l} & 100\\
\text{Linear velocity} & v & 1.0\\
\text{Mass diffusivity} & D & 0.1\\
\text{Thermal diffusivity} & k/(\rho C_{p}) & 0.1\\
\text{Liquid density} & \rho & 1.0\\
\text{Liquid specific heat} & C_{p} & 0.5\\
\text{Heat of reaction} & \Delta H & -10,000\\
\text{Specific rate constant} & k_{0} & 1.5\times10^{9}\\
\text{Activiation energy} & E & 15,000\\
\text{Gas constant} & R & 1.987
\end{array}
\]

\vspace{0.5cm}

Eq. (1) requires one initial condition (IC), two boundary conditions (BCs) in
$r$ and one BC in $z$%

\begin{equation}
c_{a}(r,z,0)=0\tag{4}%
\end{equation}%

\begin{equation}
\frac{\partial c_{a}(0,z,t)}{\partial r}=0,\frac{\partial c_{a}(r_{0}%
,z,t)}{\partial r}=0\tag{5)(6}%
\end{equation}%

\begin{equation}
c_{a}(r,0,t)=c_{a0}\tag{7}%
\end{equation}

Eq. (2) also requires one IC, two BCs in $r$ and one BC in $z$%

\begin{equation}
T(r,z,0)=T_{0}\tag{8}%
\end{equation}%

\begin{equation}
\frac{\partial T(0,z,t)}{\partial r}=0,T(r_{0},z,t)=T_{w}\tag{9)(10}%
\end{equation}%

\begin{equation}
T(r,0,t)=T_{0}\tag{11}%
\end{equation}

The solution to this problem (and in general, to problems in ODEs and PDEs)
are the dependent variables ($c_{a}(r,z,t)$,$T(r,z,t)$) as a function of the
independent variables ($r,z,t$).

\vspace{0.5cm}

Since we are intertested in the exiting conditions, we will compute
$c_{a}(r,z_{l},t),T(r,z_{l},t)$. Also, there may be significant radial
profiles at the exit, so we will also compute the integrals%

\begin{equation}
c_{a,avg}(t)=\int_{0}^{r_{0}}2\pi rc_{a}(r,z_{l},t)dr/(\pi r_{0}^{2})\tag{12}%
\end{equation}%

\begin{equation}
T_{avg}(t)=\int_{0}^{r_{0}}2\pi rT(r,z_{l},t)dr/(\pi r_{0}^{2})\tag{13}%
\end{equation}

Note that this problem includes 2-D PDEs, transcendental equations and
integrals. Because of its complexity (number of equations and their
nonlinearity), an analytical solution is precluded, and we must use numerical methods.

\vspace{0.5cm}

A MOL code for the solution of the preceding problem is listed below
\begin{verbatim}

% Dynamic analysis of a plug flow reactor
%
% The following equations model a plug flow reactor (PFR)
%
%   Material balance
%
%   ca  = -v*ca  + D*(ca   + (1/r)ca ) - k*ca^2                 (1)
%     t        z        rr          r
%
%   Energy balance
%
%   T  = -v*T  + k*/(rho*Cp)*(T   + (1/r)*T )                   (2)
%    t       z                 rr          r
%
%        - (dH*k/(rho*Cp)*ca^2
%
%   kr = k0*exp(-E/(R*T))                                       (3)
%
% The variables and parameters for this model are (in cgs units)
%
%   Reactant concentration          ca        (eq. (1))
%
%   Temperature                     T         (eq. (3))
%
%   Reaction rate constant          kr        (eq. (3))
%
%   Time                            t
%
%   Radial position                 r
%
%   Axial position                  z
%
%   Entering concentration          ca0          0.01
%
%   Entering temperature            T0           305
%
%   Wall temperature                TW           305
%                                                355
%
%   Reactor radius                  r0           2
%
%   Reactor length                  zl           100
%
%   Linear velocity                 v            1
%
%   Mass diffusivity                D            0.1
%
%   Thermal diffusivity     alpha = k/(rho*Cp)   0.1
%
%   Liquid density                  rho          1.0
%
%   Liquid specific heat            Cp           0.5
%
%   Heat of reaction                dH          -10,000
%                                         
%   Specific rate constant          k0           1.5e+09
%                                                2.0e+09
%
%   Activation energy               E            15,000
%
%   Gas constant                    R            1.987
%
% Note that there are two cases, corresponding to the two
% values of the specific rate constant, k0.  This is a
% sensitive parameter since it multiplies a stongly nonlinear
% term (Arrhenius temperature dependency) in eq. (3).  Thus,
% two values are programmed k0 = 1.5e+09, 2.0e+09, to indicate
% the effect of k0.
%
% The solution to this system, ca(r,z,t) (from eq. (1)) and
% T(r,z,t) (from eq. (2)) is computed by a fixed step Euler
% integration in t.  Also, the choice of an Euler step is
% sensitive because of the nonlinearity of eq. (3).
%
% Open output file
  fid1=fopen('reactor.out','w');
%
% Model parameters (defined in the comments above)
  ca0=0.01;
  T0=305.0;
%
% Case 1: Cooled reactor wall
% Tw=305.0;
%
% Case 2: Heated reactor wall
  Tw=355.0;
  r0=2.0;
  zl=100.0;
  v=1.0;
  D=0.1;
  alpha=0.1;
  rho=1.0;
  Cp=0.5;
  dH=-10000.0;
  E=15000.0;
  R=1.987;
%
% Grid in axial direction
  nz=21;
  dz=zl/(nz-1);
  for j=1:nz
    z(j)=(j-1)*dz;
  end
%
% Grid in radial direction
  nr=5;
  dr=r0/(nr-1);
  for i=1:nr
    r(i)=(i-1)*dr;
  end
  drs=dr^2;
%
% Case 1: Low reaction rate constant
% rk0=1.5e+09;
%
% Case 2: High reaction rate constant
  rk0=2.0e+09;
%
% Series of solutions for different Euler steps
  for ncase=1:1
%
% Initial, final times
  t=0.0; tf=300.0;
%
% Integration step, number of Euler steps for each output
  if ncase==1 h=0.1; nout=100; end
%
% Initial conditions
  for j=1:nz
  for i=1:nr
    ca(i,j)=0.0;
     T(i,j)=T0;
  end
  end
%
% Write h
  fprintf(fid1,'\n\n ncase = %5d    h = %10.3e\n\n',ncase,h);
%
% Write heading
  fprintf(fid1,'     t   ca(r,zl,t)  T(r,zl,t)\n');
%
% Integrate until t = tf
  while t < tf*1.0001
%
% Monitor solution by displaying t
  t
%
% Write selected output (a quadrature, e.g., Simpson's
% rule, could be added here for the average output ca, T)
  for i=1:nr
    fprintf(fid1,'%6.1f%12.6f%10.2f\n',t,ca(i,nz),T(i,nz));
  end
  fprintf(fid1,'\n');
%
% Take nout Euler steps
  for iout=1:nout
%
%   Temporal derivatives
%
%     Cover the nr x nz grid
%
%       Entering conditions (z = 0)
        for i=1:nr
          ca(i,1)=ca0;
          cat(i,1)=0.0;
          T(i,1)=T0;
          Tt(i,1)=0.0;
        end
%
%       Rest of reactor
        for j=2:nz
%
%         Centerline (r = 0)
          rk=rk0*exp(-E/(R*T(1,j)));
          cat(1,j)=-v*(ca(1,j)-ca(1,j-1))/dz...
                  +4.0*D*(ca(2,j)-ca(1,j))/drs...
                  -rk*ca(1,j)^2;
          Tt(1,j)=-v*(T(1,j)-T(1,j-1))/dz...
                  +4.0*alpha*(T(2,j)-T(1,j))/drs...
                  -dH*rk/(rho*Cp)*ca(1,j)^2;
%
%         Wall (r = r0)
          rk=rk0*exp(-E/(R*Tw));
          cat(nr,j)=-v*(ca(nr,j)-ca(nr,j-1))/dz...
                   +2.0*D*(ca(nr-1,j)-ca(nr,j))/drs...
                   -rk*ca(nr,j)^2;
           T(nr,j)=Tw;
          Tt(nr,j)=0.0;
%
%         Interior, r ~= 0 and r ~= r0
          for i=2:nr-1
            rk=rk0*exp(-E/(R*T(i,j)));
            cat(i,j)=-v*(ca(i,j)-ca(i,j-1))/dz...
                   +D*(ca(i+1,j)-2.0*ca(i,j)+ca(i-1,j))/drs...
                   +D*(1.0/r(i)*(ca(i+1)-ca(i-1)))/(2.0*dr)...
                   -rk*ca(i,j)^2;
            Tt(i,j)=-v*(T(i,j)-T(i,j-1))/dz...
                   +alpha*(T(i+1,j)-2.0*T(i,j)+T(i-1,j))/drs...
                   +alpha*(1.0/r(i)*(T(i+1,j)-T(i-1,j)))/(2.0*dr)...
                   -dH*rk/(rho*Cp)*ca(i,j)^2;
%
%         Next r
          end
%
%       Next z
        end
%
%   All temporal derivatives are computed
%
%   Take Euler step
    for j=1:nz
    for i=1:nr
      ca(i,j)=ca(i,j)+cat(i,j)*h;
       T(i,j)= T(i,j)+ Tt(i,j)*h;
    end
    end
    t=t+h;
%
% Next Euler step
  end
%
% Next output
  end
%
% Next case
  end

\end{verbatim}

\vspace{0.5cm}

We can note the following points about this program:

\vspace{0.5cm}

(1) The model parameters are first defined, including the variation in the
wall temperature, Tw (initially, consideration was given to cooling the wall
by keeping it at the initial temperature Tw = 305, but execution of the
program indicated that very little reaction occurred; thus, the wall in the
second case is actually heated to Tw = 355 to produce a significant reaction
as reflected in the consumption of the reactant, i.e., reduction of the
entering concentration $c_{a} (r,0,t) = 0.01$).
\begin{verbatim}

% Model parameters (defined in the comments above)
  ca0=0.01;
  T0=305.0;
%
% Case 1: Cooled reactor wall
% Tw=305.0;
%
% Case 2: Heated reactor wall
  Tw=355.0;
  r0=2.0;
  zl=100.0;
  v=1.0;
  D=0.1;
  alpha=0.1;
  rho=1.0;
  Cp=0.5;
  dH=-10000.0;
  E=15000.0;
  R=1.987;

\end{verbatim}

(2) The spatial grids in z and r are then defined numerically over the
intervals $0 \leq z \leq z_{l}$ and $0 \leq r \leq r_{0}$ (with nz = 21 and nr
= 5 points, respectively)
\begin{verbatim}

% Grid in axial direction
  nz=21;
  dz=zl/(nz-1);
  for j=1:nz
    z(j)=(j-1)*dz;
  end
%
% Grid in radial direction
  nr=5;
  dr=r0/(nr-1);
  for i=1:nr
    r(i)=(i-1)*dr;
  end
  drs=dr^2;

\end{verbatim}

(3) Variation of the specific rate constant, $k$ is also included (this can be
considered as a type of parameter sensitivity analysis, perhaps to compare the
model output with experimental data)
\begin{verbatim}

% Case 1: Low reaction rate constant
% rk0=1.5e+09;
%
% Case 2: High reaction rate constant
  rk0=2.0e+09;

\end{verbatim}

(4) The code provides for multiple runs, but only one run is coded. Within
each run, the parameters controlling the integration in $t$ are set, including
the initial and final values of $t$, the euler integration step, h, and the
number of Euler steps between outputs, nout:
\begin{verbatim}

% Series of solutions for different Euler steps
  for ncase=1:1
%
% Initial, final times
  t=0.0; tf=300.0;
%
% Integration step, number of Euler steps for each output
  if ncase==1 h=0.1; nout=100; end

\end{verbatim}

(5) The initial conditions for the two dependent variables, $c_{a}$ (eq. (5))
and $T$ (eq. (8), are set
\begin{verbatim}

% Initial conditions
  for j=1:nz
  for i=1:nr
    ca(i,j)=0.0;
     T(i,j)=T0;
  end
  end

\end{verbatim}

Note that there are nz x nr = 21 x 5 = 105 ODEs in the MOL approximation of
the two PDEs, eqs. (1) and (2)

\vspace{0.5cm}

(6) Headings are printed for the output
\begin{verbatim}

% Write h
  fprintf(fid1,'\n\n ncase = %5d    h = %10.3e\n\n',ncase,h);
%
% Write heading
  fprintf(fid1,'     t   ca(r,zl,t)  T(r,zl,t)\n');

\end{verbatim}

(7) The integration in $t$ starts with a test if t has reached the final
value, tf
\begin{verbatim}

% Integrate until t = tf
  while t < tf*1.0001
%
% Monitor solution by displaying t
  t

\end{verbatim}

(8) The solution is printed (including the initial conditions in the first
pass through this code)
\begin{verbatim}

% Write selected output (a quadrature, e.g., Simpson's
% rule, could be added here for the average output ca, T)
  for i=1:nr
    fprintf(fid1,'%6.1f%12.6f%10.2f\n',t,ca(i,nz),T(i,nz));
  end
  fprintf(fid1,'\n');

\end{verbatim}

If the average exiting concentration and temperature, according to eqs. (12)
and (13), are to be included in the output, a quadrature calculation, e.g.,
Simpson's rule, could be included at this point (the radial grid has an odd
number of points, nr = 5, as required by Simpson's rule)

\vspace{0.5cm}

(9) nout Euler steps are taken (between outputs), starting with the definition
of the boundary conditions
\begin{verbatim}

% Take nout Euler steps
  for iout=1:nout
%
%   Temporal derivatives
%
%     Cover the nr x nz grid
%
%       Entering conditions (z = 0)
        for i=1:nr
          ca(i,1)=ca0;
          cat(i,1)=0.0;
          T(i,1)=T0;
          Tt(i,1)=0.0;
        end

\end{verbatim}

Note that the entering conditions, $c_{a}(r,0,t)$ and $T(r,0,t)$ are set
according to BCs (7) and (11). Since these entering conditions are constant,
their time derivatives are set to zero. All four of these conditions are
specified with the second index set to one which is the first grid point in
$z$ (the entrance to the reactor).

\vspace{0.5cm}

(10) The remaining points in $z$, $2 \leq j \leq nz (=21)$ are then covered.
The calculations for $r=0$ (with the first index set to one) are done first

\vspace{0.5cm}
\begin{verbatim}

%
%       Rest of reactor
        for j=2:nz
%
%         Centerline (r = 0)
          rk=rk0*exp(-E/(R*T(1,j)));
          cat(1,j)=-v*(ca(1,j)-ca(1,j-1))/dz...
                  +4.0*D*(ca(2,j)-ca(1,j))/drs...
                  -rk*ca(1,j)^2;
          Tt(1,j)=-v*(T(1,j)-T(1,j-1))/dz...
                  +4.0*alpha*(T(2,j)-T(1,j))/drs...
                  -dH*rk/(rho*Cp)*ca(1,j)^2;

\end{verbatim}

The Arrhenius temperature dependency is first calculated (note the use of the
centerline temperature, $T(0,z,t)$)
\begin{verbatim}

          rk=rk0*exp(-E/(R*T(1,j)));

\end{verbatim}

The temporal derivative from eq. (1)
\[
\frac{\partial c_{a}}{\partial t}
\]

is then calculated
\begin{verbatim}

          cat(1,j)=-v*(ca(1,j)-ca(1,j-1))/dz...
                  +4.0*D*(ca(2,j)-ca(1,j))/drs...
                  -rk*ca(1,j)^2;

\end{verbatim}

We can note the following points in this coding of the temporal derivative:

\vspace{0.5cm}

(10.1) The convective derivative in eq. (1)%

\[
-v\frac{\partial c_{a}}{\partial z}
\]

is approximated by a two point upwind FD
\begin{verbatim}

                   -v*(ca(1,j)-ca(1,j-1))/dz

\end{verbatim}

The calculation of $\dfrac{\partial c_{a}}{\partial z}$ could also have been
done by a call to function dss012 as illustrated by the preceding heat
exchanger example. Note also that since the index in $z$ starts at j = 2, the
value of ca(1,j-1) is defined (by BC (7), discussed in (9) above).

\vspace{0.5cm}

(10.2) The diffusive derivatives in eq. (1)%

\[
D(\dfrac{\partial^{2}c_{a}}{\partial r^{2}}+\frac{1}{r}\dfrac{\partial c_{a}%
}{\partial r})
\]

are approximated by three point, centered FDs, including BC (5)
\begin{verbatim}

                  +4.0*D*(ca(2,j)-ca(1,j))/drs

\end{verbatim}

There is a complication at $r=0$ that was taken into account in the preceding
coding. The second term at $r=0$ is indetermine as a consquence of BC (5)%

\[
\lim_{r\rightarrow0}\frac{1}{r}\frac{\partial c_{a}}{\partial r}=\frac{0}{0}
\]

If we apply l'Hospital's rule to the indeterminant form by differentiating the
numerator and denominator with respect to $r$%

\[
\lim_{r\rightarrow0}\frac{1}{r}\frac{\partial c_{a}}{\partial r}%
=\lim_{r\rightarrow0}\frac{1}{1}\frac{\partial^{2}c_{a}}{\partial r^{2}}%
=\frac{\partial^{2}c_{a}}{\partial r^{2}}
\]

Thus, at $r=0,$ the diffusion group becomes%

\[
D(\dfrac{\partial^{2}c_{a}}{\partial r^{2}}+\frac{1}{r}\dfrac{\partial c_{a}%
}{\partial r})=2D\frac{\partial^{2}c_{A}}{\partial r^{2}}
\]

with the FD approximation%

\[
2D\frac{\partial^{2}c_{a}}{\partial r^{2}}\cong2D\frac{c_{a}(2)-2c_{a}%
(1)+c_{a}(0)}{\Delta r^{2}}
\]

We then eliminate the fictitious value $c_{a}(0)$ using BC (5)%

\[
\frac{\partial c_{a}(0,z,t)}{\partial r}\cong\frac{c_{a}(2)-c_{a}(0)}{2\Delta
r}=0
\]

or%

\[
c_{a}(0)=c_{a}(2)
\]

The final result for the diffusion group at $r=0$ is then%

\[
D(\dfrac{\partial^{2}c_{a}}{\partial r^{2}}+\frac{1}{r}\dfrac{\partial c_{a}%
}{\partial r})\cong2D\frac{c_{a}(2)-2c_{a}(1)+c_{a}(0)}{\Delta r^{2}}%
\cong4D\frac{c_{a}(2)-c_{a}(1)}{\Delta r^{2}}
\]

which is programmed as indicated above.

\vspace{0.5cm}

(10.3) The reaction term%

\[
-k_{r}c_{a}^{2}
\]

is programmed as
\begin{verbatim}

                  -rk*ca(1,j)^2;

\end{verbatim}

Note the ease with which the nonlinear, second order reaction term can be
programmed (illustrating the power of numerical methods, i.e., nonlinearities
are easily included in the analysis).

\vspace{0.5cm}

(10.4) The temporal derivative of temperature from eq. (2),%

\[
\frac{\partial T}{\partial t}
\]

is programmed in the same way
\begin{verbatim}

          Tt(1,j)=-v*(T(1,j)-T(1,j-1))/dz...
                  +4.0*alpha*(T(2,j)-T(1,j))/drs...
                  -dH*rk/(rho*Cp)*ca(1,j)^2;

\end{verbatim}

Again, the ease of programming the nonlinear Arrhenius temperature dependency
in the reaction term (discussed above)
\begin{verbatim}

          rk=rk0*exp(-E/(R*T(1,j)));

\end{verbatim}

is apparent.

\vspace{0.5cm}

(11) The temporal derivatives at the wall are then programmed
\begin{verbatim}

%         Wall (r = r0)
          rk=rk0*exp(-E/(R*Tw));
          cat(nr,j)=-v*(ca(nr,j)-ca(nr,j-1))/dz...
                   +2.0*D*(ca(nr-1,j)-ca(nr,j))/drs...
                   -rk*ca(nr,j)^2;
           T(nr,j)=Tw;
          Tt(nr,j)=0.0;

\end{verbatim}

The Arrhenius temperature is first programmed using the wall temperature, Tw
\begin{verbatim}

          rk=rk0*exp(-E/(R*Tw));

\end{verbatim}

The temporal derivative from eq. (1)%

\[
\frac{\partial c_{a}}{\partial t}
\]

is then calculated
\begin{verbatim}

         cat(nr,j)=-v*(ca(nr,j)-ca(nr,j-1))/dz...
                   +2.0*D*(ca(nr-1,j)-ca(nr,j))/drs...
                   -rk*ca(nr,j)^2;

\end{verbatim}

We can note the following points in this coding of the temporal derivative:

\vspace{0.5cm}

(11.1) The convective derivative in eq. (1)%

\[
-v\frac{\partial c_{a}}{\partial z}
\]

is again approximated by a two point upwind FD
\begin{verbatim}

                   -v*(ca(nr,j)-ca(nr,j-1))/dz
                 
\end{verbatim}

Note the use of the subscript nr since we are now considering eq. (1) at the
wall ($r=r_{0}$).

\vspace{0.5cm}

(11.2) The diffusive derivatives in eq. (1)%

\[
D(\dfrac{\partial^{2}c_{a}}{\partial r^{2}}+\frac{1}{r}\dfrac{\partial c_{a}%
}{\partial r})
\]

are approximated by three point, centered FDs, including BC (6)
\begin{verbatim}

                   +2.0*D*(ca(nr-1,j)-ca(nr,j))/drs          

\end{verbatim}

To review the details, since at the wall (from BC (6)),%

\[
\dfrac{\partial c_{a}(r_{0},z,t)}{\partial r}=0
\]

the second term in the diffusion group is zero and only the first term has to
be considered%

\[
D\dfrac{\partial^{2}c_{a}}{\partial r^{2}}\cong D\frac{c_{a}(nr+1)-2c_{a}%
(nr)+c_{a}(nr-1)}{\Delta r^{2}}
\]

But $c_{a}(nr+1)$ is a fictitious value that can be eliminated by a FD
approximation of BC (6)%

\[
\dfrac{\partial c_{a}(r_{0},z,t)}{\partial r}\cong\frac{c_{a}(nr+1)-c_{a}%
(nr-1)}{2\Delta r}=0
\]

or%

\[
c_{a}(nr+1)=c_{a}(nr-1)
\]

and%

\[
D\dfrac{\partial^{2}c_{a}}{\partial r^{2}}\cong D\frac{c_{a}(nr+1)-2c_{a}%
(nr)+c_{a}(nr-1)}{\Delta r^{2}}=2D\frac{c_{a}(nr-1)-c_{a}(nr)}{\Delta r^{2}}
\]

which is programmed as indicated above.

\vspace{0.5cm}

(11.3) The reaction term%

\[
-k_{r}c_{a}^{2}
\]

is programmed as
\begin{verbatim}

                  -rk*ca(nr,j)^2;

\end{verbatim}

(again, using the subscript nr).

(11.4) The energy balance is easily programmed since the wall temperature is constant
\begin{verbatim}

           T(nr,j)=Tw;
          Tt(nr,j)=0.0;

\end{verbatim}

Note that the temporal derivative is zeroed since the wall temperature is constant.

\vspace{0.5cm}

(12) The temporal derivatives of eqs. (1) and (2) at the interrior points (
with $r\neq0$ and $r\neq r_{0})$ are then programmed
\begin{verbatim}

%         Interior, r ~= 0 and r ~= r0
          for i=2:nr-1
            rk=rk0*exp(-E/(R*T(i,j)));
            cat(i,j)=-v*(ca(i,j)-ca(i,j-1))/dz...
                   +D*(ca(i+1,j)-2.0*ca(i,j)+ca(i-1,j))/drs...
                   +D*(1.0/r(i)*(ca(i+1)-ca(i-1)))/(2.0*dr)...
                   -rk*ca(i,j)^2;
            Tt(i,j)=-v*(T(i,j)-T(i,j-1))/dz...
                   +alpha*(T(i+1,j)-2.0*T(i,j)+T(i-1,j))/drs...
                   +alpha*(1.0/r(i)*(T(i+1,j)-T(i-1,j)))/(2.0*dr)...
                   -dH*rk/(rho*Cp)*ca(i,j)^2;
%         Next r
          end
%
%       Next z
        end

\end{verbatim}

This programming appears to be relatively complicated, but actually it is just
the systematic application of two point upwind FD approximations for the
convective derivative in $z$, and three point centered approximations for the
diffusive derivatives in $r$ (no special cases have to be considered as with
$r=0$ and $r=r_{0}$). Note that this coding concludes the loop in $z$ (index
j) initiated at the beginning of the programming of the temporal derivatives,
so that two end statements are required.

\vspace{0.5cm}

(13) Finally, all nz x nr dependent variables are advanced in time through the
application of Euler's method
\begin{verbatim}

%   All temporal derivatives are computed
%
%   Take Euler step
    for j=1:nz
    for i=1:nr
      ca(i,j)=ca(i,j)+cat(i,j)*h;
       T(i,j)= T(i,j)+ Tt(i,j)*h;
    end
    end
    t=t+h;

\end{verbatim}

(14) The time stepping is concluded, followed by the looping for output (after
nout Euler steps) and finally the loop for different cases is concluded
\begin{verbatim}

% Next Euler step
  end
%
% Next output
  end
%
% Next case
  end

\end{verbatim}

In summary, although the preceding coding is relatively detailed, it follows
the usual MOL two-step process:

\vspace{0.5cm}

\begin{itemize}
\vspace{0.5cm}

\item  Evaluation of the initial value derivatives (RHSs of the PDEs), which
includes the calculation of the spatial derivatives and application of the
BCs, to produce a set of approximating ODEs

\vspace{0.5cm}

\item  Integration of the ODEs to produce a new set of dependent variable
values that are then used to evaluate now initial value derivatives
\end{itemize}

\vspace{0.5cm}

The output from this program is:
\begin{verbatim}

 ncase =     1    h =  1.000e-01

     t   ca(r,zl,t)  T(r,zl,t)
   0.0    0.000000    305.00
   0.0    0.000000    305.00
   0.0    0.000000    305.00
   0.0    0.000000    305.00
   0.0    0.000000    305.00

  10.0    0.000000    333.05
  10.0    0.000000    334.51
  10.0    0.000000    338.58
  10.0    0.000000    344.22
  10.0    0.000000    350.00
           .             .
           .             .
  50.0    0.000030    349.94
  50.0    0.000030    349.94
  50.0    0.000030    349.96
  50.0    0.000030    349.98
  50.0    0.000030    350.00
           .             .
           .             .
 100.0    0.003540    352.84
 100.0    0.003541    352.65
 100.0    0.003544    352.08
 100.0    0.003547    351.18
 100.0    0.003549    350.00
           .             .
           .             .
 150.0    0.004622    355.44
 150.0    0.004624    355.06
 150.0    0.004632    353.96
 150.0    0.004643    352.23
 150.0    0.004649    350.00
           .             .
           .             .
 200.0    0.004624    355.45
 200.0    0.004626    355.07
 200.0    0.004634    353.96
 200.0    0.004645    352.23
 200.0    0.004651    350.00
           .             .
           .             .
 250.0    0.004624    355.45
 250.0    0.004626    355.07
 250.0    0.004634    353.96
 250.0    0.004645    352.23
 250.0    0.004651    350.00
           .             .
           .             .
 300.0    0.004624    355.45
 300.0    0.004626    355.07
 300.0    0.004634    353.96
 300.0    0.004645    352.23
 300.0    0.004651    350.00

\end{verbatim}

\vspace{0.5cm}

This output would be a good candidate for plotting, or better, visualization.
For example, centerline profiles could be plotted, ca(0,z,t), T(0,z,t).

\vspace{0.5cm}

Also, some experimentation with the temporal and spatial integration should be
done, for example:

\begin{itemize}
\item  Use a higher (than Euler) order ODE integrator such as RKF45. Also, the
Euler coding could be replaced by a library integrator such as Matlab ode23,
ode45. This basically amounts to p refinement in time

\item  The ODE integration step $h$ could also be varied; this amounts to h
refinement in time

\item  Use a higher (than two point upwind) spatial differentiator. FDs would
probably be ok for the same reason that they were satisfactory in the heat
exchanger, Case = 1, 2; i.e., the solution is probably sufficiently smooth for
FDs to give give accuracy, free of numerical diffusion and oscillation

\begin{itemize}
\item  This amounts to p refinement in space

\item  The spatial differentiation can be programmed with library routines
rather than explicitly

\item  Centered approximations should be used for the diffusion spatial
derivatives in r (there is no preferred direction)

\item  Noncentered approximations should be used for the convection
derivatives in z (some form of upwinding is required; centered approximations
should not be used)
\end{itemize}

\item  The spatial integration intervals in r and z could be varied (the
number of grid points in r and z); the amounts to h refinement in space
\end{itemize}

\vspace{0.5cm}

The files pertaining to this tubular reactor example are in:

\vspace{0.5cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace}wes1/apci/reactor.ppt\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/reactor.m\\[0pt]
\end{center}

\vspace{0.5cm}

Finally,

\vspace{0.5cm}

``The method of lines has become the most widely used solution technique for
large-scale time-dependent partial differential equations", G. Fairweather and
I. Gladwell, Preface to the special issue on the method of lines, Applied
Numerical Mathematics, v. 20, 1996, pp 1-2

\vspace{0.5cm}

In the preceding examples, we have used FD approximations that were either
preprogrammed, e.g., in dss002.m, vanl2.m, or were merely stated, then used,
e.g., in the reactor application. We now consider a systematic procedure for
generating FD approximations for virtually any situation, that is, FDs of any
order for derivatives of any order, including BCs.

\vspace{0.5cm}

The files pertaining to this heat exchanger example are in:

\vspace{0.5cm}

\begin{center}
http://www.lehigh.edu/\~{\thinspace}wes1/apci/htex1.ppt\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/htex1.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/dss002.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/dss012.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/dss020.m\\[0pt]%
http://www.lehigh.edu/\~{\thinspace}wes1/apci/vanl2.m\\[0pt]
\end{center}

\vspace{0.5cm}\textbf{Finite Differences}

\vspace{0.5cm}

If we applying the MOL to the one-dimensional Fourier's second law in
cylindrical coordinates,%

\begin{equation}
\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial r^{2}}+\frac{1}%
{r}\frac{\partial u}{\partial r}\tag{1}%
\end{equation}

or in subscript notation%

\[
u_{t}=u_{rr}+(1/r)u_{r}
\]

we need a FD approximation for the first and second derivatiuves, $u_{r}$ and
$u_{rr}.$

\vspace{0.5cm}

Considering the first derivative, we define a spatial grid in $r$ with $N$
grid points and index $i=1,2,\cdots,N$, the left grid point corresponds to
$i=1$. We can use two Taylor series centered at grid point $i=1$%

\[
u_{2}=u_{1}+\frac{du_{1}}{dr}\Delta r+\frac{d^{2}u_{1}}{dr^{2}}\frac{\Delta
r^{2}}{2!}+\frac{d^{3}u_{1}}{dr^{3}}\frac{\Delta r^{3}}{3!}+\cdots
\]%

\[
u_{0}=u_{1}+\frac{du_{1}}{dr}(-\Delta r)+\frac{d^{2}u_{1}}{dr^{2}}%
\frac{(-\Delta r)^{2}}{2!}+\frac{d^{3}u_{1}}{dr^{3}}\frac{(-\Delta r)^{3}}%
{3!}+\cdots
\]

Subtraction of the second series from the first gives
\[
u_{2}-u_{0}=2\frac{du_{1}}{dr}\frac{\Delta r}{1!}+2\frac{d^{3}u_{1}}{dr^{3}%
}\frac{\Delta r^{3}}{3!}\cdots
\]

or, since we are interested in an approximation for the first derivative,%

\[
\frac{du_{1}}{dr}\cong\frac{u_{2}-u_{0}}{2\Delta r}-\frac{d^{3}u_{1}}{dr^{3}%
}\frac{\Delta r^{2}}{3!}
\]

and we therefore have the well known centered, second order FD approximation
for a first derivative%

\begin{equation}
\frac{du_{1}}{dr}\cong\frac{u_{2}-u_{0}}{2\Delta r}+O(\Delta r^{2})\tag{2}%
\end{equation}

Note, however, that eq. (2) uses $u_{0}$ which is outside the grid, i.e., it
is $u$ at a fictitious point.

\vspace{0.5cm}

We, however, can use another approach that eliminates the use of fictitious
points, starting with Taylor series expansions based on only the boundary and
interior points. Thus,%

\[
u_{2}=u_{1}+\frac{du_{1}}{dr}\Delta r+\frac{d^{2}u_{1}}{dr^{2}}\frac{\Delta
r^{2}}{2!}+\frac{d^{3}u_{1}}{dr^{3}}\frac{\Delta r^{3}}{3!}+\cdots
\]%

\[
u_{3}=u_{1}+\frac{du_{1}}{dr}2\Delta r+\frac{d^{2}u_{1}}{dr^{2}}\frac{(2\Delta
r)^{2}}{2!}+\frac{d^{3}u_{1}}{dr^{3}}\frac{(2\Delta r)^{3}}{3!}+\cdots
\]

The derivative of interest is $\dfrac{du_{1}}{dr}$, and we wish to drop as
many terms past this derivative as possible to achieve maximum accuracy.

\vspace{0.5cm}

If we multiply the first Taylor series by $4$, and subtract the second Taylor
series from it, we can drop out the second derivative terms (involving
$\dfrac{d^{2}u_{1}}{dr^{2}}$)%

\[
4\left(  u_{2}=u_{1}+\frac{du_{1}}{dr}\Delta r+\frac{d^{2}u_{1}}{dr^{2}}%
\frac{\Delta r^{2}}{2!}+\frac{d^{3}u_{1}}{dr^{3}}\frac{\Delta r^{3}}%
{3!}+\cdots\right)
\]
\[
u_{3}=u_{1}+\frac{du_{1}}{dr}2\Delta r+\frac{d^{2}u_{1}}{dr^{2}}\frac{(2\Delta
r)^{2}}{2!}+\frac{d^{3}u_{1}}{dr^{3}}\frac{(2\Delta r)^{3}}{3!}+\cdots
\]

or%

\[
4u_{2}-u_{3}=3u_{1}+2\frac{du_{1}}{dr}\Delta r+0\frac{d^{2}u_{1}}{dr^{2}}%
\frac{\Delta r^{2}}{2!}+(-4)\frac{d^{3}u_{1}}{dr^{3}}\frac{\Delta r^{3}}%
{3!}+\cdots
\]

or, truncating after the third derivative terms, then solving for the
derivative of interest, $\dfrac{du_{1}}{dr},$%

\[
\dfrac{du_{1}}{dr}=\frac{-3u_{1}+4u_{2}-u_{3}}{2\Delta r}+\frac{4}{2}%
\frac{d^{3}u_{1}}{dr^{3}}\frac{\Delta r^{2}}{3!}
\]

Note that this second order approximation for the first derivative requires
only boundary ($u_{1}$) and interior ($u_{2},u_{3}$) values.

\vspace{0.5cm}

A similar analysis at the right boundary (grid point $i=n$) gives
\begin{equation}
\dfrac{du_{N}}{dr}=\frac{3u_{N}-4u_{N-1}+u_{N-2}}{2\Delta r}+O(\Delta
r^{2})\tag{3}%
\end{equation}

At the interior points, we can use the centered approximation (this follows
from eq. (2))%

\begin{equation}
\dfrac{du_{i}}{dr}=\frac{u_{i+1}-u_{i-1}}{2\Delta r}+O(\Delta r^{2})\tag{4}%
\end{equation}

Then for the entire domain, we can form a differentiation matrix (or
computational stencil)

\vspace{0.5cm}%

\[
\frac{d\mathbf{u}}{dr}=\frac{1}{2\Delta r}\left[
\begin{array}
[c]{ccc}%
-3 & 4 & -1\\
-1 & 0 & 1\\
1 & -4 & 3
\end{array}
\right]  \mathbf{u}+O(\Delta r^{2})
\]

\vspace{0.5cm}

where row 1 is used for grid point $i=1,$ row 2 for $i=2$ to $i=N-1$, and row
3 for $i=N$. Note that the matrix is antisymmetric with respect to the center
element (antisymmetric around $0$).

\vspace{0.5cm}

To make explicit the points at which each row of the differentiation is used,
we can write it as%

\begin{equation}
\frac{d\mathbf{u}}{dr}=\frac{1}{2\Delta r}\left[
\begin{array}
[c]{cccc}%
-3 & 4 & -1 & i=1\\
-1 & 0 & 1 & i=2,\cdots,N-1\\
1 & -4 & 3 & i=N
\end{array}
\right]  \mathbf{u}+O(\Delta r^{2})\tag{5}%
\end{equation}

The differentiation matrix of eq. (5) is programmed in function dss002.m.

\vspace{0.5cm}

Again, note that FD approximations are linear, i.e., the derivative at each
grid point is calculated as a weighted sum of the neighboring dependent
variable values, as illustrated by eqs. (3) and (4).

\vspace{0.5cm}

To obtain a FD approximation for the second derivative in eq. (1), we can use
the same two Taylor series as before (but now written at a grid point $i$ to
give a general result)%

\[
u_{i+1}=u_{i}+\frac{du_{i}}{dr}\Delta r+\frac{d^{2}u_{i}}{dr^{2}}\frac{\Delta
r^{2}}{2!}+\frac{d^{3}u_{i}}{dr^{3}}\frac{\Delta r^{3}}{3!}+\cdots
\]%

\[
u_{i-1}=u_{i}+\frac{du_{i}}{dr}(-\Delta r)+\frac{d^{2}u_{i}}{dr^{2}}%
\frac{(-\Delta r)^{2}}{2!}+\frac{d^{3}u_{i}}{dr^{3}}\frac{(-\Delta r)^{3}}%
{3!}+\cdots
\]

Adding these two series (rather than subtracting them as before) gives%

\[
u_{i+1}+u_{i-1}=2u_{i}+2\frac{d^{2}u_{i}}{dr^{2}}\frac{\Delta r^{2}}%
{2!}+2\frac{d^{4}u_{i}}{dr^{4}}\frac{\Delta r^{4}}{4!}
\]

Solving for the second derivative, we have%

\[
\frac{d^{2}u_{i}}{dr^{2}}\cong\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta r^{2}%
}-2\frac{d^{4}u_{i}}{dr^{4}}\frac{\Delta r^{2}}{4!}
\]

or%

\begin{equation}
\frac{d^{2}u_{i}}{dr^{2}}\cong\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta r^{2}%
}+O(h^{2})\tag{6}%
\end{equation}

which is the well known centered, second order, FD approximation for the
second derivative.

\vspace{0.5cm}

Note again that eq. (6) references $u_{0}$ when $i=1$, so we would have to
somehow provide this fictitious value, generally through a BC at $i=1.$ An
alternate approach is to use a noncentered difference for the second
derivative at $i=1$ that requires only the boundary and interior points. The
derivation is a bit more involved than for first derivatives, but it has been
done and programmed in library spatial differentiators.

\vspace{0.5cm}

For example, the code is given below for a spatial differentiation routine (dss042.m):

\vspace{0.5cm}
\begin{verbatim}

      function [uxx]=dss042(xl,xu,n,u,ux,nl,nu)
%...
%...  FUNCTION DSS042 COMPUTES A SECOND-ORDER APPROXIMATION OF A
%...  SECOND-ORDER DERIVATIVE, WITH OR WITHOUT THE NORMAL DERIVATIVE
%...  AT THE BOUNDARY.
%...
%...  ARGUMENT LIST
%...
%...     XL      LEFT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
%...
%...     XU      RIGHT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
%...
%...     N       NUMBER OF SPATIAL GRID POINTS, INCLUDING THE END
%...             POINTS (INPUT)
%...
%...     U       ONE-DIMENSIONAL ARRAY OF THE DEPENDENT VARIABLE TO BE
%...             DIFFERENTIATED (INPUT)
%...
%...     UX      ONE-DIMENSIONAL ARRAY OF THE FIRST DERIVATIVE OF U.
%...             THE END VALUES OF UX, UX(1) AND UX(N), ARE USED IN
%...             NEUMANN BOUNDARY CONDITIONS AT X = XL AND X = XU,
%...             DEPENDING ON THE ARGUMENTS NL AND NU (SEE THE DE-
%...             SCRIPTION OF NL AND NU BELOW)
%...
%...     UXX     ONE-DIMENSIONAL ARRAY OF THE SECOND DERIVATIVE OF U
%...             (OUTPUT)
%...
%...     NL      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
%...             X = XL (INPUT).  THE ALLOWABLE VALUES ARE
%...
%...                1 - DIRICHLET BOUNDARY CONDITION AT X = XL
%...                    (UX(1) IS NOT USED)
%...
%...                2 - NEUMANN BOUNDARY CONDITION AT X = XL
%...                    (UX(1) IS USED)
%...
%...     NU      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
%...             X = XU (INPUT).  THE ALLOWABLE VALUES ARE
%...
%...                1 - DIRICHLET BOUNDARY CONDITION AT X = XU
%...                    (UX(N) IS NOT USED)
%...
%...                2 - NEUMANN BOUNDARY CONDITION AT X = XU
%...                    (UX(N) IS USED)
%...
%...  GRID SPACING
      dx=(xu-xl)/(n-1);
%...
%...  CALCULATE UXX AT THE LEFT BOUNDARY, WITHOUT UX
      if nl==1
      uxx(1)=((     2.)*u(  1)...
             +(    -5.)*u(  2)...
             +(     4.)*u(  3)...
             +(    -1.)*u(  4))/(dx^2);
%...
%...  CALCULATE UXX AT THE LEFT BOUNDARY, INCLUDING UX
      elseif nl==2
      uxx(1)=((    -7.)*u(  1)...
             +(     8.)*u(  2)...
             +(    -1.)*u(  3))/(2.*dx^2)...
             +(    -6.)*ux( 1) /(2.*dx);
      end
%...
%...  CALCULATE UXX AT THE RIGHT BOUNDARY, WITHOUT UX
      if nu==1
      uxx(n)=((     2.)*u(n  )...
             +(    -5.)*u(n-1)...
             +(     4.)*u(n-2)...
             +(    -1.)*u(n-3))/(dx^2);
%...
%...  CALCULATE UXX AT THE RIGHT BOUNDARY, INCLUDING UX
      elseif nu==2
      uxx(n)=((    -7.)*u(n  )...
             +(     8.)*u(n-1)...
             +(    -1.)*u(n-2))/(2.*dx^2)...
             +(     6.)*ux(n ) /(2.*dx);
      end
%...
%...  CALCULATE UXX AT THE INTERIOR GRID POINTS
      for i=2:n-1
         uxx(i)=(u(i+1)-2.*u(i)+u(i-1))/dx^2;
      end

\end{verbatim}

\vspace{0.5cm}

Similar routines are available for fourth, sixth, eighth and tenth order
approximations of a second derivative, including Dirichlet and Neumann
boundary conditions at the boundaries.

\vspace{0.5cm}

We now all of the spatial derivative approximations required to convert eq.
(1) into a system of ODEs (using the approach of the MOL)%

\begin{equation}
\frac{du_{i}}{dt}=(\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta r^{2}})+\frac{1}%
{r_{i}}(\frac{u_{i+1}-u_{i-1}}{2\Delta r})\tag{7}%
\end{equation}

where $i=2,3,\cdots,N-1$. For $i=1$ and $i=N$, we use noncentered
approximations so that fictitious points are avoided. Also, we have to provide
two BCs in the programming of the ODEs, one at $i=1$ and one at $i=N$. The
preprogramming of the FD approximations for first and second derivatives in
library routines (such as dss002.m and dss042.m) makes the implementation of
BCs straightforward.

\vspace{0.5cm}

We can use more grid points in the approximations, buy starting with more
Taylor series for neighboring grid points, to arrive at higher order FD
approximations. For example, if we use five points (rather than three) in the
FD approximations, the differentiation matrix is (note: the spatial variable
is now $z$ rather than $r$)

\vspace{0.5cm}%

\[
\frac{d\mathbf{u}}{dz}=\left[
\begin{array}
[c]{ccccc}%
-50 & 96 & -72 & 32 & -6\\
-6 & -20 & 36 & -12 & 2\\
2 & -16 & 0 & 16 & -2\\
-2 & 12 & -36 & 20 & 6\\
6 & -32 & 72 & -96 & 50
\end{array}
\right]  \mathbf{u}+O(\Delta z^{4})
\]

\vspace{0.5cm}

We can now assign various combinations of grid points to the rows. For
example, if we use

\vspace{0.5cm}%

\[%
\begin{array}
[c]{cc}%
\text{Grid point(s)} & \text{Weighting coefficients}\\
1 & -50,96,-72,32,-6\\
2 & -6,-20,36,-12,2\\
3,\cdots N-2 & 2,-16,0,16,-2\\
N-1 & -2,12,-36,20,6\\
N & 6,-32,72,-96,50
\end{array}
\]

then we have used primarily centered approximations, which would be
appropriate for diffusion problems, such as for the first derivative in eq. (1).

\vspace{0.5cm}

If we use%

\[%
\begin{array}
[c]{cc}%
\text{Grid point(s)} & \text{Weighting coefficients}\\
1 & -50,96,-72,32,-6\\
2 & -6,-20,36,-12,2\\
3 & 2,-16,0,16,-2\\
4,\cdots,N-1 & -2,12,-36,20,6\\
N & 6,-32,72,-96,50
\end{array}
\]

these are the biased upwind approximations in dss020.m discussed in the heat
exchanger example. In other words, the differentiation matrix is the basis for
a spectrum of FD approximations (and their coding follows directly from the
differentiation matrix).

\vspace{0.5cm}

Finally, all of the preceding discussion has been for equally spaced grids. If
the grid is unequally spaced, for example, to concentrate the grid points in
regions of large spatial variation of the solution, analogous approximations
can be derived which have the same mathematical form, i.e., weighted sums of
the dependent variable at neighboring points. The derivation of these
nonuniform approximations is a bit more involved (the Taylor series is not as
easy to use for this case, and the starting point is therefore usually the
Lagrange interpolation polynomial).

\vspace{0.5cm}

Fortunately, all of this has been done by Bengt Fornberg in a very compact
algorithm that is programmed in a variety of languages, e.g., Fortran,
Mathematica, Maple: Fornberg, B., \textit{Calculation of Weights in Finite
Difference Formulas}, SIAM Review, v. 40, no. 3, pp 685-691, Sept., 1998. The
small Fortran subroutine, weights.for, that accompanies this algorithm is
widely used in PDE numerical analysis. It can produce the FD weighting
coefficients for

\vspace{0.5cm}

\begin{itemize}
\item  Derivatives of any order $m$ ($\dfrac{d^{m}}{dz^{m}}$)

\item  Approximations of any order $p$ ($O(\Delta z^{p})$)

\item  Approximations at any grid point, $i=1,...,n$ (centered and noncentered
to avoid fictitious points)

\item  Uniform and nonuniform grids ($\Delta z$ does not have to be constant),
so that the grid points can be concentrated where the solution changes rapidly
in $z$; further if the code does this automatically, we have adaptive mesh
refinement (\textit{r refinement} or \textit{amr}).
\end{itemize}

\vspace{0.5cm}

To conclude, we select approximations based on the general characteristics of
the PDEs (hyperbolic, parabolic, elliptic). Here are some general guidelines

\vspace{0.5cm}

\begin{itemize}
\item  Centered approximations are appropriate for elliptic and parabolic PDEs
since there is no preferred direction

\item  Parabolic PDEs will generally ``diffuse'' away discontinuities so that,
for example, incompatibilites between initial and boundary conditions will not
cause numerical difficulties

\item  Elliptic problems, which have no initial value derivatives, can be
accommodated by converting them to parabolic PDEs (by adding pseudo time
derivatives, then integrating until these derivatives essentially go to zero)

\item  Hyperbolic PDEs are the most difficult to solve numerically because
they propagate discontinuities

\item  Centered approximations for hyperbolic problem generally don't work
(they usually produce unrealistic oscillations)

\item  Noncentered approximations must be used for hyperbolic PDEs that take
into account the preferred direction (the direction of flow); in general,
hyperbolic PDEs require some form of upwinding

\item  Thus, we have to know the direction of flow so we can then select the
direction for the upwinding

\item  Using higher order approximations does not guarantee better accuracy;
for example, higher order centered approximations applied to hyperbolic
problems actually produces more oscillation

\item  The choice of an approximation(s), particularly for a high order,
nonlinear PDE problem, is as much art as science; thus we should be aware of
and consider the use of a variety of approximations, e.g., FDs, FVs, FEs,
least squares, flux limiters, ENO (essentially nonoscillatory), TVD (total
variation diminishing), etc. approximations

\item  These approximations originated in the CFD field, primarily to minimize
numerical distortion and error such as numerical diffusion and oscillation;
typically, CFD applicatiuons require the solution of strongly hyperbolic PDEs
systems, e.g., the Euler and Navier Stokes equations

\item  Numerical methods for PDEs is a very active field of research, with new
methods and associated software continuously being announced

\item  Much of the impetus for new, high speed (supercomputers) comes from
applications requiring the solution of systems of nonlinear, multidimensional,
time dependent PDEs, e.g., CFD, aircraft design, geophysics, atmospheric
modeling, numerical general relativity.
\end{itemize}

\vspace{0.5cm}

We now consider briefly some other forms of spatial discretization.

\vspace{0.5cm}

\textbf{Weighted Residuals}

\vspace{0.5cm}

The method of weighted residuals requires more analytical work than the FD
approach to PDEs, and the analytical work is rather specific for the problem
at hand (and thus, generally has to be repeated for each new problem). To give
an idea of this approach, if we consider Fourier's second law in Cartesian coordinates%

\begin{equation}
\frac{\partial u}{\partial t}=\frac{\partial^{2}u}{\partial x^{2}}\tag{1}%
\end{equation}%

\begin{equation}
u(x_{0},t)=0,\frac{\partial u(x_{L},t)}{\partial x}=0\tag{2)(3}%
\end{equation}%

\begin{equation}
u(x,0)=1\tag{4}%
\end{equation}

we assume a series solution of the form%

\begin{equation}
u(x,t)\cong\sum_{i=1}^{N}c_{i}(t)\phi_{i}(x)\tag{5}%
\end{equation}

Note that the RHS of eq. (5) is a \textit{product} or \textit{separated,
}which makes substituion into the problem equations relatively easy.

\vspace{0.5cm}

For example, if we compute the derivatives in eq. (1),%

\[
\frac{\partial u}{\partial t}\cong\sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x)
\]
\[
\frac{\partial^{2}u}{\partial x^{2}}\cong\sum_{i=1}^{N}c_{i}(t)\phi
_{i}^{\prime\prime}(x)
\]

and then substitute these derivatives in eq. (1), we have%

\begin{equation}
\sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x)=\sum_{i=1}^{N}c_{i}(t)\phi
_{i}^{\prime\prime}(x)\tag{6}%
\end{equation}

Since eq. (5) is only an approximate solution, it generally will have an error
or a \textit{residual}, $R(x,t).$ Eq. (6) can be written in residual form as%

\begin{equation}
\sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x)-\sum_{i=1}^{N}c_{i}(t)\phi
_{i}^{\prime\prime}(x)=R(x,t)\tag{7}%
\end{equation}

If we could make the residual zero everywhere, we would have the exact
solution. Generally, this is not possible, so we try to approximate this
condition by applying a weighting function, $w(x)$, to the residual according
to the integral%

\begin{equation}
\int_{x_{l}}^{x_{u}}w(x)R(x,t)dx=0\tag{8}%
\end{equation}

This is the basis of the \textit{method of weighted residuals}. Several well
known methods then follow from eq. (7), depending on the form of the weighting function.

\vspace{0.5cm}

\begin{itemize}
\item  Setting the residual to zero at a series of values of $x$ is the
\textit{collocation method.} The expectation, then, is that the residual will
be small at points in between those points where it is set to zero. Fort
example, if we select points $x_{l}<x_{1}<x_{2}<x_{3}<x_{u}$ as collocation
points, then from eq. (7), with $N=3$

\vspace{0.5cm}%

\[
\sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x_{1})-\sum_{i=1}^{N}c_{i}(t)\phi
_{i}^{\prime\prime}(x_{1})=0
\]
\[
\sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x_{2})-\sum_{i=1}^{N}c_{i}(t)\phi
_{i}^{\prime\prime}(x_{2})=0
\]%

\[
\sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x_{3})-\sum_{i=1}^{N}c_{i}(t)\phi
_{i}^{\prime\prime}(x_{3})=0
\]

The first equation is%

\[
c_{1}^{\prime}(t)\phi_{1}(x_{1})+c_{2}^{\prime}(t)\phi_{2}(x_{1}%
)+c_{3}^{\prime}(t)\phi_{3}(x_{1})-c_{1}(t)\phi_{1}^{\prime\prime}%
(x_{1})-c_{2}(t)\phi_{2}^{\prime\prime}(x_{1})-c_{3}(t)\phi_{3}^{\prime\prime
}(x_{1})=0
\]

This is a linearly implicit ODE with coupling coefficients $\phi_{1}%
(x_{1}),\phi_{2}(x_{1}),\phi_{3}(x_{1})$. This, we can use the Cash code
discussed previously for the $3\times3$ problem.

\vspace{0.5cm}

The BCs are accommodated according to the choice of the \textit{basis
functions}, $\phi_{i}(x)$. Typically, these basis functions are chosen to
identically satisfy the boundary conditions. If the collocation points are
chosen as the roots of an orthogonal polynomial, we then have
\textit{orthogonal collocation} (OC). This choice of the collocation points
gives some error minimizing qualities to the solution.

\vspace{0.5cm}

The weighting function for collocation in eq. (8) can be considered a delta
function, i.e.,%

\[
\int_{x_{l}}^{x_{u}}\delta(x-x_{1})R(x,t)dx=R(x_{1},t)=0
\]

so the residual is set to zero at $x=x_{1}$ as above.

\vspace{0.5cm}

\item  If the weighting function is taken as a basis function,%

\begin{equation}
\int_{xl}^{x_{u}}\phi_{i}(x)R(x,t)dx=0\tag{10}%
\end{equation}

This is the \textit{Galerkin method}. The essential idea is that the
\textit{basis functions are orthogonal to the residual}. Again, this has
certain error minimizing qualities for the solution of the PDE through eq. (5)
(it does not make the residual zero for all $x$ but rather, makes it small, on
average over the interval $x_{l}\leq x\leq x_{u}$). If the Galerkin criterion
(10) is applied to the $N$ weighting functions, we again arrive at linearly
implicit ODEs (or DAEs depending on how the BCs are handled).
\end{itemize}

\vspace{0.5cm}

\textbf{Finite Elements}

\vspace{0.5cm}

Finite elements (FE) is a vast subject, and we can only mention here a few
basic ideas. If in the method of weighted residuals as summarized by eq. (8),
we choice the basis functions, $\phi_{i}(x)$, to have \textit{compact
support}, that is, to have nonzero values only over a finite interval in $x$,
these ``elements'' then serve as the basis of the assumed solution such as eq.
(5). The elements can have any mathematical form (and many have been proposed
and used); a common choice is the so-called ``hat function'' or ``chapeau
function'', which makes the evaluation of the integrals in the FE method
relatively easy. Thus, if we substitute the residual for the preceding PDE
problem in eq. (10)%

\[
\int_{x_{l}}^{x_{u}}\phi_{j}(x)R(x,t)dx=\int_{x_{l}}^{x_{u}}\phi
_{j}(x)\left\{  \sum_{i=1}^{N}c_{i}^{\prime}(t)\phi_{i}(x)-\sum_{i=1}^{N}%
c_{i}(t)\phi_{i}^{\prime\prime}(x)\right\}  dx=0
\]

Note now that when the integration in $x$ is performed (e.g., by interchanging
the orders of integration and summation, and using integration by parts), the
result is a function of $t$ only, and we again arrive at linearly implicit
ODEs. For example, using hat functions as the FEs, we obtain the implicit ODE%

\begin{equation}
(1/6)c_{i+1}^{\prime}+(4/6)c_{i}^{\prime}+(1/6)c_{i-1}^{\prime}=D\left\{
\frac{c_{i+1}-2c_{i}+c_{i-1}}{\Delta x^{2}}\right\} \tag{11}%
\end{equation}

Eq. (11) is of the same form as the FD approximation of eq. (1) using second
order, central FDs. Note, however, that in using the FE method, we

\vspace{0.5cm}

\begin{itemize}
\item  Had to do a substantial amount of analytical work, even for the modest
problem of eqs. (1) to (4).

\item  Used FEs which are only linear (rather than the higher order
polynomials that are the basis of the higher order FDs)
\end{itemize}

\vspace{0.5cm}

Thus, we might expect we would have to use a relatively large number of FEs
(and therefore grid points and ODEs) to achieve accuracy in the PDE solution
comparable to a higher order FD. These conclusions also apply to orthgonal
collocation (OC is relatively difficult to apply, and its performance is
probably not any better than when using carefully selected FDs).

\vspace{0.5cm}

The principal utility of the FE method is the spatial discretization on
irregular geometries, for which the FD method becomes diffiicult to apply.
Thus, body fitted coordinates are a very active area of research, and the grid
that results from this analysis (based on triangles, for example), is then
used within the FE framework. Computer codes that accommodate irregular
geometries are therefore generally based on the FE method, and a key
consideration in designing and using these codes is the choice of the basic
finite element shapes. A spectrum of finite elements has been developed and
analyzed, particularly for 2- and 3-D PDE systems. This is the basis of the
Matlab PDE toolbox, for example (it is a 2-D FE code).

\vspace{0.5cm}

\textbf{Finite Volumes}

\vspace{0.5cm}

Another approach to spatial discretization (in addition to FDs and FEs) is the
finite volume (FV) formulation. Briefly, this method is based on the use of
conservation principles such as the conservation of mass and conservation of
energy, as applied to a finite or incremental volume of the system to be
analyzed. To illustrate this approach, we again consider the analysis of a
system modeled by the one-dimensional heat conduction equation, and we then
subdivide the spatial coordinate into a grid, corresponding to a series of
contiguous volumes. An energy balance wtitten on one of the finite volumes is then%

\begin{equation}
V_{i}\rho_{i}C_{pi}\frac{du_{i}}{dt}=A_{i-1}q_{i-1}-A_{i}q_{i}\tag{12}%
\end{equation}

which is interpreted to conserve eneregy exactly by the proper choice of the
fluxes, $q_{i}$ and $q_{i-1}.$

\vspace{0.5cm}

For example, the fluxes could be computed as%

\[
q_{i}=k_{i,i+1}(\frac{u_{i}-u_{i+1}}{\Delta x_{i}/2+\Delta x_{i+1}/2})
\]

where $u_{i}$ is interpreted as the mid-cell (or mid-volume) temperature.
Thus, we end up again with a system of ODEs, so that the FV method can be
considered as a variant of the MOL, with a different spatial discretization
(than FDs or FEs). The challenge, in general, in using the FV approach is to
calculate the intercell fluxes in such a way that conservation is preserved.
This becomes more difficult for multidimensions, irregular geometries and
coordinate systems other than Cartesian coordinates. Thus, a substantial
effort has recently been devoted to the FV method in cylindrical and spherical
coordinates so that conservation is maintained, e.g., in the solution of
Maxwell's equations for calculating EM fields in communications systems with
arrays (antennas) with non-Cartesian shapes.

\vspace{0.5cm}

In summary, the use of the FV method is a very active field, with the
attractive feature of maintaining conservation of basic physical quantities,
so that the discretizied (approximated) PDEs in fact maintain the conservation
implied by the PDEs themselves.

\vspace{0.5cm}

\textbf{Adaptive Grids}

\vspace{0.5cm}

The equidistribution principle uses a monitor function based on the curvature
of the solution.

\vspace{1cm}

At the time level $t_{k+1}$, the nodes $x_{i}^{k+1},i=1,\ldots,N$ are located
such that%

\[
\int_{x_{i-1}^{k+1}}^{x_{i}^{k+1}}m(u(x,t_{k+1}))dx=\int_{x_{i-1}^{k+1}%
}^{x_{i}^{k+1}}\left\{  \alpha+\left\|  \frac{\partial^{2}u(x,t_{k+1}%
)}{\partial x^{2}}\right\|  _{\infty}\right\}  ^{1/2}dx=\text{constant}
\]

where $u(x,t_{k+1})$ is the PDE solution which has been advanced to the time
$t_{k+1}$ using the fixed grid $x_{i}^{k},i=1,\ldots,N$

\vspace{0.5cm}

MOL Formulation:

\begin{itemize}
\item  Finite difference approximations (Fornberg)

\item  Time integration with the implicit BDF solver LSODI (Hindmarsh) or the
implicit RK solver RADAU5 (Hairer, Wanner)

\item  Parameters

\begin{itemize}
\item  Scaling factor $\alpha$

\item  Limiting factor $\beta$

\item  Updating frequency $N_{adapt}$
\end{itemize}

\item  Adaptive grids require parameter tuning \vspace{0.5cm}

\item  There is no universal method ... yet! \vspace{0.5cm}

\item  Our Fortran codes are generally available, e.g, subroutine AGE \vspace{0.5cm}

\item  The contributions of Dr. Alain Vande Wouwer and Dr. Philippe Saucez,
Polytechnique Mons (Belgium) are gratefully acknowledged
\end{itemize}

\vspace{0.5cm}

\textbf{References and Software}

\vspace{0.5cm}

Sources of mathematical software include:

\vspace{0.5cm}

Name: Transactions on Mathematical Software (TOMS), Association for Computing
Machinery (ACM)

\begin{itemize}
\item  Features: An extensive library of high quality, public domain
mathematical software

\item  Reference: http://www.acm.org/dl/search.html
\end{itemize}

\vspace{0.5cm}

Name: Netlib (Inter\textbf{net} \textbf{Lib}rary) - accessed 74,280,550 times
by 22APR00

\begin{itemize}
\item  Features: An extensive library of high quality, public domain PDE
mathematical software

\item  Reference: \newline http://www.netlib.org/liblist.html \newline Select
pde\newline \newline netlib@ornl.gov (whois ...., send index, send index from
toms) \newline 
\end{itemize}

\vspace{0.5cm}

Name: NA-Digest (\textbf{N}umerical \textbf{A}nalysis Digest; an electronic newsletter)

\begin{itemize}
\item  Features: A free, weekly newsletter edited by Dr. Cleve Moler, the
author of Matlab; to subscribe and get information:\newline
na.digest@na-net.ornl.gov\newline na.whois@na-net.ornl.gov\newline na.help@na-net.ornl.gov

\item  Reference: All past issues of NA-Digest are in Netlib (and they can be
searched by subject and author)

\item  Example:

From: DEAL
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
deal@hermes.iwr.uni-heidelberg.de%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
\newline Date: Tue, 18 Apr 2000 10:42:32 +0200 (MET DST)\newline Subject:
DEAL, C++ Finite Element Library

\vspace{0.5cm}

Version 3.0 of the deal.II object-oriented finite element library is available
on the deal.II home-page at

\vspace{0.5cm}

http://gaia.iwr.uni-heidelberg.de

\vspace{0.5cm}

deal.II is a C++ program library targeted at adaptive finite elements and
error estimation. It uses state-of-the-art programming techniques of the C++
programming language to offer you a modern interface to the complex data
structures and algorithms required for adaptivity and enables you to use a
variety of finite elements in one, two, and three space dimensions, as well as
support for time-dependent problems.

\vspace{0.5cm}

The library is written for research purposes and offers many features:

\vspace{0.5cm}

- Support for one, two, and three space dimensions, using a unified interface
that enables writing programs almost dimension independent.

\vspace{0.5cm}

- Handling of locally refined grids, including different adaptive refinement
strategies based on local error indicators and error estimators.

\vspace{0.5cm}

- Support for a variety of finite elements, including Lagrange elements of
order one through four, and discontinuous elements.

\vspace{0.5cm}

- 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 about 200 pages of tutorials, several reports, and far more than
1,000 pages of programming interface documentation with explanations of all
classes, functions, and variables.

\vspace{0.5cm}

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

\vspace{0.5cm}

- Fast algorithms that enable you to solve problems with up to several
millions of degrees of freedom quickly. As opposed to programming symbolic
algebra packages the penalty for readability is low.

\vspace{0.5cm}

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

\vspace{0.5cm}

- Support for a variety of computer platforms, including multi- processor machines.

\vspace{0.5cm}

- Free source code under an Open Source license, and the invitation to
contribute to further development of the library.

\vspace{0.5cm}

Wolfgang Bangerth, Guido Kanschat, the deal.II team

\vspace{0.5cm}

From: Tom Goodale
%TCIMACRO{\TEXTsymbol{<}}%
%BeginExpansion
$<$%
%EndExpansion
goodale@aei-potsdam.mpg.de%
%TCIMACRO{\TEXTsymbol{>}}%
%BeginExpansion
$>$%
%EndExpansion
\newline Date: Thu, 20 Apr 2000 16:13:55 +0200 (CEST)\newline Subject: Cactus,
Parallel Code for PDEs

\vspace{0.5cm}

This is to announce the release of beta 7 of Cactus 4.0.

\vspace{0.5cm}

Cactus is a general, modular, parallel code for solving systems of partial
differential equations. The code has been developmented over many years by a
large international collaboration of numerical relativity and computational
science research groups and can be used to provide a portable platform for
solving any system of partial differential equations. The code compiles and
runs on a variety of different platforms, from laptops running Linux or NT to
clusters of workstations, to large T3Es, with no modifications needing to be
made to the application codes.

\vspace{0.5cm}

As an example, one set of modules (thorns) can be used to solve problems in
numerical relativity, for example black hole or neutron star collisions. The
design of the code and the thorns means that general routines can be written,
for example to locate the position of the event horizon, which can then be
made generally available and used for other problems.

\vspace{0.5cm}

The code is available from

\vspace{0.5cm}

http://www.cactuscode.org

\vspace{0.5cm}

Notes for this release can be found at

\vspace{0.5cm}

http://www.cactuscode.org/Development/beta7.txt

\vspace{0.5cm}

For further information, please email

\vspace{0.5cm}

cactusmaint@cactuscode.org
\end{itemize}

\vspace{0.5cm}

The NAG library has a spectrum of PDE software, much of it developed by Martin
Berzins, University of Leeds

\vspace{0.5cm}

Chapter D03 - Partial Differential Equations \newline D03EAF - Elliptic PDE,
Laplace's equation, 2-D arbitrary domain \newline D03EBF - Elliptic PDE,
solution of finite difference equations by SIP, five-point 2-D molecule,
iterate to convergence \newline D03ECF - Elliptic PDE, solution of finite
difference equations by SIP for seven- point 3-D molecule, iterate to
convergence \newline D03EDF - Elliptic PDE, solution of finite difference
equations by a multigrid technique \newline D03EEF - Discretize a 2nd order
elliptic PDE on a rectangle \newline D03FAF - Elliptic PDE, Helmholtz
equation, 3-D Cartesian co-ordinates \newline D03MAF - Triangulation of a
plane region \newline D03PCF - General system of parabolic PDEs, method of
lines, finite differences, one space variable \newline D03PDF - General system
of parabolic PDEs, method of lines, Chebyshev C0 collocation, one space
variable \newline D03PEF - General system of first order PDEs, method of
lines, Keller box discretisation, one space variable \newline D03PFF - General
system of PDEs, convection-diffusion in conservative form, method of lines,
one space variable \newline D03PHF - General system of parabolic PDEs, coupled
DAEs, method of lines, finite differences, one space variable \newline D03PLF
- PDEs, general system, one space variable, method of lines, compute
interpolant in spatial direction \newline D03PJF - General system of parabolic
PDEs, coupled DAEs, method of lines, Chebyshev C0 collocation, one space
variable \newline D03PKF - General system of first order PDEs, coupled DAEs,
method of lines, Keller box discretisation, one space variable \newline D03PPF
- General system of parabolic PDEs, coupled DAEs, method of lines, finite
differences, remeshing, one space variable \newline D03PRF - General system of
first order PDEs, coupled DAEs, method of lines, Keller box discretisation,
remeshing, one space variable \newline D03PSF - General system of PDEs,
convection-diffusion in conservative form, method of lines, coupled DAEs,
remeshing, comprehensive one space variable \newline D03PUF - General system
of PDEs, appropriate Riemann solver for Euler equations, method of lines,
Roe's scheme, one space variable \newline D03PVF - General system of PDEs,
appropriate Riemann solver for Euler equations, method of lines, Osher's
scheme, one space variable \newline D03PWF - Modified HLL Riemann solver for
Euler equations in conservative form, for use with D03PFF, D03PLF and D03PSF
\newline D03PXF - Exact Riemann Solver for Euler equations in conservative
form, for use with D03PFF, D03PLF and D03PSF \newline D03PYF - PDEs, spatial
interpolation with D03PDF or D03PJF \newline D03PZF - PDEs, spatial
interpolation with D03PCF or D03PHF \newline D03RAF - General system of second
order PDEs, method of lines, finite differences, remeshing, two space
variables, rectangular region \newline D03RBF - General system of second order
PDEs, method of lines, finite differences, remeshing, two space variables,
rectilinear region \newline D03RYF - Check initial grid data in D03RBF
\newline D03RZF - Extract grid data from D03RBF \newline D03UAF - Elliptic
PDE, solution of finite difference equations by SIP, five-point 2-D molecule,
one iteration \newline D03UBF - Elliptic PDE, solution of finite difference
equations by SIP, seven-point 3-D molecule, one iteration Chapter D04 -
Numerical Differentiation \newline D04AAF - Numerical differentiation,
derivatives up to order 14, function of one real variable
\end{document}
