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

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.

$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})+hy'(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})+hy'(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).

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{aligned}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{aligned}}$

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

${\begin{aligned}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{aligned}}$

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