Chapter 14
Initial Value Problem
(Numerical Solution of Ordinary Differential
Equations)
1
Outline
Introduction
Euler’s method
Huen’s and Modified Euler’s methods
Runge-Kutta methods
Matlab tools
Example problems
2
Introduction
Differential equations arise in many science and
engineering applications.
The order of the differential equation is
determined by the highest-order derivative.
They can be linear or non-linear.
Linear 1st-order dv(t ) b h
differential equation = g − v(t )
dt m
2
d θ (t ) g
+ sin θ = 0 Non-linear 2nd-order
dt l differential equation
3
Euler’s method
We wish to solve ordinary differential equations of
the form dx
= f (t , x) x′(t )
dt
Recall Taylor’s formula, which is used throughout
this subject
x′′(t ) ∆t 2 x′′′(t ) ∆t 3
x(t + ∆t ) = x(t ) + x′(t ) ∆t + + +
2! 3!
and ignoring derivative terms greater than 1,
x(t + ∆t ) = x(t ) + x′(t ) ∆t
4
Euler’s method
Substitution gives a way to iteratively compute an
approximate solution x(t) given our differential
equation
dx x′(t )
= f (t , x)
dt
Taylor’s series x(t + ∆t ) = x(t ) + x′(t ) ∆t
Euler’s
method
x(t + ∆t ) = x(t ) + f (t , x) ⋅ ∆t
5
Euler’s method
Euler’s
method
x(t + ∆t ) = x(t ) + f (t , x) ⋅ ∆t
Euler’s method basically says that the new value
is approximately
New value = old value + slope X step size
x Generally, as the step
size decreases, our
solution improves
x(t + ∆t )
x(t ) slope = f (t , x)
t
∆t 6
Euler’s method
Recall the ball bearing problem. h
dv(t ) b
= g − v(t )
dt m
We implement the expression in a for-loop
old value
b
v(t + ∆t ) = v(t ) + g − v(t ) ⋅ ∆t
m
See Matlab example slope
7
Error
Numerical solution of ODEs involves two
types of error
Truncation or discretization error caused by
the approximations used to compute x(t)
Round-off error caused by the limited number
of significant digits that can be retained on a
computer
For Euler’s method truncation error is of the
order ∆t2
8
Modifications and improvements
to Euler’s method
Euler’s method is easily implemented and
can give satisfactory results.
However, to get better results we have to
decrease the step size, increasing the
required number of computations.
Modifications of Euler’s method can give
better results for a given step size.
Huen’s method
Modified Euler method
9
Huen’s method
Euler’s method assumes that the slope at
the start of the interval applies across the
entire interval.
Euler’s method assumes the same slope
across the entire interval
x
x(t + ∆t )
x(t ) slope = f (t , x)
t
∆t 10
Huen’s method
Huen’s method introduces another slope at
the end of the interval.
This new slope is averaged with the old one
to improve the result.
new slope = f (t + ∆t , x(t + ∆t ))
x
x(t + ∆t )
x(t ) old slope = f (t , x)
t
∆t 11
Huen’s method
Huen’s method has two steps.
Predictor Really just Euler’s
method here
0
x (t + ∆t ) = x(t ) + f (t , x) ⋅ ∆t
Corrector
f (t , x) + f (t + ∆t , x 0 (t + ∆t ))
x(t + ∆t ) = x(t ) + ⋅ ∆t
2
Average of the slope at the beginning
and end of the interval
12
Modified Euler’s method
Uses Euler’s method to compute a solution at the
midpoint of the interval x(t+∆t/2).
∆t ∆t
x(t + 2 ) = x(t ) + f (t , x) ⋅
2
Next this value is used to estimate a slope at the
midpoint.
′ ∆t ∆t ∆t
x (t + 2 ) = f (t + 2 , x(t + 2 ))
Finally, the resulting slope is used to extrapolate
to t+∆t
x(t + ∆t ) = x(t ) + x′(t + ∆2t ) ⋅ ∆t
13
Modified Euler’s method
First calculate x(t+∆
∆t/2)
then, x’(t+ ∆t/2)
x
and lastly, x(t+ ∆t)
x(t + ∆t )
x(t )
t
∆t
2
∆t
14
Comparison of Huen’s and
Modified Euler’s methods
Ball Bearing in Viscous Fluid
4.5
True
Huens
4
Modified Euler
h 3.5
Velocity m/sec
2.5
1.5
1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time (seconds)
15
Runge-Kutta methods
RK methods offer the advantage of higher-
order terms of the Taylor series, but
without having to compute these terms
explicitly.
Increment function
Has the generalized form
x(t + ∆t ) = x(t ) + φ (t , x(t ), ∆t ) ⋅ ∆t
Increment function is a representative
slope over the interval ∆t
16
Runge-Kutta methods
The increment function φ has the form
φ = a1k1 + a2 k2 + + an kn
where the a terms are constants and the k terms
are
Note use of k terms in
k1 = f ( t , x(t ) ) subsequent formulas.
k2 = f ( t + p1∆t , x(t ) + q11k1∆t )
k3 = f ( t + p2 ∆t , x(t ) + q21k1∆t + q22 k2 ∆t )
kn = f ( t + pn −1∆t , x(t ) + qn −1,1k1∆t + qn −1,2 k2 ∆t + + qn −1,n −1kn −1∆t )
17
Runge-Kutta
To find the a, p and q terms of the RK
method, we compare a RK formula to terms
in the Taylor series.
Algebraically very intensive.
It can be shown that for n=2, and with
proper selection of terms, the RK method is
identical to Huen’s and modified Euler
methods.
This is called a 2nd-order RK method.
18
Runge-Kutta
A 3rd-order RK method is shown here
1
x(t + ∆t ) = x(t ) + ( k1 + 4k2 + k3 ) ⋅ ∆t
6
with
k1 = f ( t , x(t ) )
∆t ∆t
k2 = f t + , x(t ) + k1
2 2
k3 = f ( t + ∆t , x(t ) − ∆t k1 + 2∆t k2 )
19
Runge-Kutta
Ball bearing problem…
Ball Bearing in Viscous Fluid
4.5
True
4 RK3
Euler
3.5
RK3 Huens
∆t = 0.2 sec
Velocity m/sec
2.5
1.5
0.5
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time (seconds)
20
Matlab tools for ODEs
21
Ball Bearing
problem once again h
Matlab recommends using ode45 as a first
try for your solution.
4.5
First write a function handle that
4
describes the ODE
3.5
function dv = ballBearing(t,v)
Velocity (m/sec)
3
g=9.81; b=1; m=0.15; 2.5
dv = g - (b/m)*v; 2
1.5
Next use one of the ODE solvers 0 0.2 0.4 0.6 0.8 1
Time (secs)
1.2 1.4 1.6 1.8 2
[tOut vOut]=ode45(@ballBearing, [0 2],[4.5]);
plot(tOut, vOut)
22
Non-linear ODE
Consider the simple pendulum with non-
linear 2nd –order ODE
2
d θ (t ) g
+ sin θ = 0
dt l
To use Matlab solvers, we must
write our 2nd-order ODE as
a system of 1st-order ODEs.
23
Non-linear ODE
To write the 2nd-order ODE as a system of
1st-order ODEs, introduce state variables
y = θ= y
y1 = θ and y2 = θ Also have
that
1 2
y = θ
2
We then write a system of 1st-order ODEs
as y1 = y2
d 2θ (t ) g
+ sin θ = 0 g
dt l y 2 = − sin y1
l
24
Non-linear ODE
Then we use one of our Matlab ODE solvers
Pendulum angle vs Time
First write a function handle 0.8
Θ
0.6 Θdot
function dydt = pendulum(t,y)
0.4
g=9.81; len = 1;
Pendulum angle (rad/sec)
0.2
dydt = [y(2);
-g/len*sin(y(1))]; 0
-0.2
-0.4
-0.6
Next use a solver that calls the
-0.8
function handle we created 0 1 2 3 4 5
Time (secs)
6 7 8 9 10
[tOut vOut]=ode45(@pendulum, [0 10],[0.2;0]);
Time interval Initial conditions
[y1(0), y2(0)] 25
Summary
Understand that the numerical methods used for solving
ODEs are based on Taylor’s formula.
Truncation and round-off error always present in the
numerical solution.
Higher-order methods give better performance at the
expense of greater computational effort.
A number of Matlab functions available.
Need to use a function handle script for your differential
equation.
Higher order ODEs need to be cast into a system of 1st-
order ODEs.
26