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


\documentclass{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TCIDATA{Created=Sat Jan 22 09:53:41 2000}
%TCIDATA{LastRevised=Sun Jan 23 13:40:24 2000}

\parindent0cm
\input{tcilatex}
\textheight220mm\textwidth160mm
\oddsidemargin=-2pt
\evensidemargin=-2pt
\topmargin-2cm

\begin{document}

%\large

\begin{center}
{\Huge \textbf{Introduction to Scientific Computation}}
\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)\\
(610) 758-5057 (fax)\\

\vspace{0.5cm}

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

\newpage

\begin{center}
{\Huge \textbf{Course Orientation and Content}}\\

\vspace{0.2cm}

http://www.lehigh.edu/\~{\thinspace}wes1/apci/28jan00.tex
\end{center}

\vspace{0.5cm}

\begin{enumerate}
\item  Scientific computation

\item  Quality software and software libraries

\item  Programming languages

\item  Numerical considerations

\begin{enumerate}
\item  Computer word length

\item  Roundoff and truncation

\item  Conditioning

\item  Convergence and error tolerances
\end{enumerate}

\item  Interpolation

\begin{enumerate}
\item  Polynomials

\item  Splines
\end{enumerate}

\item  Integration

\begin{enumerate}
\item  Gaussian quadrature

\item  Splines

\item  Adaptive
\end{enumerate}

\item  Differentiation

\begin{enumerate}
\item  Finite differences

\item  Symbolic computation
\end{enumerate}

\item  References and software
\end{enumerate}

\vspace{0.5cm}

In this course, we will emphasize:

\begin{itemize}

\item  Terminology

\item  Modeling of engineering systems

\item  Mathematical forms resulting from the modeling

\item  Established algorithms (not the detailed
underlying mathematics, but the principal conclusions)

\item  Computational resources (e.g., software libraries)
for the application of the algorithms to the mathematical
forms

\end{itemize}

\newpage

Why LaTeX (.tex files)?

\vspace{0.5cm}

Latex:

\begin{itemize}

\item  Has excellent mathematical formatting
capabilities (better than the equation editor of
MS Word, for example).

\item  Is entirely ASCII-based, and therefore LaTeX files
can be composed and read with an ASCII editor and can easily be
transmitted over the Internet.

\item  Is in the public domain.

\item  Is the established text formatter for scientific
publishing (rather than MS Word, for example).  Most major
scientific publishers now specify LaTeX files as the preferred
format for publication.

\item  Produces relatively small files (the file for all of
the notes to be discussed today is under 70,000 bytes).

\end{itemize}

As an incidental point, LaTex code is easily written using
Scientific Workplace (SW);
SW also has a full Maple interface, and the results of Maple
computations can be embedded in LaTeX documents.

\newpage

\begin{center}
{\Huge \textbf{Scientific Computation}}
\end{center}

\vspace{0.5cm}

Scientific computing is the collection of tools, techniques
and theories required to solve on a computer mathematical
models of problems in science and engineering.

\vspace{0.5cm}

The case for scientific computation:

\vspace{0.5cm}

\begin{itemize}
\item  Reduced time

\vspace{0.5cm}

\item  Reduced cost

\vspace{0.5cm}

\item  Increasd spectrum of problem systems and parameters

\vspace{0.5cm}

\item  Enhancement of experimental results

\vspace{0.5cm}

\begin{itemize}
\item  Framework for interpretation of experimental data

\vspace{0.5cm}

\item  Evaluation of model parameters
\end{itemize}
\end{itemize}

\vspace{1.0cm}
Golub and Ortega (1992)

\newpage

\begin{center}
{\Huge \textbf{Unique Features of Scientific Computation}}
\end{center}

\vspace{0.5cm}

Results cannot be quaranteed in advance

\vspace{0.5cm}

Example 1:

\vspace{0.5cm}

Consider the $2x2$ linear algebraic system:

\[
\begin{array}{c}
1.5x_{1}+3.5x_{2}=2 \\
3.0x_{1}+7.0x_{2}=5
\end{array}
\]

To eliminate $x_{2}$, multiply the first equation by $2$ and subtract from
the second equation

\[
0=1
\]

A rather astounding result! How did this happen?

\vspace{0.5cm}

The coefficient matrix is singular, i.e.,

\[
\left|
\begin{array}{cc}
1.5 & 3.5 \\
3.0 & 7.0
\end{array}
\right| =(1.5)(7.0)-(3.0)(3.5)=0
\]

Note that

\begin{center}
2(row 1) = row 2
\end{center}

i.e., the two rows are \textit{linearly dependent}.

\vspace{0.5cm}

Thus, the $2x2$ system is \textit{singular}. Addtionally, because of the
meaningless result ($0=1$), it is \textit{inconsistent} (does not have a
solution).

\begin{itemize}
\item  Another closely related system is
\end{itemize}

\[
\begin{array}{c}
1.5x_{1}+3.5x_{2}=2.5 \\
3.0x_{1}+7.0x_{2}=5
\end{array}
\]

which leads to

\[
0=0
\]

Thus, we at least arrive at something that makes sense mathematically.
However, the coefficient matrix is still singular so that we cannot solve
for unique values of $x_{1}$ and $x_{2}$. This system is therefore termed
\textit{consistent} but \textit{singular.}

\vspace{0.5cm}

Since these equations are essentially the same, we can really only use one
or the other. Thus, we in effect have one equation and two unknowns, and we
therefore must select one unknown arbitrarily, then solve for the other
unknown. A (nonunique) solution to this system is

\[
\begin{array}{c}
x_{1}=c \\
x_{2}=(2.5-1.5c)/3.5
\end{array}
\]

where $c$ is an arbitrary constant.

\begin{itemize}
\item  More typical situation:
\end{itemize}

\[
\begin{array}{c}
1.499999x_{1}+3.5x_{2}=2 \\
3.0x_{1}+7.0x_{2}=5
\end{array}
\]

Solving for $x_{1}$ by eliminating $x_{2}$,

\[
(2)(1.499999)x_{1}-3.0x_{1}=1
\]

or

\[
(2.999998-3.0)x_{1}=1
\]

Our computer had better carry at least seven figures in its calculations.
But this is all that single precision Fortran, for example, can do (with
32-bit arithmetic).

\vspace{0.5cm}

The system is termed \textit{nearly singular}, \textit{ill-conditioned} \ or
\textit{numerically singular} with 32-bit arithmetic.

\vspace{0.5cm}

Thus, the precision determined by the word length of the computer in
scientific calculations can be very important, and we should have some
quantitative measure of the precision of any new computer that we use (to be
addressed subsequently).

\vspace{0.5cm}

For example, with single precision Fortran, if the first coefficient had
been $1.49999999$, the problem would have been singular, and a solution
could not be calculated. However, with double precision Fortran (64-bit
arithmetic), we could still expect a solution with an accuracy of several
(i.e., $4-5$) significant figures.

\begin{itemize}
\item  We conclude that depending on the characteristics of the problem, and
the precision of the computer arithmetic, we may not be able to compute a
solution to even a $2x2$ linear algebraic system.

\item  The problems with these small ($2x2$) systems were rather obvious.
But similar problems with $1000x1000$ systems, for example, will not be so
obvious, yet can be just as troublesome.

\item  More generally, we may not be able to produce a ``solution'' because
there is none (the problem is not properly posed mathematically), or the
mathematical model is so ill-conditioned that a numerical solution is
precluded.

\item  Ideally, the software that we develop and use should include some
safeguards against such potential conceptual and numerical problems, but in
general we cannot count on this. In other words, there are no guarantees of
a useful result in scientific computation.

\item  Typically, by a process of trial and error, we develop a numerical
method that is efficient and accurate, but this is an experimental process,
and the success of this effort is directly tied to the characteristics of
the problem at hand.

\item  This contrasts with calculations that never fail (e.g., financial
calculations). Also, having a faster computer might help, but it does not
guarantee success. Understanding the source of the computational
difficulties, and choosing an effective algorithm is more important than raw
speed.

\item  Fortunately, the iterative process in scientific computation
generally leads to a successful result, but this is a function of the
experience and insight of the analyst.
\end{itemize}

\newpage

Example 2:

\vspace{0.5cm}

Consider the third order polynomial

\[
f(x)=x^{3}-2x-5=0
\]

This polynomial may be the first application of Newton's method, since
Newton used it to illustrate his method (with a geometric argument, but not
calculus!).

\vspace{0.5cm}

If we apply Newton's method

\[
x^{i+1}=x^{i}-\frac{f(x^{i})}{df(x^{i})/dx},i=0,1,\cdots
\]

(where $i$ is an iteration index or counter) and we make the choice of the
initial starting point a root of

\[
3x^{2}-2=0
\]

or

\[
x^{0}=\pm \sqrt{2/3}
\]

the method will fail because this is a \textit{singular point} (the \textit{%
\ Jacobian matrix}, or in this case, the first derivative, is singular).
Also, if we choose a starting point near the singular point, the method will
most likely fail, i.e., the system is \textit{near-singular}.

\vspace{0.5cm}

Thus, we face the same difficulties as with the linear $(2x2)$ system. The
choice of an initial point that is not singular or near-singular becomes
increasingly difficult as the size of the problem increases (more
equations), and we should have a test in our software for the \textit{%
condition number} of the Jacobian matrix.

\vspace{0.5cm}

\begin{itemize}
\item  In conclusion, mathematical computation, in general, is not a
completely straightforward process.

\item  Additionally, the mathematical forms are numerous, e.g.,

\begin{itemize}
\item  Linear algebra and eigenvalue problems

\item  Nonlinear algebraic and transcendental equations

\item  1-D and multidimensional integrals

\item  Ordinary differential equations (ODEs)

\item  Partial differential equations (PDEs)

\item  Optimization (e.g., max/min) problems
\end{itemize}

\item  All of these mathematical forms require some knowledge of at least
basic algorithms and software.

\item  The field of scientific computation is large, challenging and rapidly
moving; not surprisingly, professional numerical analysts devote full time
effort to keeping abreast of this field.
\end{itemize}

\vspace{0.5cm}

How do we arrive at these mathematical forms when applying mathematical
analysis of engineering systems? Here's an example:

\vspace{0.5cm}

A tubular reactor, as depicted in the accompanying diagram, 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 analayis 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\rightarrow 0,\Delta z\rightarrow 0$,

\[
\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\rightarrow 0,\Delta z\rightarrow 0$,

\[
\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}{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,375 \\
\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\times 10^{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) requires 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 interested 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}

We will consider the numerical solution of this problem later in the section
on PDEs, and we will observe that with $T_{w}=305^{o}K$, essentially no
reaction takes place. If the wall is heated so that $T_{w}=375^{o}K$,
significant reaction takes place, but there are no hotspots or a runaway
reaction. Thus, the reactor should give a reasonable production of the
product if the wall is heated, not cooled, as originally suggested.

\vspace{0.5cm}

The validity of these conclusions, of course, depends on the validity of the
mathematical model and the numerical values of the parameters. But assuming
these requirements are satisfied, we can gain valuable insights and useful
information about the performance of the reactor from the computer-based
study.

\vspace{0.5cm}

``Off-the-shelf'' software for the solution of eqs. (1) to (13) is not
available. Rather, we have to construct a code, perhaps using subroutines
that will calculate terms and special parts of this set of equations. For
example, we could use library routines to calculate the partial derivatives
in the PDEs, and the integrals in eqs. (12) and (13). Thus, software
construction and testing must, in general, precede the analysis of the
problem system.

\newpage

\vspace{0.5cm}

\begin{center}
\textbf{\Huge Quality Software and Software Libraries}
\end{center}

\vspace{0.5cm}

Quality mathematical software is readily available with the following
properties:

\begin{itemize}
\item  Based on algorithms proven to be effective through use over time

\item  Written to include special cases and identify problems that might not
occur to an analyst without the experience of the software author(s)

\item  Tranportable to most computers, and easily called, typically with a
one line command or a single parameter in a call to a compiler

\item  Thoroughly tested and documented
\end{itemize}

There is also the consideration of software based on the use of a compiler
with a national or international standard vs a commercial product whose
existence is dependent on a single company.

\vspace{0.5cm}

The standard compilers for mathematical software include:

\begin{itemize}
\item  Fortran 77/90/95/2000

\item  C

\item  C++ Includes object oriented programming (OOP)

\item  Java and extensions
\end{itemize}

Commercial products include:

\begin{itemize}
\item  Maple

\item  Mathcad

\item  Mathematica

\item  Matlab

\item  Visual Fortran

\item  Visual Basic (VB)
\end{itemize}

Note that Visual Fortran and Visual Basic are not languages.

\vspace{0.5cm}

Collections/Libraries of mathematical routines:

\begin{itemize}
\item  Forsythe et. al. (1977), Computer Methods for Computational
Mathematics 

\item  Kahaner et. al. (1989), Numerical Methods and Software 

\item  Press et. al. (1992), Numerical Recipes (Fortran, C)

\item  IMSL (Fortran, C)

\item  NAG (Fortran, C; available at Air Products)

\item  Netlib (ORNL, AT\&T, e.g., ODEPACK, DASSL)

\item  ACMS TOMS library

\item  NA Digest (sent weekly)

\end{itemize}

Also, there are combinations, e.g., VB/NAG using Fortran dll's.

\vspace{0.5cm}

As an example of what is available (taken from)

\begin{center}
http://www.nag.co.uk/numeric/FLOLCH/mk18.html
\end{center}

\vspace{0.5cm}

Fortran Library Mark 18 Introduction \\

\vspace{0.5cm}

Chapter A00 - Library Identification\\
Chapter A02 - Complex Arithmetic\\
Chapter C02 - Zeros of Polynomials\\
Chapter C05 - Roots of Transcendental Equations\\
Chapter C06 - Summation of Series\\
Chapter D05 - Integral Equations\\
Chapter E01 - Interpolation\\
Chapter E02 - Curve and Surface Fitting\\
Chapter E04 - Minimizing or Maximizing a Function\\
Chapter F01 - Matrix Factorizations\\
Chapter F02 - Eigenvalue and Eigenvectors\\
Chapter F03 - Determinants\\
Chapter F04 - Simultaneous Linear Equations\\
Chapter F05 - Orthogonalisation\\
Chapter F06 - Linear Algebra Support Routines\\
Chapter F07 - Linear Equations (LAPACK)\\
Chapter F08 - Least Squares and Eigenvalue Problems\\
Chapter F11 - Sparse Linear Algebra\\
Chapter G01 - Simple Calculations and Statistical Data\\
Chapter G02 - Correlation and Regression Analysis\\
Chapter G03 - Multivariate Methods\\
Chapter G04 - Analysis of Variance\\
Chapter G05 - Random Number Generators\\
Chapter G07 - Univariate Estimation\\
Chapter G08 - Nonparametric Statistics\\
Chapter G10 - Smoothing in Statistics\\
Chapter G11 - Contingency Table Analysis\\
Chapter G12 - Survival Analysis\\
Chapter G13 - Time Series Analysis\\
Chapter H - Operations Research\\
Chapter M01 - Sorting\\
Chapter P01 - Error Trapping\\
Chapter S - Approximations of Special Functions\\
Chapter X01 - Mathematical Constants\\
Chapter X02 - Machine Constants\\
Chapter X03 - Inner Products\\
Chapter X04 - Input/Output Utilities\\
Chapter X05 - Date and Time Utilities\\

\vspace{0.5cm}

The Numerical Algorithms Group Ltd, Oxford UK, 1999\\

\vspace{0.2cm}

Thus, there are many sources and combinations, and standardization is
therefore a problem in the sense that the analyst must select software items
that are compatible, or can be made compatible with reasonable effort.

\newpage

\begin{center}
\textbf{\Huge Programming Languages}
\end{center}

\vspace{0.5cm}

\begin{itemize}
\item  Fortran 77/90/95/2000

\item  C

\item  C++ Includes object oriented programming (OOP)

\item  Java and extensions
\end{itemize}

All have international standards that insure tranportability, but not
necessarily longevity. For example, what happened to PASCAL, APL, ALGOL, to
name a few? Longevity is an important issue in terms of preserving and
extending the investment in programming over time.

\vspace{0.5cm}

Fortran, for all of its purported limitations, has preserved the investment
made in it over an extraordinarily long time (for the computer field), and
it is evolving, and thus improving so that it's lifetime is unpredictable,
but certainly extended.

\vspace{0.5cm}

There is the basic question of what newer languages will do in scientific
computing that Fortran cannot do, especially when considering new releases
of Fortran.

\vspace{0.5cm}

The same questions of longevity and investments in programming can be
applied to reliance on a commercial product that depends on the viability of
one company.

\newpage

\begin{center}
\textbf{\Huge Numerical considerations}
\end{center}

\begin{enumerate}
\item  Computer word length

\item  Roundoff and truncation

\item  Conditioning

\item  Convergence and error tolerances
\end{enumerate}

We considered briefly the effect of computer word length on the solution of
a $2x2$ linear system. The system

\[
\begin{array}{c}
1.499999x_{1}+3.5x_{2}=2 \\
3.0x_{1}+7.0x_{2}=5
\end{array}
\]

is ill-conditioned if the computer carries seven significant figures
(corresponding to 32 bits). This conclusion follows from the subtraction $%
2.999998-3.0$ when solving for $x_1$.

The system

\[
\begin{array}{c}
1.49999999x_{1}+3.5x_{2}=2 \\
3.0x_{1}+7.0x_{2}=5
\end{array}
\]

is numerically singular for a computer with seven figure arithmetic.

\vspace{0.5cm}

In general, the precision of the numerical calculations as determined by the
computer word length is an important factor in determining the accuracy of
the final computed result. In fact, the single most important number that
characterizes the arithmetic performed by a specific computer with specifc
software is the machine epsilon or unit roundoff, defined as the smallest
number $eps$ such that

\vspace{0.5cm}

\[
1+eps>1
\]

To determine the value of the machine epsilon, consider the following small
Fortran program

\vspace{0.5cm}
\begin{verbatim}
      program epssp
      eps=1.0
      do while(1.0+eps.gt.1.0)
         eps=eps/2.0
      end do
      eps=2.0*eps
      write(*,1)eps
 1    format(/,' eps = ',e10.3)
      stop
      end
\end{verbatim}

\vspace{0.5cm}

The output from this program is eps = 0.119E-06. Thus, single precision
Fortran has a precision of about one part in $10^7$.

\vspace{0.5cm}

This might seem more than adequate for scientific and engineering
applications, e.g., rarely are experimental data better than one part in $%
10^4$. However, there are at least two reasons why better precision is
required:

\vspace{0.5cm}

\begin{itemize}
\item  In solving large sets of complicated equations, $%
10^{6},10^{9},10^{12} $ ... floating point operations are performed, each
with a small roundoff, and these errors can accumulate,
and may eventually invalidate the computed result. The best protection
against this accumulated error is to use higher precision arithmetic.

\item  The problem may be poorly conditioned as we observed, even with a $%
2x2 $ linear algebraic system (the equations in linear regression are well
known for being poorly conditioned). Again, the way around this problem may
be greater precision (although better algorithms will probably be more
effective).
\end{itemize}

\vspace{0.5cm}

Thus, if we consider the same program written in double precision Fortran,

\vspace{0.5cm}
\begin{verbatim}
      program epsdp
      implicit double precision (a-h,o-z)
      eps=1.0d0
      do while(1.0d0+eps.gt.1.0d0)
         eps=eps/2.0d0
      end do
      eps=2.0d0*eps
      write(*,1)eps
 1    format(/,' eps = ',d10.3)
      stop
      end
\end{verbatim}

\vspace{0.5cm}

The output from this program is eps = 0.222D-15. Thus, double precision
Fortran has a precision of about one part in $10^{16}$, which experience has
demonstrated is usually adequate for most applications in science and
engineering.

\vspace{0.5cm}

The follow program does the same calculation in Matlab:

\vspace{0.5cm}
\begin{verbatim}
      eps=1.0;
      while(1.0+eps>1.0)
         eps=eps/2.0;
      end
      eps=2.0*eps;
      eps
      end
\end{verbatim}

\vspace{0.5cm}

The output from this program is ans = 2.2204e-16. Thus, double precision
Fortran and Matlab have about the same precision.

\vspace{0.5cm}

Programming in double precision Fortran is straightforward and generally
requires only:

\vspace{0.5cm}

\begin{itemize}
\item  Declaring all floating (real) variables as double precision, \\
e.g., IMPLICIT DOUBLE PRECISION (A-H,O-Z)

\item  Converting all constants to double precision, e.g., $1.0D0$

\item  Declaring all in-line functions as double precision, e.g., DSIN (in
place of SIN)
\end{itemize}

\vspace{0.5cm}

Other numbers that characterize the arithmetic of a computer are:

\vspace{0.5cm}

\begin{itemize}
\item  The largest number the computer can accommodate (which defines the
\textit{dynamic range})

\item  The smallest number a computer can accommodate that is different than
zero
\end{itemize}

\vspace{0.5cm}

All of the key numbers are determined by the computer word length (number of
bits assigned to each word), and how many bits are assigned to the mantissa
and exponent. Typically, for a 32-bit word,

\vspace{0.5cm}

\begin{itemize}
\item  23 bits for mantissa

\item  1 bit for sign

\item  8 bits for exponent
\end{itemize}

\vspace{0.5cm}

In the case of Fortran, some typical numbers are

\vspace{0.5cm}

$
\begin{tabular}{cccc}
Precision & Bits per word & Machine epsilon & Dynamic range \\
Single & 32 & $10^{-7}$ & $10^{-37}-10^{37}$ \\
&  &  &  \\
Double & 64 & $10^{-14}$ & $10^{-300}-10^{300}$%
\end{tabular}
$

\vspace{0.5cm}

To conclude, the calculation of the machine epsilon is easy (consider the
preceding programs), and is probably the first thing we should do when using
a new computer/software.

\vspace{0.5cm}

We have considered briefly \textbf {computational roundoff}, due to the
finite wordlength e.g., what
happens in the last figure of the subtraction $2.999998-3.0$ will clearly
have a major effect on the answer (but we should not be working at this
threshold of precision).

\vspace{0.5cm}

Roundoff can also occur during storage of the number.  For example,
if two numbers with $n$ figures are multiplied, the full accuracy
product is of length $2n$.  However, when this product is stored,
only the first $n$ figures can be retained, and the least $n$ figures
are therefore dropped.

\vspace{0.5cm}

Finally, accumulation of roundoff
over many calculations can lead to computational failure unless we carry
enough significant figures (use a sufficiently long word length).

\vspace{0.5cm}

{\bf Truncation error}, as the name applies results from using 
a truncated series representation of a function.  For example, 
consider, the well-known infinite series expansion of the exponential
function around $x=0$

\[
e^{x}=1+x+x^{2}/2!+x^{3}/3!+\cdots 
\]

A series representation is required since the computer can only do
arithmetic, e.g., addition, subtraction, multiplication, division. If we
truncate this series (since the computer cannot sum an infinite series),

\[
e^{x}\approx 1+x+x^{2}/2!+x^{3}/3!
\]

Note that this truncated series is a third order polynomial. Generally,
polynomial approximations are used throughout scientific
computation, and they can be considered as truncated series approximations
of the functions of interest.

\vspace{0.5cm}

If we consider the Taylor series expansion of $e^{x}$ about the point $a$%

\[
e^{x}=e^{a}+e^{a}(x-a)+e^{a}(x-a)^{2}/2!+e^{a}(x-a)^{3}/3!+\cdots 
\]

or with $\Delta x=x-a$

\[
e^{x}=e^{a}(1+\Delta x+\Delta x^{2}/2!+\Delta x^{3}/3!+\cdots )
\]

Thus, we can consider the third order polynomial approximation as

\[
e^{x}\approx e^{a}(1+\Delta x+\Delta x^{2}/2!+\Delta x^{3}/3!+O(\Delta
x^{4}))
\]

where $O(\Delta x^{4})$ is interpreted as \textit{of order} $\Delta x^{4}$.
The order of the approximation is a common way of indicating its
truncation error.

\vspace{0.5cm}

Concerning the \textbf{condition}, we observed that the system

\[
\begin{array}{c}
1.49999999x_{1}+3.5x_{2}=2 \\
3.0x_{1}+7.0x_{2}=5
\end{array}
\]

is so ill-conditioned for a 32-bit computer, that the calculation of a
solution will fail. One way we might check for a numerically singular system
is to use the definition, i.e., effectively zero determinant.

\vspace{0.5cm}

However, the determinant is actually a poor practical measure of the
condition of a linear algebraic system because it scales poorly. What might
appear to be a small determinant (e.g., $10^{-40}$) may actually be for a
well-conditioned system. Thus, condition numbers have been developed that
are a better reflection of the actual condition of the linear algebraic
system.

\vspace{0.5cm}

For example, if the condition number, $cond$, is calculated, and if

\[
cond>1/eps
\]

then the system is numerically singular, and the computer with the given $eps
$ will not be able to solve the linear algebraic system.

\vspace{0.5cm}

The preceding discussion of condition number pertains to linear
algebraic systems.  More generally, we can consider the condition number
as the sensitivity to the problem system parameters.  If the condition
number is high, the system is very sensitive to the values of the
parameters, and therefore a solution to the model equations may
be precluded.  For example, extreme
sensitivity of differential equations to the initial conditions is
a condition for chaos.

\vspace{0.5cm}

To conclude the discussion of the first three topics

\vspace{0.5cm}

\begin{enumerate}
\item  Computer word length

\item  Roundoff and truncation

\item  Conditioning
\end{enumerate}

consider what is available in the NAG library. The following chapters
contain one or two line listings for each routine, with a total of more than
800 lines. Therefore, only the first listing or a small number of listings
in each chapter is indicated below:

\vspace{0.5cm}

Chapter F01 - Matrix Factorizations \\
F01ABF - Inverse of real symmetric positive-definite matrix using iterative
refinement (simplified parameter list) \\

Chapter F02 - Eigenvalue and Eigenvectors \\
F02BBF - Selected eigenvalues and eigenvectors of real symmetric matrix
(Black Box) \\

Chapter F03 - Determinants \\
F03AAF - Determinant of real matrix (Black Box) \\

Chapter F04 - Simultaneous Linear Equations \\
F04AAF - Solution of real simultaneous linear equations with multiple
right-hand sides (Black Box) \\

Chapter F05 - Orthogonalisation \\
F05AAF - Gram-Schmidt orthogonalisation of n vectors of order m \\

Chapter F06 - Linear Algebra Support Routines \\
F06BLF - Compute quotient of two real scalars, with overflow flag \\

Chapter F07 - Linear Equations (LAPACK) \\
F07ADF - LU factorization of real m by n matrix \\
F07AGF - Estimate condition number of real matrix, matrix already
factorized by F07ADF (SGECON/DGECON)\\

Chapter F08 - Least Squares and Eigenvalue Problems \\
F08AEF - QR factorization of real general rectangular matrix \\

Chapter F11 - Sparse Linear Algebra \\
F11BAF - Real sparse nonsymmetric linear systems, set-up for F11BBF \\

Chapter X02 - Machine Constants \\
X02AHF - Largest permissible argument for SIN and COS \\
X02AJF - Machine precision \\
X02AKF - Smallest positive model number \\
X02ALF - Largest positive model number

\vspace{0.5cm}

These chapters contain many routines pertaining to linear algebra, including
the calculation of the condition of the coefficient matrix, and basic
constants that determine the performance of the routines (e.g., the
precision reflected in the machine epsilon). A more complete discussion of
linear algebra will be given later.

\vspace{0.5cm}

Considering the final topic in this section, \textbf{Convergence and error
tolerances} the convergence of an iterative procedure, such as Newton's
method, is of course required to arrive at a solution with acceptable
accuracy. For example,

\[
x^{i+1}=x^{i}-\frac{f(x^{i})}{df(x^{i})/dx},i=0,1,\cdots
\]

is continued until

\[
\left| x^{j}-x^{j-1}\right| <\epsilon
\]

The question, is what is a reasonable value for $\epsilon$?

\vspace{0.5cm}

The choice of an error tolerance to determine when the iteration should be
terminated, is crucial, and is one of the most common causes for failures in
scientific computation.

\vspace{0.5cm}

Error tolerances are of two types: {\it absolute} and {\it relative}. Some
considerations in selecting these are illustrated with the following
examples:

\[
\begin{array}{cccc}
\text{Variable} & \text{Range} & \text{Type of error} & \text{Error tolerance%
} \\
\text{temperature (}^{o}K\text{)} & 300-400 & \text{abs} & 0.1 \\
\text{concentration} & 0-1 & \text{abs} & 0.001
\end{array}
\]

For temperature, an error tolerance of $0.1$ is reasonable, but it is
too large for concentration.  A tolerance of $0.001$ is appropriate
for concentration, but it is too stringent for temperature.  Thus,
using an absolute error tolerance for variables with widely different
magnitudes causes computational problems.

\vspace{0.5cm}

The solution to this problem is generally to use a relative tolerance.
For example,

\[
\begin{array}{cccc}
\text{Variable} & \text{Range} & \text{Type of error} & \text{Error tolerance%
} \\
\text{temperature (}^{o}K\text{)} & 300-400 & \text{rel} & 0.001 \\
\text{concentration} & 0-1 & \text{rel} & 0.001
\end{array}
\]

A tolerance of $0.001$ is appropriate for both variables.  However, a
relative error tolerance is not defined for a variable when it has a
zero value.

\vspace{0.5cm}

Thus, a tolerance that is a combination of the two is often used, e.g.,

\[
tol = relerr*y + abserr
\]

$tol$ is then compared with the estimated error from the numerical algorithm
to determine if the solution has achieved a reasonable accuracy.

\vspace{0.5cm}

As another example of the distinction between absolute and relative error,
in the least squares analysis of experimental data, if the data are
measured with constant absolute error, the error will be large for small
values of the measured variable.  If the data are measured with constant
relative error, the actual error is small for small values of the
measured data.  Thus, depending on the type of error, the reliability
of the least squares fit will be significantly different for small values
of the measured variable.

\vspace{0.5cm}

In summary, the careful selection of an error tolerance is often the
determining factor for successfully computing a numerical solution.
The thoughtless selection of a tolerance, which is tempting particularly
when using a library routine, generally leads to a failed calculation.

\newpage

\vspace{0.5cm}

\begin{center}
\textbf{\Huge Interpolation}
\end{center}

\begin{enumerate}
\item  Polynomials

\item  Splines
\end{enumerate}

\vspace{0.5cm}

Frequently in scientific and engineering analysis, we start with a
set of numbers (e.g., experimental data, output from a complicated
mathematical model), and we may then wish
to perform some mathematical analysis on the numbers (e.g., approximate
mathematical function, differentiation, integration).

\vspace{0.5cm}

But this implies that we have a mathematical function, when in fact,
we have only numbers.  Thus, is the mathematical analysis possible?

\vspace{0.5cm}

Yes, we can still do a mathematical analysis if we first develop a
functional representation or approximation of the numbers, that is,
develop a function that approximately fits the numbers in some sense.
In this discussion, we will consider standard and spline polynomials
(although there are many other possibilities).

\vspace{0.5cm}

Polynomials are selected because they are: (1) quite general and (2) can
be manipulated mathematically relatively easily.

\vspace{0.5cm}

Consider a set of three data pairs

\[
\begin{array}{cc}
x & f(x) \\
x_{1} & f(x_{1}) \\
x_{2} & f(x_{2}) \\
x_{3} & f(x_{3})
\end{array}
\]

We can consider a functional approximation of these of the form

\[
f(x)=(x-x_{2})(x-x_{3})a_{1}
\]
\[
+(x-x_{1})(x-x_{3})a_{2}
\]
\begin{equation}
+(x-x_{1})(x-x_{2})a_{3}
\tag{1}
\end{equation}

where $a_{1},a_{2}$ and $a_{3}$ are constants to be determined. If we
substitute $x=x_{1}$ in eq. (1),

\[
f(x_{1})=(x_{1}-x_{2})(x_{1}-x_{3})a_{1}+0a_{2}+0a_{3}
\]

we can solve for $a_{1}$,

\[
a_{1}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}
\]

Similarly, substituting $x=x_{2}$ in eq. (1) gives

\[
a_{2}=\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}
\]

and finally, substituting $x=x_{3}$ in eq. (1) gives

\[
a_{3}=\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}
\]

Substituting these coefficients in eq. (1), we arrive at the {\it
second order Lagrange interpolation polynomial}

\[
f(x)=\frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})}f(x_{1})
\]
\[
+\frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}f(x_{2})
\]
\begin{equation}
+\frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}f(x_{3})  \tag{2}
\end{equation}

Note that eq. (2) has three important properties:

\vspace{0.5cm}

\begin{itemize}
\item  It passes through the points $
(x_{1},f(x_{1}))$,$(x_{2},f(x_{2}))$,$(x_{3},f(x_{3}))$ exactly, i.e., it
returns the original data exactly.

\item  It can be applied to unequally spaced data, i.e., for $%
(x_{1}-x_{2})\neq (x_{1}-x_{3})\neq (x_{2}-x_{3})$.

\item  Since it is a polynomial, it can be easily integrated,
differentiatiated, and generally manipulated relatively easily
mathematically.
\end{itemize}

\vspace{0.5cm}

However, there is no guarantee that it is well behaved between the data
points, and of course, it is limited to just three data pairs.

\vspace{0.5cm}

If we have more than three data pairs, we can still use a polynomial.
For example, if we start with the model

\[
y(x)=a_{0}+a_{1}x+a_{2}x^{2}
\]

then a sum of squared errors for $n$ data pairs gives

\[
E(a_{0},a_{1}a_{2})=\sum_{i=1}^{n}[y_{i}-y(x_{i})]^{2}=%
\sum_{i=1}^{n}[y_{i}-(a_{0}+a_{1}x_{i}+a_{2}x_{i}^{2})]^{2}
\]

To minimize $E$

\[
\frac{\partial E}{\partial a_{0}}=2%
\sum_{i=1}^{n}[y_{i}-(a_{0}+a_{1}x_{i}+a_{2}x_{i}^{2})](-1)=0
\]

\[
\frac{\partial E}{\partial a_{1}}=2%
\sum_{i=1}^{n}[y_{i}-(a_{0}+a_{1}x_{i}+a_{2}x_{i}^{2})](-x_{i})=0
\]
\[
\frac{\partial E}{\partial a_{2}}=2%
\sum_{i=1}^{n}[y_{i}-(a_{0}+a_{1}x_{i}+a_{2}x_{i}^{2})](-x_{i}^{2})=0
\]

or

\[
\begin{array}{c}
a_{0}\sum_{i=1}^{n}1+a_{1}\sum_{i=1}^{n}x_{i}+a_{2}\sum_{i=1}^{n}x_{i}^{2}=%
\sum_{i=1}^{n}y_{i} \\
a_{0}\sum_{i=1}^{n}x_{i}+a_{1}\sum_{i=1}^{n}x_{i}^{2}+a_{2}%
\sum_{i=1}^{n}x_{i}^{3}=\sum_{i=1}^{n}x_{i}y_{i} \\
a_{0}\sum_{i=1}^{n}x_{i}^{2}+a_{1}\sum_{i=1}^{n}x_{i}^{3}+a_{2}%
\sum_{i=1}^{n}x_{i}^{4}=\sum_{i=1}^{n}x_{i}^{2}y_{i}
\end{array}
\]

which is a $3x3$ linear system for the three coefficients $a_{0},a_{1},
a_{2}$.  The number of data pairs, $n$, is arbitrary, and the resulting
third order polynomial will pass through the data 11smoothly'' (in the
least squares sense).

\vspace{0.5cm}

Other models are possible, for example:

\[
y(x)=\sum_{j=0}^{m}a_j x^j\text{ (polynomial of degree $m$)}
\]

\[
y(x)=a_{1}x+a_{2}x^{2}\text{ (linear, includes the origin)}
\]

\[
y(x)=a_{0}+a_{1}x+a_{2}e^{b_{0}x}\text{ (linear, $b_0$ fixed)}
\]

\[
y(x)=a_{0}+a_{1}x+a_{2}e^{b_{0}x}\text{ (nonlinear, $b_0$ variable)}
\]

For the last model, we would have to solve nonlinear equations (to be
considered subsequently).

\vspace{0.5cm}

{\bf Splines}

\vspace{0.5cm}

Splines are piecewise polynomials, selected to ensure continuity not only in
the dependent variable where the pieces connect, but also in their
derivatives. For example, when using cubic splines of the form

\[
f(x)=a_{0}+a_{1}(x-x_{i})+a_{2}(x-x_{i})^{2}+a_{3}(x-x_{i})^{3}
\]

we require that two splines coming together at point $x_{i}$ have continuity
in $f(x),df(x)/dx$ and $d^{2}f(x)/dx^{2}$. These three continuity
requirements are applied at each of the points, $x_{1},x_{2},\cdots ,
x_{N}$ to evaluate the coefficients $a_{1},a_{2}$ and $a_{3}$
at each of the points (also, we have $a_{0}=f(x_{i})$ to ensure the splines
go through the given points). This procedure leads to a system of
tridiagonal equations that is solved for $a_{1},a_{2},a_{3}$ at each
of the points, $x_{1},x_{2},\cdots ,x_{N\text{ }}.$

\vspace{0.5cm}

The solution of the tridiagonal system can be done very efficiently because
the coefficient matrix is:

\begin{itemize}

\item  Tridiagonal (so storage is only $3xn$)

\item  Symmetric

\item  Nonsingular

\item  Well conditioned

\item  Diagonally dominant (pivoting not required with
Gaussian elimination)

\end{itemize}

Thus, updating the spline as the calculation proceeds can be
done efficiently, e.g., in the solution of PDEs.

\vspace{0.5cm}

NAG includes the following routines for functional approximation:

\vspace{0.5cm}

Chapter E01 - Interpolation\\
\\        E01AAF - Interpolated values, Aitken's technique, unequally
spaced data, one variable\\
\\       E01ABF - Interpolated values, Everett's formula, equally
spaced data, one variable\\
\\       E01AEF - Interpolating functions, polynomial interpolant,
data may include , derivative values, one variable\\
\\       E01BAF - Interpolating functions, cubic spline interpolant,
one variable\\
\\       E01BEF - Interpolating functions, monotonicity-preserving,
piecewise cubic  Hermite, one variable\\
\\       E01BFF - Interpolated values, interpolant computed by E01BEF,
function only, one variable\\
\\       E01BGF - Interpolated values, interpolant computed by E01BEF,
function and 1st derivative, one variable\\
\\       E01BHF - Interpolated values, interpolant computed by E01BEF,
definite integral, one variable\\
\\       E01DAF - Interpolating functions, fitting bicubic spline, data
on rectangular grid\\
\\       E01RAF - Interpolating functions, rational interpolant, one
variable\\
\\       E01RBF - Interpolated values, evaluate rational interpolant
computed by E01RAF, one variable\\
\\       E01SAF - Interpolating functions, method of Renka and Cline,
two variables\\
\\       E01SBF - Interpolated values, evaluate interpolant computed by
E01SAF, two variables\\
\\       E01SEF - Interpolating functions, modified Shepard's method,
two variables\\
\\       E01SFF - Interpolated values, evaluate interpolant computed
by E01SEF, two variables\\
\\       E01SGF - Interpolating functions, modified Shepard's method,
two variables\\
\\       E01SHF - Interpolated values, evaluate interpolant computed
by E01SGF, functiona nd first derivatives, two variables\\
\\       E01TGF - Interpolating functions, modified Shepard's method,
three variables\\
\\       E01THF - Interpolated values, evaluate interpolant computed
by E01TGF, function and first derivatives, three variables\\

\vspace{0.5cm}

Chapter E02 - Curve and Surface Fitting
\\       E02ACF - Minimax curve fit by polynomials\\
\\       E02ADF - Least-squares curve fit, by polynomials, arbitrary
data points\\
\\       E02AEF - Evaluation of fitted polynomial in one variable
from Chebyshev series form (simplified parameter list)\\
\\       E02AFF - Least-squares polynomial fit, special data points
(including interpolation)\\
\\       E02AGF - Least-squares polynomial fit, values and derivatives
may be constrained, arbitrary data points,\\
\\       E02AHF - Derivative of fitted polynomial in Chebyshev series
form\\
\\       E02AJF - Integral of fitted polynomial in Chebyshev series form\\
\\       E02AKF - Evaluation of fitted polynomial in one variable, from
Chebyshev series form\\
\\       E02BAF - Least-squares curve cubic spline fit (including
interpolation)\\
\\       E02BBF - Evaluation of fitted cubic spline, function only\\
\\       E02BCF - Evaluation of fitted cubic spline, function and
derivatives\\
\\       E02BDF - Evaluation of fitted cubic spline, definite integral\\
\\       E02BEF - Least-squares cubic spline curve fit, automatic knot
placement\\
\\       E02CAF - Least-squares surface fit by polynomials, data on lines\\
\\       E02CBF - Evaluation of fitted polynomial in two variables\\
\\       E02DAF - Least-squares surface fit, bicubic splines\\
\\       E02DCF - Least-squares surface fit by bicubic splines with
automatic knot placement, data on rectangular grid\\
\\       E02DDF - Least-squares surface fit by bicubic splines with
automatic knot placement, scattered data\\
\\       E02DEF - Evaluation of a fitted bicubic spline at a vector of
points\\
\\       E02DFF - Evaluation of a fitted bicubic spline at a mesh of points\\
\\       E02GAF - L(1)-approximation by general linear function\\
\\       E02GBF - L(1)-approximation by general linear function subject
to linear inequality constraints\\
\\       E02GCF - L(infinity)-approximation by general linear function\\
\\       E02RAF - Pade-approximants\\
\\       E02RBF - Evaluation of fitted rational function as computed by
E02RAF\\
\\       E02ZAF - Sort 2-D data into panels for fitting bicubic splines

\vspace{0.5cm}

Chapter E04 - Minimizing or Maximizing a Function\\
(routines not listed)

\newpage

\vspace{0.5cm}

\begin{center}
{\bf \Huge Integration}
\end{center}

\begin{enumerate}

\item  Polynomials

\item  Gaussian quadrature

\item  Splines

\item  Adaptive
\end{enumerate}

\vspace{0.5cm}

Returning to the second order Lagrange interpolation polynomial
integration of $f\left( x\right) $,

\[
I=\int_{x_{1}}^{x_{3}}f(x)dz
\]

\[
=\int_{x_{1}}^{x_{3}}\frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})}%
f(x_{1})
\]
\[
+\frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}f(x_{2})
\]
\[
+\frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}f(x_{3})dx
\]

Integrating the second order polynomials, we arrive at a quadrature rule
(numerical integration rule) for nonuniformly spaced data (which can be
applied over successive sets of three points or two intervals).

\vspace{0.5cm}

For the special case of a uniform grid, $(x_2 - x_1) = (x_3 - x_2) =
\Delta x$, $I$ becomes

\[
I=\int_{x_{1}}^{x_{3}}f(x)dx=\frac{\Delta x}{3}(f(x_{1})+4f(x_{2})+f(x_{3}))
\]

which is {\bf Simpson's rule.}

\vspace{0.5cm}

Some important details about Simpson's rule:

\begin{itemize}
\item  An even number of ``panels'' (odd number of grid points) is required.
In other words, Simpson's rule can be extended to $n$ points as

\[
I=\int_{x_{1}}^{x_{n}}f(x)dx=\frac{\Delta x}{3}(f(x_{1})+4f(x_{2})+f(x_{3}))
\]
\[
+\frac{\Delta x}{3}(f(x_{3})+4f(x_{4})+f(x_{5}))
\]
\[
+\cdots +\frac{\Delta x}{3}(f(x_{n-2})+4f(x_{n-1})+f(x_{n}))
\]

\[
=\frac{\Delta x}{3}(f(x_{1})+4f(x_{2})+2f(x_{3})+4f(x_{4})+\cdots
\]
\begin{equation}
+2f(x_{n-2})+4f(x_{n-1})+f(x_{n}))  \tag{1}
\end{equation}

\item  Eq. (1) is surprisingly accurate. In fact, it is {\bf exact for
third order polynomials}, even thought it is derived using a second order
polynomial approximation for the integrand (in general, for equally
spaced points, quadrature formulas based on a polynomial of order $2n$ are
exact for polynomials of order $2n+1$.

\item  Although eq. (1) is generally sufficiently accurate (this, of course,
depends on $\Delta x$) it is based on equally spaced data (constant $\Delta x
$).

\end{itemize}

\vspace{0.5cm}

Even more accurate integrals can be computed using unequally
spaced data or points; the best example of this idea is
{\bf Gaussian quadrature}:

\begin{itemize}

\item  The spacing of the points is selected to give an exact intergal for
a polynomial of maximum order.  Specifically, integrating a polynomial
of order $n-1$ gives
an integral that is exact for polynomials of order $2n-1$.  For example,
with $n=2$ $n-1 = 2 - 1 = 1$ is a linear approximation of the integrand,
but the computed integral is exact for a polynomial of order $2n-1 = 4-1
=3$!

\item The spacing
of the points generally corresponds to the roots of an orthogonal
polynomial.  The location of the points is determined by the number
of points and the limits of integration, i.e., $a\leq x\leq b,
0\leq x\leq \infty , -\infty \leq x\leq \infty $.

\end{itemize}

\vspace{0.5cm}

Gaussian quadrature, as other quadrature rules (rectangular, trapezoidal,
Simpson's rules) is based on a weighted sum

\vspace{0.5cm}

\[
\int_{a}^{b}f(x)dx\approx \sum_{i=1}^{n}w_{i}f(x_{i})
\]

Now, however, we are going to choose the points $x_{i}$ as well as the
weights $w_{i}$ to get a numerical integral that is exact for a
polynomial of as high degree as possible (degree $2n-1$ when $f(x)$ is
a polynomial of order $n-1$ - the basic idea of Gaussian quadrature). A
short table of points and weights follows (Borse, 1997):

\vspace{0.5cm}

\[
\begin{array}{cccc}
n & i & x_{i} & w_{i} \\
2 & 1 & 0.5773502692 & 1.0 \\
3 & 0 & 0.0 & 0.8888888889 \\
  & 1 & 0.7745966692 & 0.5555555556\\
5 & 0 & 0.0 & 0.5688888889\\
  & 1 & 0.5384693101 & 0.4786286705\\
  & 2 & 0.9061798459 & 0.2369268850
\end{array}
\]

\vspace{0.5cm}

The points are in the interval  $-1\leq x_{i}\leq 1$ (they are symmetric
around $x=0$). Thus, the limits of the integral must be $a=-1,$$b=1$).
This can easily be achieved by a linear change of variables

\[
I=\int_{a}^{b}f(x)dx=\int_{-1}^{1}f(\lambda )d\lambda
\]

where

\[
x=\frac{b-a}{2}\lambda +\frac{b+a}{2}\lambda
\]

As an example,

\[
I=\int_{0}^{2}\sin xdx
\]

Then, $x=\lambda +1$, and the integral in $\lambda $ is

\[
I=\int_{-1}^{1}\sin (\lambda +1)d\lambda
\]

Then, by Gaussian quadrature for $n=5,$

\[
I\approx w_{0}f(0)+w_{1}\left[ f(x_{1})+f(-x_{1})\right] +w_{2}\left[
f(x_{2})+f(-x_{2})\right]
\]

\[
=0.5688888889\sin (1.0)
\]
\[
+0.4786286705\left[ \sin (0.5384693101+1)\right] +\sin (-0.5384693101+1)
\]
\[
+0.2369268850\left[ \sin (0.9061798459+1)\right] +\sin (-0.9061798459+1)
\]

\[
=1.4161467\text{; exact answer }=1-\cos (2)=1.4161468\ldots
\]

The method can also be applied to integrals with infinite limits by
considering a product integrand

\[
\int_{a}^{b}g(x)f(x)dx\approx \sum_{i=1}^{n}w_{i}f(x_{i})
\]

Then, the methods are summarized as

\[
\begin{array}{cccc}
\text{Integral type} & g(x) & \left[ a,b\right]  & \text{Method name} \\
\int_{-1}^{1}f(x)dx & 1 & \left( -1,1\right)  & \text{Gauss-Legendre} \\
\int_{-1}^{1}\frac{f(x)}{\sqrt{1-x^{2}}}dx & (1-x^{2})^{-1/2} & (-1,1) &
\text{Gauss-Chebyshev} \\
\int_{0}^{\infty }x^{q}e^{-x}f(x)dx & x^{q}e^{-x} & (0,\infty ) & \text{%
Gauss-Laguerre} \\
\int_{-\infty }^{\infty }e^{-x^{2}}f(x)dx & e^{-x^{2}} & (-\infty ,\infty )
& \text{Gauss-Hermite}
\end{array}
\]

The choice of a method is determined primarily by the limits of the
integral. The name indicates the orthogonal polynomial for which the
roots determine the points (we previously used Gauss-Legendre, and the
tabulated values of $x_{i}$ are the roots of the Legendre polynomial
of order $n$).

\vspace{0.5cm}

Borse, G. J. (1997), \textit{Numerical Methods with Matlab; A Resource for
Scientists and Engineers}, PWS Publishing Company, Boston

\vspace{0.5cm}

We can also use the spline as an approximation for an integrand, that is
\textbf{spline quadrature}

\[
I=\int_{x_{0}}^{x_{n}}f(x)
\]

with

\[
f(x)=a_{0}+a_{1}(x-x_{i})+a_{2}(x-x_{i})^{2}+a_{3}(x-x_{i})^{3}
\]

The integration is easily carried out, and the spline coefficients can be
evaluated by the solution of tridiagonal equations, as described previously.

\vspace{0.5cm}

Another approach to improving accuracy is to have the code sense the
rate of change of $f(x)$, then insert points automatically in regions of
high rate of change. This is termed {\bf adaptive quadrature}.
Consider for example the algorithm

\vspace{0.5cm}

\[
\begin{array}{cc}
1. & \text{Select }h,I_{0}=0,i=0 \\
2. & \text{If }x_{i}>x_{n}\text{stop} \\
3. & \text{ }I_{i+1}^{p}=I_{i}+f(x_{i})h, \\
4. & \varepsilon =\left[ f(x_{i+1})-f(x_{i})\right] h/2 \\
5. & \text{If }\varepsilon >tol,h=h/2,\text{ go to 3} \\
6. & \text{If }\varepsilon \leq tol,I_{i+1}^{c}=I_{i+1}^{p}+\varepsilon ,%
\text{ }i=i+1,x_{i}=x_{i}+h\text{ go to 2}
\end{array}
\]

This is an adaptive trapezoidal rule.

\vspace{0.5cm}

Adaptive quadrature is readily available in a series of library routines.
For example, QUANC8 from Forsythe, Malcolm and Moler (1977) is named from
Quadrature, Adaptive, Newton-Cotes 8-panel. A typical call to QUANC8 is

\begin{verbatim}

CALL QUANC8(FUN,A,B,ABSERR,RELERR,RESULT,ERREST,NOFUN,FLAG)

\end{verbatim}

where

\vspace{0.5cm}

\begin{center}
\begin{tabular}{ll}
FUN      & External to compute the integrand (input)\\
A        & Lower limit of integration (input)\\
B        & Upper limit of integration (input)\\
ABSERR   & Absolute error tolerance (input)\\
RELERR   & Relative error tolerance (input)\\
RESULT   & Computed integral (output)\\
ERREST   & Estimated error for RESULT (output)\\
\end{tabular}
\end{center}

\newpage

\begin{center}
\begin{tabular}{ll}
NOFUN    & Number of integrand evaluations (output)\\
FLAG     & Error flag indicating where calculation failed (output)\\
\end{tabular}
\end{center}

\vspace{0.5cm}

Matlab has an implementation of QUANC8, as well as an adaptive
Simpson's rule routine:

\vspace{0.5cm}

\begin{itemize}

\item quad uses an adaptive, recursive Simpson's rule; the call is
\begin{verbatim}
q = quad('function',a,b,tol)
\end{verbatim}

\item quad8 implements an adaptive recursive Newton-Cotes 8 panel rule
(as in QUANC8); the call is
\begin{verbatim}
q=quad8('function',a,b,tol)
\end{verbatim}
\end{itemize}

\vspace{0.5cm}

The routines in NAG for quadrature are:

\vspace{0.5cm}

Chapter D01 - Quadrature
\\       D01AHF - 1-D quadrature, adaptive, finite interval, strategy
due to Patterson, suitable for well-behaved integrands\\
\\       D01AJF - 1-D quadrature, adaptive, finite interval, strategy
due to Piessens and deDoncker, allowing for badly-behaved integrands\\
\\       D01AKF - 1-D quadrature, adaptive, finite interval, method
suitable for oscillating functions\\
\\       D01ALF - 1-D quadrature, adaptive, finite interval, allowing
for singularities at user-specified break-points\\
\\       D01AMF - 1-D quadrature, adaptive, infinite or semi-infinite
interval\\
\\       D01ANF - 1-D quadrature, adaptive, finite interval, weight
function cos(wx) or sin(wx)\\
\\       D01APF - 1-D quadrature, adaptive, finite interval, weight
function with end-point singularities of algebraico-logarithmic type\\
\\       D01AQF - 1-D quadrature, adaptive, finite interval, weight
function 1/(x-c), Cauchy principal value (Hilbert transform)\\
\\       D01ARF - 1-D quadrature, non-adaptive, finite interval with
provision for indefinite integrals\\
\\       D01ASF - 1-D quadrature, adaptive, semi-infinite interval,
weight function cos(wx) or sin(wx)\\
\\       D01ATF - 1-D quadrature, adaptive, finite interval, variant of
D01AJF efficient on vector machines\\
\\       D01AUF - 1-D quadrature, adaptive, finite interval, variant of
D01AKF efficient on vector machines\\
\\       D01BAF - 1-D Gaussian quadrature\\
\\       D01BBF - Pre-computed weights and abscissae for Gaussian
quadrature rules, restricted choice of rule\\
\\       D01BCF - Calculation of weights and abscissae for Gaussian
quadrature rules, general choice of rule\\
\\       D01BDF - 1-D quadrature, non-adaptive, finite interval\\
\\       D01DAF - 2-D quadrature, finite region\\
\\       D01EAF - Multi-dimensional adaptive quadrature over
hyper-rectangle, multiple integrands\\
\\       D01FBF - Multi-dimensional Gaussian quadrature over
hyper-rectangle\\
\\       D01FCF - Multi-dimensional adaptive quadrature over
hyper-rectangle\\
\\       D01FDF - Multi-dimensional quadrature, Sag-Szekeres method,
general product region or n-sphere\\
\\       D01GAF - 1-D quadrature, integration of function defined
by data values, Gill-Miller method\\
\\       D01GBF - Multi-dimensional quadrature over hyper-rectangle,
Monte Carlo method\\
\\       D01GCF - Multi-dimensional quadrature, general product region,
number-theoretic method\\
\\       D01GDF - Multi-dimensional quadrature, general product region,
number-theoretic method, variant of D01GCF efficient on vector machines\\
\\       D01GYF - Korobov optimal coefficients for use in D01GCF or D01GDF,
when number of points is prime\\
\\       D01GZF - Korobov optimal coefficients for use in D01GCF or D01GDF,
when number of points is product of two primes\\
\\       D01JAF - Multi-dimensional quadrature over an n-sphere, allowing
for badly-behaved integrands\\
\\       D01PAF - Multi-dimensional quadrature over an n-simplex

\newpage

\vspace{0.5cm}

\begin{center}
{\bf \Huge Differentiation}
\end{center}

\begin{enumerate}
\item  Finite differences

\item  Symbolic computation

\item Automatic differentiation (to be discussed later)
\end{enumerate}

Returning to the Lagrange second order interpolation polynomial

\[
f(x)=\frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})}f(x_{1})
\]
\[
+\frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}f(x_{2})
\]
\begin{equation}
+\frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}f(x_{3})  \tag{2}
\end{equation}

we can differentiate eq. (2) to obtain an approximation to the
derivative of $f(x)$

\[
\frac{df(x)}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}\left\{
(x-x_{2})+(x-x_{3})\right\}
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}\left\{
(x-x_{1})+(x-x_{3})\right\}
\]

\begin{equation}
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}\left\{
(x-x_{1})+(x-x_{2})\right\}  \tag{3}
\end{equation}

Eq. (3) could then be used to calculate the derivative at any $x_{1}\leq
x\leq x_{3}$. In particular,

\vspace{0.5cm}

\[
\frac{df(x_{1})}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}\left\{
(x_{1}-x_{2})+(x_{1}-x_{3})\right\}
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}(x_{1}-x_{3})
\]
\begin{equation}
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}(x_{1}-x_{2})  \tag{4a}
\end{equation}

\vspace{0.5cm}

\[
\frac{df(x_{2})}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}%
(x_{2}-x_{3})
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}\left\{
(x_{2}-x_{1})+(x_{2}-x_{3})\right\}
\]

\begin{equation}
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}(x_{2}-x_{1})  \tag{4b}
\end{equation}

\vspace{0.5cm}

\[
\frac{df(x_{3})}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}%
(x_{3}-x_{2})
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}(x_{3}-x_{1})
\]

\begin{equation}
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}\left\{
(x_{3}-x_{1})+(x_{3}-x_{2})\right\}  \tag{4c}
\end{equation}

Note that the RHSs of eqs. (4a) to (4c) are just constants, and therefore
eqs. (4) can be considered as a \textbf{three point differentiation matrix
on a nonuniform grid}.

\vspace{0.5cm}

If we consider the special case of a uniform grid, $\Delta
x=x_{2}-x_{1}=x_{3}-x_{2}$; $2\Delta x = x_{3}-x_{1}$, eqs. (4) become

\[
\frac{df(x_{1})}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}\left\{
(x_{1}-x_{2})+(x_{1}-x_{3})\right\}
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}(x_{1}-x_{3})
\]
\[
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}(x_{1}-x_{2})
\]

\vspace{0.5cm}

\[
=\frac{f(x_{1})}{(-\Delta x)(-2\Delta x)}\left\{ -\Delta x-2\Delta x\right\}
\]
\[
+\frac{f(x_{2})}{(\Delta x)(-\Delta x)}(-2\Delta x)
\]
\[
+\frac{f(x_{3})}{(2\Delta x)(\Delta x)}(-\Delta x)
\]

\vspace{0.5cm}

\[
=\frac{-3}{2\Delta x}f(x_{1})
\]
\[
+\frac{4}{2\Delta x}f(x_{2})
\]
\[
+\frac{-1}{2\Delta x}f\left( x_{3}\right)
\]

\vspace{0.5cm}

\[
=\frac{1}{2!\Delta x}(
\begin{array}{ccc}
-3 & 4 & -1
\end{array}
)\mathbf{f}
\]

\vspace{0.5cm}

\[
\frac{df(x_{2})}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}%
(x_{2}-x_{3})
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}\left\{
(x_{2}-x_{1})+(x_{2}-x_{3})\right\}
\]
\[
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}(x_{2}-x_{1})
\]

\vspace{0.5cm}

\[
=\frac{f(x_{1})}{(-\Delta x)(-2\Delta x)}(-\Delta x)
\]
\[
+\frac{f(x_{2})}{(\Delta x)(-\Delta x)}(\Delta x-\Delta x)
\]
\[
+\frac{f(x_{3})}{(2\Delta x)(\Delta x)}(\Delta x)
\]

\vspace{0.5cm}

\[
=-\frac{1}{2\Delta x}f(x_{1})+0f(x_{2})+\frac{1}{2\Delta x}f(x_{3})
\]

\vspace{0.5cm}

\[
=\frac{1}{2!\Delta x}(
\begin{array}{ccc}
-1 & 0 & 1
\end{array}
)\mathbf{f}
\]

\vspace{0.5cm}

\[
\frac{df(x_{3})}{dx}=\frac{f(x_{1})}{(x_{1}-x_{2})(x_{1}-x_{3})}%
(x_{3}-x_{2})
\]
\[
+\frac{f(x_{2})}{(x_{2}-x_{1})(x_{2}-x_{3})}(x_{3}-x_{1})
\]
\[
+\frac{f(x_{3})}{(x_{3}-x_{1})(x_{3}-x_{2})}\left\{
(x_{3}-x_{1})+(x_{3}-x_{2})\right\}
\]

\vspace{0.5cm}

\[
=\frac{f(x_{1})}{(-\Delta x)(-2\Delta x)}(\Delta x)
\]
\[
+\frac{f(x_{2})}{(\Delta x)(-\Delta x)}(2\Delta x)
\]
\[
+\frac{f(x_{3})}{(2\Delta x)(\Delta x)}(2\Delta x+\Delta x)
\]

\vspace{0.5cm}

\[
=\frac{1}{2\Delta x}f(x_{1})-\frac{4}{2\Delta x}f(x_{2})+\frac{3}{2\Delta x}%
f(x_{3})
\]

\vspace{0.5cm}

\[
=\frac{1}{2!\Delta x}(
\begin{array}{ccc}
1 & -4 & 3
\end{array}
)\mathbf{f}
\]

Combining results, we arrive at the three point, uniformly spaced
differentiation matrix

\[
\frac{d\mathbf{f}}{dx}=\frac{1}{2!\Delta x}\left[
\begin{array}{ccc}
-3 & 4 & -1 \\
-1 & 0 & 1 \\
1 & -4 & 3
\end{array}
\right] \mathbf{f+}O(\Delta x^{2})
\]

This matrix will be used later in the numerical integration of PDEs.

\vspace{0.5cm}

A Fortran subroutine, weights.for, by Bengt Fornberg, is available for
the convenient calculation of the weights in differentiation formulas.
For example, the following fourth order, differentiation matrix is easily
obtained using weights.for

\[
\frac{d\mathbf{u}}{dx}=
\]
\[
\frac{1}{4!\Delta x}\left[
\begin{array}{cccccc}
-50 & 96 & -72 & 32 & -6 & i=1 \\
-6 & -20 & 36 & -12 & 2 & i=2 \\
2 & -16 & 0 & 16 & -2 & i=3,\cdots ,N-2 \\
-2 & 12 & -36 & 20 & 6 & i=N-1 \\
6 & -32 & 72 & -96 & 50 & i=N
\end{array}
\right] \mathbf{u}
\]
\[
+O(\Delta x^{4})
\]

Note that in these two differentiation matrices, the weighting
coefficients in any give row sum to zero, which is required
to correctrly differentiate a constant (to zero).

\vspace{0.5cm}

Also, the second order matrix differentiates a second order polynomial
exactly, and the fourth order matrix differentiates a fourth order
matrix exactly.  For example, using the second row of the second
order matrix,

\[
\frac{df(x)}{dx}=\frac{d\left( a_{0}+a_{1}x+a_{2}x^{2}\right) }{dx}
\]
\[
=\frac{-1[a_{0}+a_{1}(x-\Delta x)+a_{2}(x-\Delta
x)^{2}]}{2\Delta x}
\]
\[
+\frac{
1[a_{0}+a_{1}(x+\Delta x)+a_{2}(x+\Delta x)^{2}]}{2\Delta x}
\]

\[
=\frac{2a_{1}\Delta x+4a_{2}x\Delta x}{2\Delta x}=a_{1}+2a_{2}x
\]

as expected.

\vspace{0.5cm}

We will later use the formulas defined by these differentiation matrices
for first order derivatives, as well as for second derivatives, in
the numerical integration of PDEs.

\vspace{0.5cm}

Subroutine weights.for can be used to compute the finite difference
(FD) weights for:

\begin{itemize}

\item  Derivatives of any order

\item  Approximations of any order

\item  Any grid point

\item  Equal and unequally spaced grids

\end{itemize}

\vspace{0.5cm}

NAG has the following routine for the calculation of derivatives:

\vspace{0.5cm}

D04AAF - Numerical differentiation, derivatives up to order 14,
function of one real variable

\vspace{0.5cm}

{\bf Symbolic computation} (analytical differentiation) is also possible
using any of a number of computer algebra systems (CAS) such as Maple
and Mathematica.  Symbolic (exact) mathematics of other types is also
possible with these systems, partial differentiation, integration,
linear algebra, etc.

\vspace{0.5cm}

As an example of symbolic differentiation with Maple, consider (Nicolaides
and Walkington (1996))

\begin{verbatim}

>u:=x*tan(x*y);

                     u:=x tan(x y)

>diff(u,x);
                                      2
            tan(x y) + x (1 + tan(x y) ) y

>diff(u,y):
                   2               2
                  x  (1 + tan( x y) )

>diff(u,x,y);
                   2         2                       
    2 (1 + tan(x y) ) x + 2 x  tan(x y)

                 2
    (1 + tan(x y) ) y

>diff(x^2*y^3,y$3);
                         2
                      6 x
\end{verbatim}

Briefly, here are some examples of integration:

\begin{verbatim}

Indefinite integration:

>int(sin(ln(x)),x);

           1/2 x sin(ln(x)) - cos(ln(x)))

Definite integration:

>int(x*ln(sin(x)),x=0..Pi);

                         2
                 - 1/2 Pi   ln(2)

Numerical integration:

>evalf(int(x*ln(sin(x)),x=0..Pi);

                  -3.420544232

>evalf(-Pi^2*ln(2)/2);

                  -3.420544233

\end{verbatim}

The integration is done by adaptive quadrature for 10-figure
accuracy (the default for all Maple calculations).  This can
lead to some very long computation times.  If the accuracy is
reduced to 4-5 figures, the integration can be much faster.

\vspace{0.5cm}

In other words, Maple has variable precision arithmetic, so
in principle, calculations of arbitrary precision can be
performed (at the price of greater computation time).  Thus,
Maple has, in effect, a variable machine epsilon.

\vspace{0.5cm}

Additionally, Maple and NAG have recently been combined.  Matlab
has a ``NAG foundation toolbox'' with about 50 basic NAG routines.

\vspace{0.5cm}

Nicolaides, R., and N. Walkington (1996), {\it Maple: A Comprehensive
Introduction}, Cambridge University Press, Cambridge

\newpage

\begin{center}
{\Huge \textbf{References and Software}}
\end{center}

\vspace{0.5cm}

\begin{enumerate}

\item  Forsythe, et al (1977)

\item  Kahaner, et al (1989)

\item  Golub, et al (1992)

\item  Press, et al (1992)

\item  NAG Library, Mark 18 Release

\end{enumerate}

\vspace{0.5cm}

\end{document}
