Consider the situation in which the solution, *y(t)*, to a differential equation

is to be approximated by computer (starting from some known initial condition, *y(t₀)=y₀*; also, note that the tick mark denotes differentiation). The following text develops an intuitive technique for doing so, and presents some examples. This technique is known as "Second Order Runge-Kutta".

The first order Runge-Kutta method used the derivative at time *t₀* (*t₀*=0 in the graph below) to estimate the value of the function at one time step in the future. If you are not familiar with it, you should read the section entitled: *A First Order Linear Differential Equation with No Input*. We repeat the central concept of generating a step forward in time in the following text.

with the initial condition set as *y(0)=3*. The exact solution in this case is *y(t)=3e ^{-2t}*,

In the graph below, the slope at *t=0* is called *k _{1}*, and the estimate is called

This obviously leads to some error in the estimate, and we would like to reduce this error. One way we could do this, conceptually, is to use the derivative at the halfway point between *t=0* and *t=h=0.2*. The slope at this point (*t=½h=0.1*) is shown below (and is labeled* k _{2}*). Note the line (orange) is tangent to the curve (blue) at

Now if we use this intermediate slope, *k _{2}*, as we step ahead in time then we get better estimate,

This seems like a very nice solution, and obviously generates a significantly more accurate approximation than the first order technique that uses a line with slope, *k _{1}*, calculated at

What we do instead is use the First Order Runge-Kutta to generate an approximate value for *y(t)* at *t=½h=0.1*, call it y_{1}(*½h*). We then use this estimate to generate k_{2} (which will be an approximation to the slope at the midpoint), and then use k_{2} to find y*(h). To step from the starting point at *t=0* to an estimate at *t=h*, follow the procedure below.

y'(0) &= - 2y(0)\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = 0 \cr

{k_1} &= - 2y(0)\quad &{\rm{derivative}}\;{\rm{at}}\;t = 0 \cr

{y_1}\left( {{h \over 2}} \right) &= y(0) + {k_1} {h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = h/2 \cr

{k_2} &= - 2{y_1}\left( {{h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = h/2 \cr

y(h) &= y(0) + y'(0)h + y''(0){{{h^2}} \over 2} + \cdots \quad &{\rm{Taylor}}\;{\rm{Series}}\;{\rm{around}}\;t = 0 \cr

y(t) &\approx y(0) + y'(0)h\quad &{\rm{Truncate}}\;{\rm{Taylor}}\;{\rm{Series}} \cr

{y^*}(h) &= y(0) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y(h) \cr} $$

In general, to go from the estimate *t=t₀* to an estimate at *t=t₀+h*

y'({t_0}) &= - 2y({t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{k_1} &= - 2{y^*}({t_0})\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + {h \over 2}} \right) &= y^*({t_0}) + {k_1} {h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{k_2} &= - 2{y_1}\left( {{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

y({t_0} + h) &= y({t_0}) + y'({t_0})h + y''(0){{{h^2}} \over 2} + \cdots \quad &{\rm{Taylor}}\;{\rm{Series}}\;{\rm{around}}\;t = {t_0} \cr

y({t_0} + h) &\approx y({t_0}) + y'({t_0})h\quad &{\rm{Truncated}}\;{\rm{Taylor}}\;{\rm{Series}} \cr

{y^*}({t_0} + h) &= y({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y(t_0+h) \cr} $$

We can use MATLAB to perform the calculation described above. A simple loop accomplishes this:

%% Example 1 % Solve y'(t)=-2y(t) with y0=3, midpoint method y0 = 3; % Initial Condition h=0.2; % Time step t = 0:h:2; % t goes from 0 to 2 seconds. yexact = 3*exp(-2*t) % Exact solution (in general we won't know this ystar = zeros(size(t)); % Preallocate array (good coding practice) ystar(1) = y0; % Initial condition gives solution at t=0. for i=1:(length(t)-1) k1 = -2*ystar(i) % Approx for y gives approx for deriv y1 = ystar(i)+k1*h/2; % Intermediate value k2 = -2*y1 % Approx deriv at intermediate value. ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y end plot(t,yexact,t,ystar); legend('Exact','Approximate');

The MATLAB commands match up easily with the steps of the algorithm. A slight variation of the code was used to show the effect of the size of *h*on the accuracy of the solution (see image below). Note that larger values of *h* result in poorer approximations, but that the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of *h*.

For a first order ordinary differential equation defined by

$${{dy(t)} \over {dt}} = f(y(t),t)$$to progress from a point at *t=t₀*, *y*(t₀)*, by one time step, *h*, follow these steps (repetitively).

{k_1} &= f({y^*}({t_0}),\;{t_0})\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + {h \over 2}} \right) &= {y^*}({t_0}) + {k_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + {h \over 2} \cr

{k_2} &= f\left( {{y_1}\left( {{t_0} + {h \over 2}} \right),\;{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + {h \over 2} \cr

{y^*}\left( {{t_0} + h} \right)& = {y^*}({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y\left( {{t_0} + h} \right) \cr} $$

Notes:

- an initial value of the function must be given to start the algorithm.
- see the MATLAB programs on this page for examples.
- this is often refered to as the "midpoint" algorithm for Second Order Runge-Kutta because it uses the slope at the midpoint,
*k*_{2.}

Adding an input function to the differential equation presents no real difficulty. You just need to ensure that you evaluate the function at *t=t₀* to find *k _{1}*, and at

So

$$\eqalign{y'({t_0}) &= - 2y({t_0}) + \cos (4{t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{k_1} &= - 2{y^*}({t_0}) + \cos (4{t_0})\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

&{y_1}\left( {{t_0} + {h \over 2}} \right) = {y^*}({t_0}) + {k_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{k_2} &= - 2{y_1}\left( {{t_0} + {h \over 2}} \right) + \cos \left( {{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{y^*}({t_0} + h) &= y({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y({t_0} + h) \cr} $$

Note: the exact solution in this case is $y(t) = 2.9{e^{ - 2{\kern 1pt} t}} + 0.1\cos (4t) + 0.2\sin (4t)$.

We can use MATLAB to perform the calculation described above. To perform this new approximation all that is necessary is to change the calculation of* k _{1}* and k

%% Example 2 % Solve y'(t)=-2y(t)-cos(4t) with y0=3 y0 = 3; % Initial Condition h=0.2; % Time step t = 0:h:2; % t goes from 0 to 2 seconds. yexact = 0.1*cos(4*t)+0.2*sin(4*t)+2.9*exp(-2*t); % Exact Solution ystar = zeros(size(t)); % Preallocate array (good coding practice) ystar(1) = y0; % Initial condition gives solution at t=0. for i=1:(length(t)-1) k1 = -2*ystar(i)+cos(4*t(i)); % Approx for y gives approx for deriv y1 = ystar(i)+k1*h/2; % Intermediate value k2 = -2*y1+cos(4*(t(i)+h/2)); % Approx deriv at intermediate value. ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y end plot(t,yexact,t,ystar); legend('Exact','Approximate');

A modified version of the program (that uses several values of *h*) generated the plot below. As before, the solution is better with smaller values of *h*, and the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of *h*.

Using a non-linear differential equation presents no real difficulty. Again, you just need to ensure that you evaluate the function at t=t₀ to find k_{1}, and at t=t₀+frac12;h to find k_{2}. Consider

with solution (described here)

$$y(t) = {e^{t - {t^2}}}$$The process follows those that came before

$$\eqalign{y'({t_0}) &= y({t_0})(1 - 2{t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{k_1} &= y^*({t_0})(1 - 2{t_0})\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + {h \over 2}} \right) &= {y^*}({t_0}) + {k_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{k_2} &= - 2{y_1}\left( {{t_0} + {h \over 2}} \right)\left( {1 - 2\left( {{t_0} + {h \over 2}} \right)} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{y^*}({t_0} + h) &= y({t_0}) + {k_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;y({t_0} + h) \cr} $$

As before, to perform this new approximation all that is necessary is to change the calculation of* k _{1}* and k

%% Example 3 % Solve y'(t)=y(t)(1-2t) with y0=1 y0 = 1; % Initial Condition h=0.2; % Time step t = 0:h:2; % t goes from 0 to 2 seconds. % Exact solution, hard to find in this case (in general we won't have it) yexact = exp(t-t.^2); ystar = zeros(size(t)); % Preallocate array (good coding practice) ystar(1) = y0; % Initial condition gives solution at t=0. for i=1:(length(t)-1) k1 = ystar(i)*(1-2*t(i)); % Approx for y gives approx for deriv y1 = ystar(i)+k1*h/2; % Intermediate value k2 = y1*(1-2*(t(i)+h/2)); % Approx deriv at intermediate value. ystar(i+1) = ystar(i) + k2*h; % Approximate solution at next value of y end plot(t,yexact,t,ystar); legend('Exact','Approximate');

A modified version of the program (that uses several values of *h)* generated the plot below. As before, the solution is better with smaller values of *h*, and the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of *h*.

If we start with a higher order differential equation

$$\displaylines{{{{d^3}y(t)} \over {dt}} + 4{{{d^2}y(t)} \over {dt}} + 6{{dy(t)} \over {dt}} + 4y(t) = \gamma (t)\quad \quad \quad \gamma (t) = {\rm{unit step function}} \cr

{\left. {{{{d^2}y(t)} \over {dt}}} \right|_{t = {0^ + }}} = 0,\quad \quad {\left. {{{dy(t)} \over {dt}}} \right|_{t = {0^ + }}} = - 1,\quad \quad y({0^ + }) = 0\quad \cr} $$

We can rewrite as a matrix differential equation (see previous example if this is unclear)

$$\eqalign{ {\bf{q}}'(t) &= \left[ {\matrix{ {{q_1}'(t)} \cr {{q_2}'(t)} \cr {{q_3}'(t)} \cr} } \right] = \left[ {\matrix{ 0 & 1 & 0 \cr 0 & 0 & 1 \cr{ - 4} & { - 6} & { - 4} \cr} } \right]\left[ {\matrix{ {{q_1}(t)} \cr {{q_2}(t)} \cr {{q_3}(t)} \cr} } \right] + \left[ {\matrix{ 0 \cr 0 \cr 1 \cr} } \right]\gamma (t) \cr

{\bf{q}}'(t) &= {\bf{Aq}}(t) + {\bf{B\gamma }}(t) \cr {\bf{A}} &= \left[ {\matrix{ 0 & 1 & 0 \cr 0 & 0 & 1 \cr

{ - 4} & { - 6} & { - 4} \cr} } \right],\quad \quad \quad {\bf{B}} = \left[ {\matrix{ 0 \cr 0 \cr 1 \cr } } \right] \cr} $$

Proceed as before (using matrices in place of scalars where appropriate

$$\eqalign{{\bf{q}}'({t_0}) &= {\bf{Aq}}({t_o}) + {\bf{B}}\gamma ({t_0})\quad &{\rm{expression}}\;{\rm{for}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{{\bf{k}}_1} &= {\bf{A}}{{\bf{q}}^*}({t_o}) + {\bf{B}}\gamma ({t_0})\quad &{\rm{approximate}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{{\bf{q}}_1}\left( {{t_0} + {h \over 2}} \right) &= {{\bf{q}}^*}({t_o}) + {{\bf{k}}_1}{h \over 2}\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{{\bf{k}}_2} &= {\bf{A}}{{\bf{q}}_1}\left( {{t_0} + {h \over 2}} \right) + {\bf{B}}\gamma \left( {{t_0} + {h \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h/2 \cr

{{\bf{q}}^*}\left( {{t_0} + h} \right) &= {{\bf{q}}^*}({t_o}) + {{\bf{k}}_2}h\quad &{\rm{estimate}}\;{\rm{of}}\;{\bf{q}}({t_0} + h) \cr} $$

Note: the exact solution to this problem is $y(t) = {1 \over 4}+{{\rm{e}}^{ - t}}{\mkern 1mu} \left( {\cos (t) - {5 \over 2}\sin (t)} \right) - {5 \over 4}{{\rm{e}}^{ - 2{\kern 1pt} t}}$, for t≥0.

To perform this new approximation all that is necessary is to change the approriate variables from scalars to vectors or matrices, and to define the **A** and **B** matrices. Note that in this case we multiply the **B** matrix by 1 since the input is a unit step (*γ(t)*=1 for *t*≥0). If the input were a sine wave, the sine wave would multiply **B**; if there is no input, it is not necessary to include the **B** matrix at all (it is multiplied by 0).

%% Example 4 % Solve y'''(t)+4y''(t)+6y'(t)+4y(t)=gamma(t), y''(0)=0, y'(0)=-1, y(0)=0 q0 = [0; -1; 0]; % Initial Condition (vector) h=0.2; % Time step t = 0:h:5; % t goes from 0 to 2 seconds. A = [0 1 0; 0 0 1; -4 -6 -4]; % A Matrix B = [0; 0; 1]; % B Matrix yexact = 1/4 + exp(-t).*(cos(t) - 5*sin(t)/2) - 5*exp(-2*t)/4; % Exact soln qstar = zeros(3,length(t)); % Preallocate array (good coding practice) qstar(:,1) = q0; % Initial condition gives solution at t=0. for i=1:(length(t)-1) k1 = A*qstar(:,i)+B*1; % Approx for y gives approx for deriv q1 = qstar(:,i)+k1*h/2; % Intermediate value k2 = A*q1+B*1; % Approx deriv at intermediate value. qstar(:,i+1) = qstar(:,i) + k2*h; % Approx solution at next value of q end plot(t,yexact,t,qstar(1,:)); % ystar = first row of qstar legend('Exact','Approximate');

A modified version of the program (that uses several values of *h*) generated the plot below. As expected, the solution is better with smaller values of *h*, and the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of *h*.

The Second Order Runge-Kutta algorithm described above was developed in a purely ad-hoc way. It seemed reasonable that using an estimate for the derivative at the midpoint of the interval between *t₀* and *t₀+h* (i.e., at *t₀+½h*) would result in a better approximation for the function at *t₀+h*, than would using the derivative at *t₀* (i.e., Euler's Method &emdash; the First Order Runge-Kutta). And, in fact, we showed that the resulting approximation was better. But, is there a way to derive the Second Order Runge-Kutta from first principles? If so, we might be able to develop even better algorithms.

In the following derivation we will use two math facts that are reviewed here. You should be familiar with this from a course in multivariable calculus.

**1) **If there is a function of two variable, *g(x,y)*, it's Taylor Series about *x₀*, *y₀* is given by:

&= g + {g_x}\Delta x + {g_y}\Delta y + \cdots \cr} $$

where *Δx* is the increment in the first dimenstion, and *Δy* is the increment in the second dimension. In the last line we use a shorthand notation that removes the explicit functional notation, and also represents the partial derivative of *g* with respect to *x* as *g _{x}* (and likewise for

**2) **If z is a function of two variables *z=g(x,y)*, where *x=r(t)* and *y=s(t)*, then by the chain rule for partial derivatives

In particular if

$$y'(t) = {{dy(t)} \over {dt}} = f(y(t),t) = f$$then

$$y''(t) = {{df(y(t),t)} \over {dt}} = {{\partial f} \over {\partial t}}{{dt} \over {dt}} + {{\partial f} \over {\partial y}}{{dy} \over {dt}} = {{\partial f} \over {\partial t}} + {{\partial f} \over {\partial y}}{{dy} \over {dt}} = {f_t} + {f_y}f$$To start, recall that we are expressing our differential equation as

$${{dy(t)} \over {dt}} = y'\left( t \right) = f(y(t),t)$$Now define two approximations to the derivative

$$\eqalign{& {k_1} = f(y^*({t_0}),{t_0}) \cr

& {k_2} = f(y^*({t_0}) + \beta h{k_1},{t_0} + \alpha h) \cr} $$

In all cases α and β will represent fractional values between 0 and 1. These equation state that *k _{1}* is the approximation to the derivative based on the estimated value of

To update our solution with the next estimate of *y(t)* at *t₀+h* we use the equation

This equation states that we get the value of *y*(t₀+h)* from the value of *y*(t₀)* plus the time step, *h*, multiplied by a slope that is a weighted sum of *k _{1}* and

Our goal now is to determine, from first principles, how to find the values *a*, *b*, *α* and *β* that result in low error. Starting with the update equation (above) we can substitue the previously given expressions for *k _{1}* and

We can use a two-dimensional Taylor Series (where the increment in the first dimension is *βhk _{1}*, and the increment in the second dimension is

f\left( {{y^*}({t_0}) + \beta h{k_1},\,{t_0} + \alpha h} \right) &= f + {f_y}\beta h{k_1} + {f_t}\alpha h + \cdots \cr

&= f + {f_y}f\beta h + {f_t}\alpha h + \cdots \cr} $$

In the last line we used the fact the *k _{1}=f* . The ellipsis denotes terms that are second order or higher. Substituting this into the previous expression for

{y^*}({t_0} + h) &= {y^*}({t_0}) + h\left( {af + b\left( {f + {f_y}f\beta h + {f_t}\alpha h + \cdots } \right)} \right) \cr

&= {y^*}({t_0}) + haf + hbf + {h^2}b\beta {f_y}f + {h^2}{f_t}b\alpha + \cdots \cr

&= {y^*}({t_0}) + h(a + b)f + {h^2}{f_t}b\alpha + {h^2}b\beta {f_y}f + \cdots \cr} $$

The ellipsis was multiplied by *h* between the first or second line and now represents terms that are third order or higher.

To finish we compare this approximation with the expressionfor a Taylor Expansion of the exact solution (going from the first line to the second we used the chain rule for partical derivatives). In this expression the ellipsis represents terms that are third order or higher.

$$\eqalign{y({t_0} + h) &= y({t_0}) + h{\left. {y'(t)} \right|_{{t_0}}} + {{{h^2}} \over 2}{\left. {y''(t)} \right|_{{t_0}}} + \cdots \cr

&= y({t_0}) + hf + {{{h^2}} \over 2}\left( {{f_t} + {f_y}f} \right) + \cdots \cr

&= y({t_0}) + hf + {{{h^2}} \over 2}{f_t} + {{{h^2}} \over 2}{f_y}f + \cdots \cr} $$

Comparing this expression with our final expression for the approximation,

$${y^*}({t_0} + h) = {y^*}({t_0}) + hf(a + b) + {h^2}{f_t}b\alpha + {h^2}b\beta {f_y}f + \cdots $$we see that they agree up to the error terms (third order and higher) represented by the ellipsis if we define the constants, *a*, *b*, *α* and *β* such that

a + b = 1 \cr

b\alpha = {1 \over 2} \cr

b\beta = {1 \over 2} \cr} $$

This system is underspecified, there are four unknowns, and only three equations, so more than one solution is possible. Clearly the choice we made* a=0* and *b=1* and *α=β=½* is one of these solutions. Because all of the terms of the approximation are equal to the terms in the exact solution, up to the error terms, the local error of this method is therefore *O(h ^{3})* (O(

The global error of the Second Order Runge-Kutta algorithm is *O(h ^{2})*.

Another common choice for the coefficients of the algorithm are *a=**b=½* and *α=β=1*. Before giving an example, let's figure out, intuitively what this is doing. We start with our equations for k_{1}, k_{2}, and y*(t₀+h).

{k_1} &= f({y^*}({t_0}),{t_0}) \cr

{k_2} &= f({y^*}({t_0}) + \beta h{k_1},{t_0} + \alpha h) \cr

{y^*}({t_0} + h) &= {y^*}({t_0}) + h(a{k_1} + b{k_2}) \cr} $$

Now put in our chosen constants

$$\eqalign{{k_1} &= f({y^*}({t_0}),{t_0}) \cr

{k_2} &= f({y^*}({t_0}) + h{k_1},{t_0} + \alpha h) \cr

{y^*}({t_0} + h) &= {y^*}({t_0}) + h\left( {{{{k_1} + {k_2}} \over 2}} \right) \cr} $$

In terms of our algorithm description

$$\eqalign{{k_1} &= f({y^*}({t_0}),\;{t_0})\quad& {\rm{estimate}}\;{\rm{of}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {k_1}h\quad &{\rm{intermediate}}\;{\rm{estimate}}\;{\rm{of}}\;{\rm{function}}\;{\rm{at}}\;t = {t_0} + h \cr

{k_2} &= f\left( {{y_1}\left( {{t_0} + h} \right),\;{t_0} + h} \right)\quad& {\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;t = {t_0} + h \cr

{y^*}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + h\left( {{{{k_1} + {k_2}} \over 2}} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{y}}\left( {{{\rm{t}}_{\rm{0}}}{\rm{ + h}}} \right)\;{\rm{using}}\;{\rm{average}}\;{\rm{of}}\;{\rm{slopes}} \cr} $$

We will again use the the differential equation

$${{dy(t)} \over {dt}} + 2y(t) = 0\quad {\rm{or }} \quad {{dy(t)} \over {dt}} = - 2y(t)$$with the initial condition set as *y(0)=3* with exact solution *y(t)=3e ^{-2t}*,

The slope at the end point is shown below (and is labeled* k _{2}*); numerically

We now take the average of these two slopes (-5.0110) and use it to move forward by one time step.

$${y^*}\left( h \right) = {y^*}(0) + h\left( {{{{k_1} + {k_2}} \over 2}} \right) = 3 + 0.2\left( { - 5.0110} \right) = 1.9978$$

This is very close to the exact answer, y(0.2)=2.0110, with an error of about 0.6%. This is close to that of the midpoint method (it is slightly worse, but that is becaus of the specific problem chosen, that won't generally be the case), and much better than that of the first order algorithm.

For a first order ordinary differential equation defined by

$${{dy(t)} \over {dt}} = f(y(t),t)$$to progress from a point at *t=t₀*, *y*(t₀)*, by one time step, *h*, follow these steps (repetitively).

{k_1} &= f({y^*}({t_0}),\;{t_0})\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{derivative}}\;{\rm{at}}\;t = {t_0} \cr

{y_1}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + {k_1}h\quad &{\rm{approximate}}\;{\rm{enpoint}}\;{\rm{value}}\;{\rm{at}}\;t = {t_0} + h \cr

{k_2} &= f\left( {{y_1}\left( {{t_0} + h} \right),\;{t_0} + h} \right)\quad &{\rm{estimate}}\;{\rm{of}}\;{\rm{slope}}\;{\rm{at}}\;{\rm{enpoint,}}\;t = {t_0} + h \cr

{y^*}\left( {{t_0} + h} \right) &= {y^*}({t_0}) + h{{\left( {{k_1} + {k_2}} \right)} \over 2}\quad &{\rm{estimate}}\;{\rm{of}}\;y\left( {{t_0} + h} \right)\;{\rm{using}}\,{\rm{average}}\;{\rm{slope}} \cr} $$

The following MATLAB code repeats Example 1 (a linear differential equation with no input). Example 1 used the "midpoint" method, this example uses the "endpoint" method. The MATLAB commands match up easily with the steps of the algorithm (only the lines that calculate y1 and k2 have changed from the midpoint method).

%% Example 5 % Solve y'(t)=-2y(t) with y0=3, endpoint method y0 = 3; % Initial Condition h=0.2; % Time step t = 0:h:2; % t goes from 0 to 2 seconds. yexact = 3*exp(-2*t) % Exact solution (in general we won't know this ystar = zeros(size(t)); % Preallocate array (good coding practice) ystar(1) = y0; % Initial condition gives solution at t=0. for i=1:(length(t)-1) k1 = -2*ystar(i) % Approx for y gives approx for deriv y1 = ystar(i)+k1*h; % enpoint approximation k2 = -2*y1 % Approx deriv at endpoint. ystar(i+1) = ystar(i) + h*(k1+k2)/2; % Approx solution at next value of y end plot(t,yexact,t,ystar); legend('Exact','Approximate');

A slight variation of the code was used to show the effect of the size of *h*on the accuracy of the solution (see image below). Note that larger values of *h* result in poorer approximations, but that the solutions are much better than those obtained with the First Order Runge-Kutta for the same value of *h* (and appears similar to results obtained with the mid-point method).

The other examples can be coded for the endpoint method in a similar way.

With a little more work we can develop an algorithm that is accurate to even higher order. The next page describes just such a method.

© Copyright 2005 to 2015 Erik Cheever This page may be freely used for educational purposes.

Erik Cheever Department of Engineering Swarthmore College