1

Numerical Solution of Ordinary Differential Equations

Consider the 1st order ordinary differential equation (ODE)

.

The initial condition can be taken as

.

Then we could use a Taylor series about and obtain the complete solution, or

Since

and

then we can find the first two terms. For the second derivative

.

Similarly the third derivative is

.

If we truncate at the third derivative

,

and

, where .

Euler’s Method

Take the Taylor series to 1st order, and let the interval , then

.

The error for a time step (the local error) is . The global error, after many steps, is . Then

where ,

where ,

where .

Example:

The exact solution can be found from

.

Let where , or . Then for all x ,

or , and . Since the right hand side is linear in x try . Then

and becomes which must hold for all x. Hence

, and B=-1, making , and since then

.

But @ , or , and . Making the complete solution

.

Using Euler’s method and taking

, since . In general,

n / xn / yn / yexact
0 / 0.00 / 1.0000 / 1.0000
1 / 0.02 / 1.0200 / 1.0204
2 / 0.04 / 1.0408 / 1.0416
3 / 0.06 / 1.0624 / 1.0637
4 / 0.08 / 1.0848 / 1.0866
5 / 0.10 / 1.1081 / 1.1103

For the error, , , can be defined as

The results plot as

It would be better to use the slope at the beginning and end of the increment (e.g., the average at each end), and although we don’t know the slope at the end we can approximate it.

Modified Euler Method

Let . Then an approximation for y at the end of the increment is

and an estimate for the slope at the end of the increment is .

We can now set

.

The error can be found from

and since

or

.

Hence the local error is and the global error is . Another way to write our results is

The previous example now can include modified Euler

n / xn / yeuler / ymodified / yexact
0 / 0.00 / 1.0000 / 1.0000 / 1.0000
1 / 0.02 / 1.0200 / 1.0204 / 1.0204
2 / 0.04 / 1.0408 / 1.0416 / 1.0416
3 / 0.06 / 1.0624 / 1.0637 / 1.0637
4 / 0.08 / 1.0848 / 1.0866 / 1.0866
5 / 0.10 / 1.1081 / 1.1104 / 1.1103

which is much better.

Runge-Kutta Methods

The modified Euler method is actually a two step (second order) Runge-Kutta algorithm.

These methods can be readily extended to eight and even tenth order. The derivations follow the same procedure. Assume for the second order method

where

, and

.

The parameters a, b,  and  are found by comparing to a Taylor series expansion.

Recall

.

But

.

Since ,

Using ,

,

or

Comparing to the Taylor series

.

Then which is three equations in four unknowns. Hence we can pick for the fourth equation any equation that is convenient.

For example we can take

or

.

Modified Euler is

.

Fourth Order Runge-Kutta

For fourth order Runge-Kutta the estimates for the changes are

and the updated value for y is found from

.

Note that all of the algorithms presented are for first order equations with only one dependent variable. These can be readily extended to systems of higher order differential equations. For those cases please see page 7.

Higher Order Differential Equations

Consider

and let the vector be defined as

where

which is now a vector first order equation, and the first order rules can be applied to a system of first order equations.

Systems of First Order Equations

Let

Euler:

Modified Euler:

where

Runge-Kutta (4th order):

where

Problem

The equation for a pendulum with a mass, m, at the end of a rod of negligible mass with length, L, is

.

For the initial conditions take

, and ,(1)

where

is the angle from the vertical, , and

Let . The equation becomes

. (2)

Now let then since

,

The equation can be written

,

and integrating

.

This is the conservation of energy. The initial condition is ,

hence or

(3)

is also the governing differential equation.

Integrate the equation of motion (2) subject to the initial conditions (1), and using Euler, modified Euler, and Runge-Kutta. Use the expression for the energy to check the accuracy of your integration. Integrate for one-quarter of a cycle and then determine the period for a complete cycle.