Numerical Methods/Solution of IVP

Suppose you are given $y'(x)=f(x,y)$ where $y$, the dependent variable, is a function of the independent variable $x$ and $y(x_{0})=y_{0}$ is given. This is an initial value problem of ODE's because it specifies the initial condition(s) and the differential equation giving $y$. The problem is to calculate the values of $y$ at points $x>x_{0}$. There are a variety of numerical methods to solve this type of problem. They yield a numerical solution which is nothing but a series of values $y_{n}$ corresponding to the points $x_{n}$.

The numerical solution is an approximate solution. The real solution denoted by $y(x)$ could either be

• impossible to get because the function f may be non integrable
• too difficult to get

Anyways, a numerical solution gives us an algorithm to compute $y_{n}$ based on already available information. So one could write a computer program to compute it. This is the basis of simulations.

Euler's Method is the simplest, and its somewhat similar to calculating the integral numerically. It is easy to see what is happening on a graph. At each point x_n you are just looking at what direction the derivative tells you to go, and you take a step in that direction.

Euler's Explicit Method

• Let $x_{n+1}=x_{n}+h$, where h is some small step value.
• Let $y_{n+1}=y_{n}+h\cdot f(x_{n},y_{n})$.

Heun's second order or explicit trapezoidal method

• Let $x_{n+1}=x_{n}+h$, where h is some small step value.
• Let ${\tilde {y}}_{n+1}=y_{n}+h\cdot f(x_{n},y_{n})$.
• Let $y_{n+1}=y_{n}+(h/2)\cdot [f(x_{n},y_{n})+f(x_{n+1},{\tilde {y}}_{n+1})]$.

Examples

Now we will go over some examples of using this algorithm.

$y'=x$

Suppose we know also that we want to calculate f(1) using n=4 iterations. Thus h=1/4=.25.

• x_0=0, y_0=5
• x_1=.25, y_1=5+.25*f'(0,5)=5+.25*0=5
• x_2=.5, y_2=5+.25*f'(.25,5)=5+.25*.25=5+.0625=5.0625
• x_3=.75, y_3=5.0625+.25*f'(.5,5.0625)=5.0625+.25*.5=5.0625+.125=5.1875
• x_4=1, y_4=5.1875+.25*f'(.75,5.1875)=5.1875+.25*.75=5.1875+.1875=5.375

So we are done, and we know that f(1)~=5.375.

We know the exact solution by integrating, y'=x, so y=x^2/2+c. In our case, to calculate c, we would substitute 0 for x and 5 for y to see that c=5. so f(1)=1^2/2+5 = 5.5.

Error analysis

The error only depends on the concavity of the function (the second derivative). Thus we can estimate the error as O(x^2). We can decrease the error by taking smaller step sizes.

Runge-Kutta Methods

Runge-Kutta methods are a family of algorithms particularly well-suited for numerically solving systems of simultaneous differential equations.

Suppose we have a pair of simultaneous equations of the form

${\frac {dx}{dt}}=f(t,x,y)$, and

${\frac {dy}{dt}}=g(t,x,y)$,

and we wish to numerically integrate these with respect to t, with a step size of $\Delta t$. We then compute the following quantities:

${\begin{array}{lcl}m_{1}&=&f(t,x,y)\\n_{1}&=&g(t,x,y)\\m_{2}&=&f(t+{\frac {\Delta t}{2}},x+{\frac {m_{1}}{2}},y+{\frac {n_{1}}{2}})\\n_{2}&=&g(t+{\frac {\Delta t}{2}},x+{\frac {m_{1}}{2}},y+{\frac {n_{1}}{2}})\\m_{3}&=&f(t+{\frac {\Delta t}{2}},x+{\frac {m_{2}}{2}},y+{\frac {n_{2}}{2}})\\n_{3}&=&g(t+{\frac {\Delta t}{2}},x+{\frac {m_{2}}{2}},y+{\frac {n_{2}}{2}})\\m_{4}&=&f(t+\Delta t,x+m_{3},y+n_{3})\\n_{4}&=&g(t+\Delta t,x+m_{3},y+n_{3})\\\Delta x&=&{\frac {1}{6}}(m_{1}+2m_{2}+2m_{3}+m_{4})\Delta t\\\Delta y&=&{\frac {1}{6}}(n_{1}+2n_{2}+2n_{3}+n_{4})\Delta t\\\end{array}}$

$x,y$ and $t$ are then incremented by $\Delta x,\Delta y$ and $\Delta t$ respectively, and all the above quantities calculated again for the next step.

This method can obviously be extended to three or more simultaneous equations, and the error is of order $\Delta t^{5}$. The high accuracy and ease of implementation make this method a popular computational tool.

Also useful is the fact that $\Delta t$ can be changed after any step of the algorithm, which allows rapid changes in the variables to be dealt with.