# Introduction to Numerical Methods/Ordinary Differential Equations

Resources:

-   notes on Ordinary Differential
    Equations
-   notes on Euler\'s
    Method
-   notes on Runge-Kutta 2nd Order
    Method
-   notes on Runge-Kutta 4th Order
    Method

\_\_TOC\_\_

A differential
equation is a
equation that relates some function with its derivatives. Such
relationship is commonly found in dynamic system modeling where the rate
of change of a quantity is affect by another quantity, which can be a
function of other quantities. Therefore differential equations are often
used to mathematically represent such relations in many disciplines
including engineering, physics, economics, and biology. We will focus on
ordinary differential equations, which involves no partial derivative.

The solution to a differential equation is the function or a set of
functions that satisfies the equation. Some simple differential
equations with explicit formulas are solvable analytically, but we can
always use numerical methods to estimate the answer using computers to a
certain degree of accuracy.

<http://martinfowler.com/bliki/TwoHardThings.html> \"Two hard things in
computer science: cache invalidation, naming things, and off-by-one
errors.\"

# Initial Value Problem

The general form of first order differential equation: $y'=f(x, y)$ or
$y'=f(x, y(x))$. The solution contains an arbitrary constant (the
constant of integration) because $z(x)=y(x)+c$ satisfies all
differential equations $y(x)$ satisfies. With an auxiliary condition,
e.g. y(a)=b, we can solve y\'=F(x, y). This is known as the initial
value problem.

An indefinite integral of a function $f(x)$ is a class of functions
$F(x)$ so that $F'(x)=f(x)$. The indefinite integral is also known as
the antiderivative. The definite integral of a function over a fixed
interval is a number.

# Connection with Integration

To calculate the integral $F(x)=\int_{a}^{x}f(t)dt$ is to find a $F(x)$
that satisfies the condition $\frac{dF(x)}{dx}=f(x)$, $F(a)=0$. In other
words, integration can be used to solve differential equations with a
initial condition.

# Euler\'s Method

!Illustration of the Euler method. The unknown curve is in blue, and
its polygonal approximation is in
red.
Euler\'s method uses a
step-wise approach to solve ordinary differential equations with an
initial condition.

Given $\frac{dy}{dx} = f(x,y)$ and $y(x_{0})=y_{0}$ ($x$ is an
independent variable and $y$ is a dependent variable) Euler\'s methods
solves for $y_{i}$ for $x=x_i$.

Euler\'s method can be derived in a number of ways. Lets look at a
geographical description.

Given that $\frac{dy}{dx} = f(x,y)$ at a particular point
$(x_{i}, y_{i})$ we can approximate $f(x_{i}, y_{i})$ with
$\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}$. The derivative represents the
instantaneous (happening continuously) change in $y$ at $x$:

$$\lim_{\Delta x \to 0}\frac{\Delta y}{\Delta x}=\frac{dy}{dx}=y'(x)$$
The average change in $y$ over $x$ is

:   $\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}=\frac{\Delta y}{\Delta x}$,

which represents a discrete version of the change (changes over
\"steps\" of time). If the change in $x$ is small enough, either of the
two can be used to approximate the other.

Assume that we know the derivative function of an unknown curve and
wants to find the shape of the curve (a function that satisfies a
differential equation). We can start from a known point, use the
differential equation as as a formula to get the slope of the tangent
line to the curve at any point on the curve, and take a small step along
that tangent line up to the next point. If the step is small enough, we
can expect the next point to be close to the curve. If we pretend that
point is still on the curve, the same process can be used to find the
next point and so on. Graphically we are using a polynomial to
approximate the unknown curve. The gap between the two curves represents
the error of the approximation, which can be reduced by shortening the
step sizes.

We just derived the formula for Euler\'s method:

$$y_{i+1} = y_{i} + f(x_{i},y_{i})(x_{i+1}-x_{i})$$

## Example

Given $\frac{dy}{dx} = y$, and $y(0)=1$, solve it for $x$ at 1, i.e.
y(1).

The Euler\'s method formula is a recursively defined function, which can
be calculated iteratively. $(x_{i+1}-x_{i})$ is called the step size. If
we pick a step size of 1 the solution is as follows: !Illustration of
numerical integration for the equation $y'=y, y(0)=1.$ Blue is the Euler
method; green, the midpoint method; red,
the exact solution, $y=e^t.$ The step size is
*h* = 1.0.=1. Blue is the Euler method; green, the midpoint method; red, the exact solution, y=e^t. The step size is h = 1.0."){width="200"}

:   {\| class=\"wikitable\"

\|- ! $i$ !! $x_{i}$ !! $y_i$ !!$dy/dx$ \|- \| 0 \|\| 0 \|\| 1 \|\| 1
\|- \| 1 \|\| 1 \|\| 2 \|\| 2 \|- \| 2 \|\| 2 \|\| 4 \|\| 4 \|- \| 3
\|\| 3 \|\| 8 \|\| 8 \|}

Another way to derive Euler\'s method is to use taylor expansion around
$x_{i}$:

$$y(x_0 + h) = y(x_0) + h y'(x_0) + \frac{1}{2}h^2 y''(x_0) + O(h^3).$$
If we ignore higher order derivative terms we get the Euler\'s method
formula:

$$y(x_0 + h) = y(x_0) + h y'(x_0)$$

# Runge-Kutta 2nd Order Method

One of the Runge-Kutta 2nd order method is the midpoint
method, which is a
modified Euler\'s method (one-step method) for numerically solving
ordinary differential equations:

$$y'(x) = f(x, y(x)), \quad y(x_0) = y_0$$. The midpoint method is given
by the formula

$$y_{i+1} = y_i + hf\left(x_i+\frac{h}{2},y_i+\frac{h}{2}f(x_i, y_i)\right)$$.
Note that $y_{i}$ is the approximation of $y(x_{i})$ using this method.
!Illustration of the midpoint method assuming that $y_n$ equals the
exact value $y(t_n).$ The midpoint method computes $y_{n+1}$ so that the
red chord is approximately parallel to the tangent line at the midpoint
(the green
line).. The midpoint method computes y_{n+1} so that the red chord is approximately parallel to the tangent line at the midpoint (the green line).")

# Runga-Kutta 4th Order Method

The Runga-Kutta 4th order
method is
often called \"RK4\", \"classical Runge--Kutta method\" or simply as
\"the Runge--Kutta method\".

$$\begin{align}
y_{i+1} &= y_i + \tfrac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4 \right)\\
x_{i+1} &= x_i + h \\
\end{align}$$

for *i* = 0, 1, 2, 3, . . . , using

$$\begin{align}
k_1 &= f(x_i, y_i),
\\
k_2 &= f(x_i + \tfrac{h}{2}, y_i + \tfrac{1}{2} k_1 h),
\\
k_3 &= f(x_i + \tfrac{h}{2}, y_i + \tfrac{1}{2} k_2 h),
\\
k_4 &= f(x_i + h, y_i + k_3 h).
\end{align}$$ This method calculate $y_{i+1}$ by adding to $y_i$) a
weighted average of four increments, each of which is the step size $h$
multiplied by an estimated slope:

-   $k_1$ is the increment based on the estimated slope at the beginning
    of the interval (Euler\'s method) ;
-   $k_2$ is the increment based on the slope at the midpoint of the
    interval;
-   $k_3$ is again the increment based on the slope at the midpoint, but
    now using $k_2$;
-   $k_4$ is the increment based on the slope at the end of the
    interval, using $k_3$.

In averaging the four increments, greater weight is given to the
increments at the midpoint. The weights are chosen such that if $f$ is
independent of $y$, so that the differential equation is equivalent to a
simple integral.

# Case Study

source

solution in a iPython
notebook

