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 / yexact0 / 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 / yexact0 / 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.