\chapter{Numerical Techniques}

The dynamic systems we analyze in this book are modeled by systems of ordinary differential equations. A main use of a mathematical model is to predict the system transient behavior. In general, finding analytical, exact, solution to the modeling equations is impossible. Thus, we must resort to approximate solutions using numerical techniques. We begin our presentation of numerical techniques that allow us to determine approximate solution to the given system of differential equations with a method based on a Taylor series expansion.

\section{The Method of Taylor Series}
The basis for many numerical methods for solving differential equations is the Taylor formula. Suppose we are given a model of the form
\[
\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t)), \quad \mb{\bm$x$}(t_0)=\mb{\bm$x$}_0.
\]
The method we discuss in this section depends on the fact that if $\mb{\bm$f$}(t,\mb{\bm$x$},\mb{\bm$u$})$ is analytic in its $n+m+1$ variables $x_1,\ldots, x_n,u_1,\ldots ,u_m, t$ in some $(n+m+1)$-dimensional region ${\cal D}$, and if $\mb{\bm$x$}(t)$ is a solution of $\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$ u$}(t)), \; \mb{\bm$x$}(t_0)=\mb{\bm$x$}_0$, where $(t_0,\mb{\bm$x$}_0)$ lies in ${\cal D}$, then $\mb{\bm$ x$}(t)$ is analytic. We expand $\mb{\bm$x$}(t)$ into a Taylor series about the initial point $\mb{\bm$x$}_0$ to obtain
\begin{eqnarray*}
\mb{\bm$x$}(t)&=&\mb{\bm$x$}_0+\dot{\mb{\bm$ x$}}(t_0)(t-t_0)+\ddot{\mb{\bm$ x$}}(t_0)\frac{(t-t_0)^2}{2}+\cdots \\
&=&\mb{\bm$x$}_0+\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$ u$}(t_0))(t-t_0)\\
{}&{}&+\left[\frac{\partial \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0))}{\partial t}+\sum_{i=1}^{n} \frac{\partial \mb{\bm$f$}(t_0,\mb{\bm$ x$}_0,\mb{\bm$u$}(t_0))}{\partial x_i} \dot{x}_i \right.\\
{}&{}&+ \left.\sum_{j=1}^{m} \frac{\partial\mb{\bm$ f$}(t_0,\mb{\bm$ x$}_0,\mb{\bm$ u$}(t_0))}{\partial u_j} \dot{u}_j \right] \frac{(t-t_0)^2}{2}+\cdots
\end{eqnarray*}
If we stop after the $q$th term of the Taylor series solution, then the remainder after $q$ terms of a Taylor series gives us an estimate of the error. We can expect the Taylor series solution to give a fairly accurate approximation of the solution of the state equations only in a small interval about $t_0$. Thus, if $\Delta t =t-t_0$ is some small time interval such that the Taylor series converges for $t \in [t_0,t_0+\Delta t]$, we can use the expansion
\[
\mb{\bm$x$}(t)=\mb{\bm$x$}_0+\frac{d\mb{\bm$ x$}(t_0)}{dt}(t-t_0)+\frac{d^2\mb{\bm$ x$}(t_0)}{dt^2}\frac{(t-t_0)^2}{2}+\ldots .
\]
to evaluate $\mb{\bm$ x$}(t_1)=\mb{\bm$ x$}(t_0+\Delta t)$. Then, we can start over again solving $\dot{\mb{\bm$ x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t))$ with the new initial condition $\mb{\bm$x$}(t_1)$.
We now illustrate the method with simple examples.
\begin{example} \rm
Consider the first order differential equation
\[
\frac{dx}{dt}=t^2-x^2,\quad x(1)=1.
\]
Expanding $x=x(t)$ into a Taylor series, yields
\[
x(t)=x(1)+\frac{dx(1)}{dt}(t-1)+\frac{d^2x(1)}{dt^2} \frac{(t-1)^2}{2}+\ldots .
\]
Since $\frac{dx}{dt}=t^2-x^2$, we obtain by successive differentiations
\begin{eqnarray*}
\frac{d^2x}{dt^2}&=&2t-2x\frac{dx}{dt}, \\
\frac{d^3x}{dt^3}&=&2-2\left( \frac{dx}{dt} \right) ^2-2x\frac{d^2x}{dt^2}, \\
\frac{d^4x}{dt^4}&=&-6\frac{dx}{dt}\frac{d^2x}{dt^2}-2x\frac{d^3x}{dt^3},\\
& & \vdots
\end{eqnarray*}
Evaluating the above at $t=1$, yields
\begin{eqnarray*}
\frac{dx(1)}{dt}&=&0,\\
\frac{d^2x(1)}{dt}&=&2,\\
\frac{d^3x(1)}{dt}&=&-2,\\
\frac{d^4x(1)}{dt}&=&4,\\
& & \vdots
\end{eqnarray*}
Thus,
\[
x(t)=1+(t-1)^2-\frac{(t-1)^3}{3}+\frac{(t-1)^4}{6}+\cdots .
\]
\end{example}
\markex
\begin{example} \rm
We consider now the second-order differential equation
\[
\frac{d^2x}{dt^2}+3\frac{dx}{dt}+6tx=0,
\]
subject to the initial conditions $x(0)=1$, and $\frac{dx(0)}{dt}=1$.
Taking into account that $\frac{d^2x}{dt^2}=-3\frac{dx}{dt}-6tx$, we obtain by successive differentiation
\begin{eqnarray*}
\frac{d^3x}{dt^3}&=&-3\frac{d^2x}{dt^2}-6t\frac{dx}{dt}-6x,\\
\frac{d^4x}{dt^4}&=&-3\frac{d^3x}{dt^3}-6t\frac{d^2x}{dt^2}-12\frac{dx}{dt},\\
\frac{d^5x}{dt^5}&=&-3\frac{d^4x}{dt^4}-6t\frac{d^3x}{dt^3}-18\frac{d^2x}{dt^2},\\
& & \vdots
\end{eqnarray*}
Taking into account the initial conditions, we evaluate all coefficients at $t=0$ to get
\begin{eqnarray*}
\frac{d^2x(0)}{dt^2}&=&-3,\\
\frac{d^3x(0)}{dt^3}&=&3,\\
\frac{d^4x(0)}{dt^4}&=&-21,\\
\frac{d^5x(0)}{dt^5}&=&117,\\
& & \vdots
\end{eqnarray*}
Hence,
\begin{eqnarray*}
x(t)&=&x(0)+\frac{dx(0)}{dt}t+\frac{d^2x(0)}{dt^2} \frac{t^2}{2}+\ldots\\
&=&1+t-\frac{3}{2}t^2+\frac{1}{2}t^3-\frac{7}{8}t^4+\frac{117}{120}t^5+\ldots .
\end{eqnarray*}
\end{example}
\markex
We now illustrate the method of Taylor series when applied to a system of first order differential equations.
\begin{example} \rm
We consider the following system of two first order equations
\begin{eqnarray*}
\frac{dx_1}{dt}&=&f_1(t,x_1,x_2)\\
\frac{dx_2}{dt}&=&f_2(t,x_1,x_2).
\end{eqnarray*}
The initial conditions are at $t=t_0$, $x_1(t_0)=x_{10}$, $x_2(t_0)=x_{20}$. Having the initial conditions, we can determine the values $\frac{dx_1(t_0)}{dt}$ and $\frac{dx_2(t_0)}{dt}$ from the original equations. We then differentiate $\frac{dx_1}{dt}$ and $\frac{dx_2}{dt}$ to obtain
\begin{eqnarray*}
\frac{d^2x_1}{dt^2}&=&\frac{d}{dt}f_1(t,x_1,x_2,\frac{dx_1}{dt},\frac{dx_2}{dt})\\
\frac{d^2x_2}{dt^2}&=&\frac{d}{dt}f_2(t,x_1,x_2,\frac{dx_1}{dt},\frac{dx_2}{dt}).
\end{eqnarray*}
Now, $\frac{d^2x_1}{dt^2}$ and $\frac{d^2x_2}{dt^2}$ can be calculated by substituting known quantities on the right-hand sides of the above equations. By successive differentiation and substitution, we can determine $\frac{d^3x_1}{dt^3}$ and $\frac{d^3x_2}{dt^3}$ as well as higher derivatives at $t=t_0$. The solution to the given system is then computed as
\begin{eqnarray*}
x_1(t)&=&x_{10}+\frac{dx_1(t_0)}{dt}(t-t_0)+\frac{d^2x_1(t_0)}{dt^2}\frac{\left(t-t_0\right)^2}{2}+\frac{d^3x_1(t_0)}{dt^3}\frac{\left(t-t_0\right)^3}{3!}+\ldots \\
x_2(t)&=&x_{20}+\frac{dx_2(t_0)}{dt}(t-t_0)+\frac{d^2x_2(t_0)}{dt^2}\frac{\left(t-t_0\right)^2}{2}+\frac{d^3x_2(t_0)}{dt^3}\frac{\left(t-t_0\right)^3}{3!}+\ldots .
\end{eqnarray*}
\end{example}
\markex
The Taylor series method can be programmed on a computer. However, it is not very efficient for high order differential equations since it requires too much storage. In the following sections we discuss other methods for numerical solution of the state equations that are related to the method of Taylor series. We first present a numerical technique known as {\em Euler's method}.

\section{Euler's Method}
Suppose that we wish to obtain the solution of the model
\[
\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t))
\]
in the time interval $(t_0,t_f)$, subject to the initial condition $\mb{\bm$x$}(t_0)=\mb{\bm$x$}_0$, and the input vector $\mb{\bm$u$}(t)$. We divide the interval $\left[t_0,t_f \right]$ into $N$ equal subintervals of width
\[
h=\frac{t_f-t_0}{N}.
\]
We call $h$ the {\em step length}. We set
\[
t_k=t_0+kh.
\]
To proceed further, we approximate the derivative $\frac{d\mb{\bm$x$}(t)}{dt}$ at time $t_k$ by
\[
\frac{\mb{\bm$x$}(t_{k+1})-\mb{\bm$x$}(t_k)}{h},
\]
and write
\[
\frac{\mb{\bm$x$}(t_{k+1})-\mb{\bm$x$}(t_k)}{h}=\mb{\bm$ f$}(t_k,\mb{\bm$x$}(t_k),\mb{\bm$u$}(t_k)).
\]
Hence, for $k=1,2,\ldots ,N$, we have
\[
\mb{\bm$x$}(t_{k+1})=\mb{\bm$x$}(t_k)+h\mb{\bm$f$}(t_k,\mb{\bm$ x$}(t_k),\mb{\bm$u$}(t_k)).
\]
Thus we see that if $h$ is sufficiently small, we can approximately determine the state at time $t_1=t_0+h$ from the initial condition $\mb{\bm$x$}_0$ to get
\[
\mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+h\mb{\bm$f$}(t_0,\mb{\bm$ x$}_0,\mb{\bm$u$}(t_0)).
\]
Once we have determined the approximate solution at time $t_1$, we can determine the approximate solution at time $t_2=t_1+h$, and so on.
Euler's method can also be arrived at when instead of considering the differential equation $\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$u$}(t))$, we rather start with its equivalent integral representation
\[
\mb{\bm$x$}(t)=\mb{\bm$x$}_0+\int_{t_0}^{t}\mb{\bm$f$}(s,\mb{\bm$x$}(s),\mb{\bm$u$}(s))ds.
\]
We then use the {\em rectangular rule\/} for numerical integration. The rectangular rule can be stated as follows. The area under the curve $z=g(t)$ between $t=t_k$ and $t=t_k+h$ is approximately equal to the area of the rectangle $ABCD$ in Figure~\ref{fig_rectint}, that is
\[
\int_{t_k}^{t_k+h}g(t)dt \approx hg(t_k).
\]
\begin{figure}
\vspace{3.5in}
\hskip1in\special{bmp:c:/mybook/rectint.bmp x=4.67in y=3.5in}
\caption{Illustration of the rectangular rule of integration}
\label{fig_rectint}
\end{figure}
Applying the rectangular rule of integration to $\mb{\bm$ x$}(t)=\mb{\bm$x$}_0+\int_{t_0}^{t}\mb{\bm$f$}(s,\mb{\bm$ x$}(s),\mb{\bm$u$}(s))ds$, we obtain Euler's method.
Euler's method is rather elementary, and thus, not as accurate as some of the more sophisticated techniques that we are going to discuss next.
\section{Predictor-Corrector Method}
The predictor-corrector method, also known as Heun's method, is based on the trapezoid rule for numerical integration. The trapezoid rule states that the area under the curve $z=g(t)$ between $t=t_k$ and $t=t_k+h$ is approximately equal to the area of the trapezoid $ABCD$ as depicted in Figure~\ref{fig_trapint}, that is
\[
\int_{t_k}^{t_k+h}g(t)dt\approx \frac{g(t_k)+g(t_k+h)}{2}h.
\]
\begin{figure}
\vspace{3.5in}
\hskip1in\special{bmp:c:/mybook/trapint.bmp x=4.67in y=3.5in}
\caption{The trapezoid rule of integration}
\label{fig_trapint}
\end{figure}
Let
\[
\dot{\mb{\bm$x$}}(t)=\mb{\bm$f$}(t,\mb{\bm$x$}(t),\mb{\bm$ u$}(t)),\quad \mb{\bm$x$}(t_0)=\mb{\bm$x$}_0.
\]
Suppose we wish to find $\mb{\bm$x$}(t)$ at $t=t_{f} > t_0$. We first divide the interval $\left[t_0,t_{f} \right]$ into $N$ equal subintervals using evenly spaced subdivision points to obtain
\[
t_0 < t_1 < t_2 < \ldots < t_{N}=t_{f}.
\]
Let
\begin{eqnarray*}
t_{k+1}-t_{k} &=&h,\\
\mb{\bm$x$}^{*}_1 &=& \mb{\bm$x$}_0+\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0))\,h,\\
\mb{\bm$m$}_0 &=& \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)),\\
\mb{\bm$m$}_1 &=& \mb{\bm$f$}(t_1,\mb{\bm$x$}^{*}_1,\mb{\bm$u$}(t_1)).
\end{eqnarray*}
Then,
\[
\mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+\frac{h}{2}\left(\mb{\bm$m$}_0 + \mb{\bm$m$}_1 \right).
\]
Once we have found $\mb{\bm$x$}(t_1)$ we can apply the same procedure to find $\mb{\bm$x$}(t_2)$, and so on. Note that $\mb{\bm$m$}_1$ we used to find $\mb{\bm$x$}(t_1)$ is a first approximation to $\mb{\bm$x$}(t_1)$ in the trapezoidal integration formula
\[
\mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+\frac{h}{2}\left(\mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)) + \mb{\bm$f$}(t_1,\mb{\bm$x$}(t_1),\mb{\bm$u$}(t_1) \right).
\]
If the function $\mb{\bm$f$}$ is nonlinear, then we will, in general, not be able to solve the above equation for $\mb{\bm$x$}(t_1)$ exactly. We attempt to obtain $\mb{\bm$x$}(t_1)$ iteratively. A first approximation to $\mb{\bm$x$}(t_1)$ is
\[
\mb{\bm$x$}^{*}_1 = \mb{\bm$x$}_0 + \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0))\,h,
\]
and is called a {\em predictor\/}. The next approximation is
\[
\mb{\bm$x$}(t_1)=\mb{\bm$x$}_0+\frac{h}{2}\left(\mb{\bm$m$}_0 + \mb{\bm$m$}_1 \right),
\]
and is called a {\em corrector}.
\begin{example}
\label{ex:pred_cor}
\rm
We use the predictor-corrector method to find $x(1)$, where the step length $h=0.5$, for the following differential equation
\[
\frac{dx}{dt}=tx^2-x,\quad x(0)=1.
\]
The calculations are summarized in Table~\ref{tab:pred_cor}.
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}\hline
$t_0$ & $x_0$ & $t_1$ & $m_0$ & $x^{*}_1$ & $m_1$ & $m=\frac{m_0+m_1}{2}$ & $\Delta x$\\ \hline 
$0$ & $1$ & $0.5$ & $-1.0$ & $0.5$ & $-0.375$ & $-0.688$ & $-0.344$\\ \hline
$0.5$ & $0.656$ & $1.0$ & $-0.441$ & $0.435$ & $-0.246$ & $-0.344$ & $-0.172$\\ \hline
$1.0$ & $0.484$ &{}&{}&{}&{}&{}&{}\\\hline 
\end{tabular}
\end{center}
\caption{Solving the differential equation using the predictor-corrector method.}
\label{tab:pred_cor}
\end{table}
\end{example}
\markex
We could use $\mb{\bm$x$}(t_1)$ obtained from the corrector formula to construct a next approximation to $\mb{\bm$x$}(t_1)$. In general, we can use the formula
\[
\mb{\bm$x$}^{(k)}(t_1)=\mb{\bm$x$}_0 + \frac{h}{2} \left( \mb{\bm$f$}(t_0,\mb{\bm$x$}_0,\mb{\bm$u$}(t_0)) + \mb{\bm$f$}(t_1,\mb{\bm$x$}^{(k-1)}(t_1),\mb{\bm$u$}(t_1) \right)
\]
to evaluate $\mb{\bm$x$}(t_1)$ until two successive iterations agree to the desired accuracy. This general algorithm is referred to as a second-order predictor-corrector method.

An $n$-dimensional system of first-order differential equations can be worked analogously. We illustrate the procedure for $n=2$, where
\begin{eqnarray*}
\frac{dx}{dt} &=& f(t,x,y), \quad x(t_0)=x_0\\
\frac{dy}{dt} &=& g(t,x,y), \quad y(t_0)=y_0.
\end{eqnarray*}
Let
\begin{eqnarray*}
m_0 &=& f(t_0,x_0,y_0),\\
n_0 &=& g(t_0,x_0,y_0).
\end{eqnarray*}
Then, the predictor equations are
\begin{eqnarray*}
x^{*}_1 &=& x_0 +m_0\,h,\\
y^{*}_1 &=& y_0 + n_0\,h.
\end{eqnarray*}
Let
\begin{eqnarray*}
m_1 &=& f(t_1,x^{*}_1,y^{*}_1),\\
n_1 &=& g(t_1,x^{*}_1,y^{*}_1).
\end{eqnarray*}
Then, the corrector equations are
\begin{eqnarray*}
x(t_1) &=& x_0 + \frac{m_0+m_1}{2}\,h,\\
y(t_1) &=& y_0+\frac{n_0+n_1}{2}\,h.
\end{eqnarray*}
Using the above one can easily generate the corresponding formulas for $t_k$.

We will now show that the predictor-corrector method yields a solution that agrees with the Taylor series solution through terms of degree two. We consider the scalar case. First observe that $m_1=f\left(t_1,x^{*}_1 \right)$ can be written as
\begin{eqnarray*}
f\left(t_1,x^{*}_1 \right) &=& f \left( t_0 + h,x_0+m_0h \right)\\
&\approx& f\left(t_0,x_0 \right) + \left[\begin{array}{cc}
f_{t} & f_{x}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]+ \frac{1}{2} \left[\begin{array}{cc}
h & m_0h
\end{array}\right] \left[\begin{array}{cc}
f_{tt} & f_{tx}\\
f_{xt} & f_{xx}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]. 
\end{eqnarray*}
Substituting the above into the corrector equations yields
\begin{eqnarray*}
x(t_1) &=& x_0 + \frac{h}{2}\left( m_0 + m_1 \right)\\
&=& x_0 + \frac{h}{2}\left( m_0 + m_0\right)+ \frac{h}{2}\left(\left[\begin{array}{cc}
f_{t} & f_{x}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]\right.\\
&+& \left.\frac{1}{2} \left[\begin{array}{cc}
h & m_0h
\end{array}\right] \left[\begin{array}{cc}
f_{tt} & f_{tx}\\
f_{xt} & f_{xx}
\end{array}\right] \left[\begin{array}{c}
h\\
m_0h
\end{array}\right]\right)\\
&=& x_0+m_0h+\left(f_{t}+f_{x}m_0\right)\frac{h^2}{2}+\left(f_{tt}+2f_{tx}m_0 + f_{xx}m_0^2 \right)\frac{h^3}{4}.
\end{eqnarray*}
The Taylor series expansion, on the other hand, gives
\begin{eqnarray*}
x\left(t_1 \right) &=& x_0 + \dot{x}\,h + \ddot{x}\, \frac{h^2}{2} + x^{(3)}\frac{h^3}{6} + \cdots\\
&=& x_0 + f(t_0,x_0) \,h + \left(f_{t} + f_{x}\dot{x} \right) \frac{h^2}{2}\\
{}&{}&+ \left(f_{tt}+2f_{tx} \dot{x} +f_{xx}\left(\dot{x}\right)^2+f_{x}\ddot{x} \right)\frac{h^3}{6} + \cdots,
\end{eqnarray*}
where $f_{t}$, $f_{x}$, $\dot{x}$, etc.\ are all evaluated at $t=t_0$. Comparing the right-hand sides of the last two equations, we see that the solution obtained using the predictor-corrector method agrees with the Taylor series expansion of the true solution through terms of degree two.

\section{Runge's Method}
Runge's method is an improvement over the predictor-corrector method in the sense that it agrees with the Taylor series solution up to terms of degree three. The method depends upon integration by means of Simpson's rule. Simpson's rule says that the area $A$ under the curve $x=g(t)$ between $t_0$ and $t_0+h$ is given approximately by
\[
A=\int_{t_0}^{t_0+h}g(t)dt \approx \frac{h}{6}\left[x_0+4x_1+x_2 \right],
\]
where $x_0=g(t_0)$, $x_1=g\left( t_0 + h/2 \right)$, and $x_2=g(t_0+h)$. The above formula is derived by passing a parabola through the three points $(t_0,x_0)$, $\left(t_0+h/2,x_1\right)$, and $(t_0+h,x_2)$, and assuming that the area under the given curve is approximately equal to the area under the parabola as it is illustrated in Figure~\ref{fig:simpson}.
\begin{figure}
\vskip3.5in
\hskip1in\special{bmp:c:/mybook/simpson.bmp x=4.67in y=3.5in}
\caption{Integration using Simpson's rule}
\label{fig:simpson}
\end{figure}
To verify Simpson's formula we translate the axes so that in the new coordinates $\left( \tilde{t},\tilde{x} \right)$ we have $t=t_0+h/2+\tilde{t}$, $x=\tilde{x}$. In the new coordinates the three points become
\[
\left(-h/2,x_0 \right), \quad \left(0,x_1 \right),\quad \left(h/2,x_2 \right).
\]
The parabola through the above points is described by the equation
\[
\tilde{x}=a\tilde{t}^2 + b\tilde{t} + c.
\]
Hence
\begin{eqnarray*}
x_0 &=& a\frac{h^2}{4}-b\frac{h}{2}+c,\\
x_1 &=& c,\\
x_2 &=& a\frac{h^2}{4}+b\frac{h}{2}+c.
\end{eqnarray*}
Using the above three equations we calculate
\begin{eqnarray*}
c &=& x_1,\\
b &=& \frac{x_2 -x_0}{h},\\
a &=& \frac{2}{h^2}\left( x_0 -2x_1 + x_2 \right).
\end{eqnarray*}
Then, the area under the parabola is
\begin{eqnarray*}
A &=& \int_{-h/2}^{h/2}\left( a\tilde{t}^2 + b\tilde{t} +c \right) d\tilde{t}\\
&=& \int_{-h/2}^{h/2}\left(\frac{2}{h^2}\left( x_0 -2x_1 + x_2 \right)\tilde{t}^2 + \frac{x_2 -x_0}{h}\tilde{t} +x_1 \right)d\tilde{t}\\
&=& \frac{h}{6}\left( x_0 +4x_1 + x_2 \right).
\end{eqnarray*}

If we are given the differential equation
\[
\frac{dx}{dt}=f(t,x),\quad x(t_0)=x_0,
\]
we integrate it to obtain
\[
x(t_1)=x(t_0)+\int_{t_0}^{t_1}f(s,x(s))ds.
\]
Then, we approximate the integral
\[
\int_{t_0}^{t_1}f(s,x(s))ds
\]
by
\[
\frac{h}{6}\left( m_0 + 4m_1 +m_3 \right),
\]
where
\begin{eqnarray*}
h &=& t_1-t_0,\\
m_0 &=& f(t_0,x_0),\\
m_1 &=& f \left( t_0+\frac{h}{2},x_0+m_0\frac{h}{2} \right),\\
m_2 &=& f\left( t_0+h,x_0+m_0h \right),\\
m_3 &=& f\left( t_0+h,x_0 + m_2h \right).
\end{eqnarray*}
Thus, we may write
\[
x(t_1)=x_0+ \Delta x,
\]
where $\Delta x = mh$, and $m=\frac{1}{6}\left( m_0 + 4m_1 + m_3\right)$.
\begin{example} \rm
\label{ex:runge}
We use Runge's method to find $x(1)$, where the step length $h=0.5$, for the differential equation 
\[
\frac{dx}{dt}=tx^2-x,\quad x(0)=1.
\]
We found $x(1)$ using the predictor-corrector method in Example~\ref{ex:pred_cor}. The calculations for Runge's method are summarized in Table~\ref{tab:runge}.
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}\hline
$t_0$ & $x_0$ & $m_0$ & $m_1$ & $m_2$ & $m_3$ & $m$ & $\Delta x$\\ \hline 
$0$ & $1$ & $-1$ & $-0.610$ & $-0.375$ & $-0.482$ & $-0.654$ & $-0.327$\\ \hline
$0.5$ & $0.673$ & $-0.446$ & $-0.325$ & $-0.248$ & $-0.248$ & $-0.332$ & $-0.166$\\ \hline
$1.0$ & $0.507$ &{}&{}&{}&{}&{}&{}\\\hline 
\end{tabular}
\end{center}
\caption{Solving the differential equation using Runge's method.}
\label{tab:runge}
\end{table}
\end{example}
\markex
Runge's method can, of course, be used for systems of first order differential equations. If, for example,
\begin{eqnarray*}
\frac{dx}{dt} &=& f(t,x,y),\quad x(t_0)=x_0,\\
\frac{dy}{dt} &=& g(t,x,y),\quad y(t_0)=y_0,
\end{eqnarray*}
then
\begin{eqnarray*}
x(t_1) &=& x_0 + \Delta x,\\
y(t_1) &=& y_0 + \Delta y,
\end{eqnarray*}
where
\begin{eqnarray*}
\Delta x &=& \frac{h}{6} \left( m_0 + 4m_1 + m_3 \right),\\
\Delta y &=& \frac{h}{6} \left( n_0 + 4n_1 + n_3 \right).
\end{eqnarray*}
In the above expression
\begin{eqnarray*}
m_0 &=& f(t_0,x_0,y_0),\\
m_1 &=& f\left( t_0 + h/2,x_0+m_0\,h/2,y_0+n_0\,h/2 \right),\\
m_2 &=& f \left( t_0+h,x_0+m_0h,y_0+n_0h \right),\\
m_3 &=& f \left( t_0 +h,x_0+m_2h,y_0+n_2h \right),\\
n_0 &=& g(t_0,x_0,y_0),\\
n_1 &=& g\left( t_0 + h/2,x_0+m_0\,h/2,y_0+n_0\,h/2 \right),\\
n_2 &=& g\left( t_0+h,x_0+m_0h,y_0+n_0h \right),\\
n_3 &=& g \left( t_0 +h,x_0+m_2h,y_0+n_2h \right).
\end{eqnarray*}
\section{The Runge-Kutta Method}
This method is an improvement over the Runge method, and agrees with the Taylor series solution up to terms of degree four. It is based on the use of the formula
\[
x(t_1)=x(t_0)+\frac{h}{6}\left( m_0+2m_1 + 2m_2 + m_3 \right),
\]
where
\begin{eqnarray*}
m_0 &=& f(t_0,x_0),\\
m_1 &=& f\left( t_0 + h/2, x_0 + m_0\,h/2 \right),\\
m_2 &=& f \left( t_0 + h/2, x_0 + m_1\,h/2 \right),\\
m_3 &=& f \left( t_0 + h, x_0 +h\,m_2 \right).
\end{eqnarray*}
Thus, we may write
\[
x(t_1)=x_0+ \Delta x,
\]
where $\Delta x = mh$, and $m=\frac{1}{6}\left( m_0 + 2m_1 + 2m_2 + m_3\right)$.
\begin{example} \rm
We use the Runge-Kutta method to find $x(1)$, where the step length $h=0.5$, for the differential equation we used in Example~\ref{ex:pred_cor} and Example~\ref{ex:runge}.
Our calculations for the Runge-Kutta method are summarized in Table~\ref{tab:runge_kutta}.
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}\hline
$t_0$ & $x_0$ & $m_0$ & $m_1$ & $m_2$ & $m_3$ & $m$ & $\Delta x$\\ \hline 
$0$ & $1$ & $-1.0$ & $-0.609$ & $-0.668$ & $-0.444$ & $-0.666$ & $-0.333$\\ \hline
$0.5$ & $0.667$ & $-0.444$ & $-0.324$ & $-0.328$ & $-0.250$ & $-0.333$ & $-0.167$\\ \hline
$1.0$ & $0.500$ &{}&{}&{}&{}&{}&{}\\\hline 
\end{tabular}
\end{center}
\caption{Solving the differential equation using the Runge-Kutta method.}
\label{tab:runge_kutta}
\end{table}
\end{example}
\markex
The above problem has been worked by the predictor-corrector, Runge's, and the Runge-Kutta methods. This problem can be solved exactly by means of the Taylor series. The Taylor series solution is
\[
x(t) = 1-t + t^2 - t^3+\cdots
\]
The series converges to $1/(1+t)$ for $-1 < t < 1$. A check on the differential equation shows that
\[
x(t) = \frac{1}{1+t}
\]
is indeed a solution for all $t> -1$. Thus, $x(1)=0.500$.

\section{Exercises}
\begin{exercise}\rm
Show
\end{exercise}
