CLP-2 Integral Calculus

Our goal is to approximate (numerically) the unknown function

for $t\ge t_0\text{.}$ We are told explicitly the value of $y(t_0)\text{,}$ namely $y_0\text{.}$ So we know $f(\tau ,y(\tau )){|}_{\tau =t_0}=f(t_0,y_0)\text{.}$ But we do not know the integrand $f(\tau ,y(\tau ))$ for $\tau >t_0\text{.}$ On the other hand, if $\tau$ is close $t_0\text{,}$ then $f(\tau ,y(\tau ))$ will remain close

²

This will be the case as long as $f(t,y)$ is continuous.

to $f(t_0,y_0)\text{.}$ So pick a small number $h$ and define

By the above argument

Now we start over from the new point $(t_1,y_1)\text{.}$ We now know an approximate value for $y$ at time $t_1\text{.}$ If $y(t_1)$ were exactly $y_1\text{,}$ then the instantaneous rate of change of $y$ at time $t_1\text{,}$ namely $y^{\prime }(t_1)=f(t_1,y(t_1))\text{,}$ would be exactly $f(t_1,y_1)$ and $f(t,y(t))$ would remain close to $f(t_1,y_1)$ for $t$ close to $t_1\text{.}$ Defining

we have

We just repeat this argument ad infinitum. Define, for $n=0,1,2,3,\cdots$

Suppose that, for some value of $n\text{,}$ we have already computed an approximate value $y_n$ for $y(t_n)\text{.}$ Then the rate of change of $y(t)$ for $t$ close to $t_n$ is $f(t,y(t))\approx f(t_n,y(t_n))\approx f(t_n,y_n)$ and

Here is a table applying a few steps of Euler’s method to the initial value problem

with step size $h=0.1\text{.}$ For this initial value problem

Of course this initial value problem has been chosen for illustrative purposes only. The exact solution is

³

Even if you haven’t learned how to solve initial value problems like this one, you can check that $y(t)=2+2t+e^t$ obeys both $y^{\prime }(t)=-2t+y(t)$ and $y(0)=3\text{.}$

$y(t)=2+2t+e^t\text{.}$

Euler’s method is one algorithm which generates approximate solutions to the initial value problem

In applications, $f(t,y)$ is a given function and $t_0$ and $y_0$ are given numbers. The function $y(t)$ is unknown. Denote by $\phi (t)$ the exact solution

⁴

Under reasonable hypotheses on $f\text{,}$ there is exactly one such solution. The interested reader should search engine their way to the Picard-Lindelöf theorem.

for this initial value problem. In other words $\phi (t)$ is the function that obeys

exactly.

Fix a step size $h$ and define $t_n=t_0+nh\text{.}$ By turning the problem into one of approximating integrals, we now derive another algorithm that generates approximate values for $\phi$ at the sequence of equally spaced time values $t_0,t_1,t_2,\cdots \text{.}$ We shall denote the approximate values $y_n$ with

By the fundamental theorem of calculus and the differential equation, the exact solution obeys

Fix any $n$ and suppose that we have already found $y_0,y_1,\cdots ,y_n\text{.}$ Our algorithm for computing $y_{n+1}$ will be of the form

In Euler’s method, we approximated $f(t,\phi (t))$ for $t_n\le t\le t_{n+1}$ by the constant $f(t_n,y_n)\text{.}$ Thus

So Euler approximates the area of the complicated region $0\le y\le f(t,\phi (t))\text{,}$ $t_n\le t\le t_{n+1}$ (represented by the shaded region under the parabola in the left half of the figure below) by the area of the rectangle $0\le y\le f(t_n,y_n)\text{,}$ $t_n\le t\le t_{n+1}$ (the shaded rectangle in the right half of the figure below).

Recall that we have already found $y_0,\cdots ,y_n$ and are in the process of finding $y_{n+1}\text{.}$ So we already have an approximation for $\phi (t_n)\text{,}$ namely $y_n\text{.}$ But we still need to approximate $\phi (t_{n+1})\text{.}$ We can do so by using one step of the original Euler method! That is

So our approximation of $\frac{1}{2}[f(t_n,\phi (t_n))+f(t_{n+1},\phi (t_{n+1}))]$ is

and

Putting everything together

⁵

Notice that we have made a first approximation for $\phi (t_{n+1})$ by using Euler’s method. Then improved Euler uses the first approximation to build a better approximation for $\phi (t_{n+1})\text{.}$ Building an approximation on top of another approximation does not always work, but it works very well here.

, the improved Euler’s method algorithm is

Here are the first two steps of the improved Euler’s method applied to

with $h=0.1\text{.}$ In each step we compute $f(t_n,y_n)\text{,}$ followed by $y_n+f(t_n,y_n)h\text{,}$ which we denote ${\tilde{y}}_{n+1}\text{,}$ followed by $f(t_{n+1},{\tilde{y}}_{n+1})\text{,}$ followed by $y_{n+1}=y_n+\frac{1}{2}[f(t_n,y_n)+f(t_{n+1},{\tilde{y}}_{n+1})]h\text{.}$

Here is a table which gives the first five steps.

The Runge-Kutta

⁶

Carl David Tolmé Runge (1856–1927) and Martin Wilhelm Kutta (1867–1944) were German mathematicians.

algorithm is similar to the Euler and improved Euler methods in that it also uses, in the notation of the last subsection,

But rather than approximating ${\int }_{t_n}^{t_{n+1}}f(t,\phi (t))\mathrm{d}t$ by the area of a rectangle, as does Euler, or by the area of a trapezoid, as does improved Euler, it approximates by the area under a parabola. That is, it uses Simpson’s rule. According to Simpson’s rule (which is derived in Section 1.11.3)

Analogously to what happened in our development of the improved Euler method, we don’t know $\phi (t_n)\text{,}$ $\phi (t_n+\frac{h}{2})$ or $\phi (t_n+h)\text{.}$ So we have to approximate them as well. The Runge-Kutta algorithm, incorporating all these approximations, is

⁷

It is well beyond our scope to derive this algorithm, though the derivation is similar in flavour to that of the improved Euler method. You can find more in, for example, Wikipedia.

That is, Runge-Kutta uses

Here are the first two steps of the Runge-Kutta algorithm applied to

with $h=0.1\text{.}$

Now, while this might look intimidating written out in full like this, one should keep in mind that it is quite easy to write a program to do this. Here is a table giving the first five steps. The data is only given to three decimal places even though the computation has been done to many more.