Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Lecture on Numerical Analysis
Dr.-Ing. Michael Dumbser
M. Dumbser
1/9
17 / 12 / 2007
Analisi Numerica
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Ordinary Differential Equations
An equation of the form
dy
 f ( y (t ), t )
dt
is called a nonlinear ordinary differential equation (ODE) of first order. The simplest
first order ODE is the equation
dy
 ay
dt
y  y (t )
with
The exact solution of this simple ODE is obtained as follows:
dy
 adt
y
y
dy
y y  t adt
0
0
ln y  ln y0  a t  t0 
M. Dumbser
2/9
t
ln y yy
0
y
ln
 a t  t0 
y0
 at t0
t
y (t )  y0ea t t0 
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Ordinary Differential Equations
A simple nonlinear ordinary differential equation (ODE) of first order is given by
dy
 ay 2
dt
with
y  y (t )
The exact solution of this simple nonlinear ODE can be computed in a similar way:
dy
 adt
2
y
y
dy
y y 2  t adt
0
0
1 1
  a t  t0 
y0 y
M. Dumbser
3/9
t



y
1
t



at
t0
y  y0
 1

y (t )    a t  t0 
 y0

1
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Numerical Methods for the Integration of ODE
The simplest numerical method for the integration of first order ODEs is the
so-called explicit (forward) Euler method, invented already in the 18th century by the
mathematician Leonhard Euler. It corresponds to the use of a first order accurate
finite difference scheme to approximate the derivative on the left hand side of the ODE
and to compute the operator on the right hand side at the current state yn and time tn:
dy
 f ( y, t )
dt
y n 1  y n
 f ( yn , tn )
t
y
n 1
 y  t  f ( y , t )
n
n
n
t n 1  t n  t
Computing the function f on the right hand side at the unknown time tn+1, we
obtain the so-called implicit (backward) Euler method:
y n1  y n  t  f ( y n1, t n1 )
M. Dumbser
4/9
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Numerical Methods for the Integration of ODE
A simple modification of the forward Euler scheme, which yields a second order
accurate method, is based on the following approach: Integrate the ODE from time
level tn to time level tn+1:
t n 1

tn
dy

dt
t n 1

y n 1  y n 
f ( y , t )dt
t n 1
tn
 f ( y, t )dt
tn
The integral on the right hand side is now evaluated using a one-point Gaussian
quadrature formula (the midpoint is the integration point):
y
n 1
 y  t  f ( y
n
n  12
,t
n  12
)
Unfortunately, the value of y is not known at tn+1/2, but we can approximate it making a
half time-step of the forward Euler scheme. This means, that first, a so-called predictor
is computed at the half time-step t/2:
y
M. Dumbser
5/9
The corrector is then:
y
n  12
n 1
t
 y   f ( yn , tn )
2
n
 y  t  f ( y
n
n  12
,t
n  12
)
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Numerical Methods for the Integration of ODE
One of the most famous and mostly commonly used schemes for ODE integration
is the explicit fourth-order scheme of Runge and Kutta:
k1  f ( y n , t n )
k2  f ( y n 
t
2
k1 , t n  2t )
k3  f ( y n 
t
2
k2 , t n  2t )
k4  f ( y n  tk3 , t n  t )
y
M. Dumbser
6/9
n 1
t
 y  k1  2k2  2k3  k4 
6
n
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 1
Write a MATLAB script ODECompare.m that integrates the following two ODE:
dy
 y
dt
y ( 0)  1
t  [0;10]
(1)
Use the explicit (forward) Euler method, using different time steps and compare
with the exact solutions.
(2)
Use the modified second order Euler method (=second order Runge-Kutta
scheme) with different time steps, and compare with the exact solutions.
(3)
Use the fourth order Runge-Kutta scheme with various time steps. Compare with
the exact solutions.
M. Dumbser
7/9
dy
  y2
dt
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 2
Wite a MATLAB script Springs.m that integrates the linear ODE system for position
and velocity for a system of N point masses, connected by springs as follows:
m1
mi-1
mi
mi+1
mN
The ODE system for the positions of the point masses is simply:
dxi
 vi
dt
The forces acting on point mass i are
Fi    xi 1  xi  Li 1 ki 1
Fi    xi 1  xi  Li ki
F1  0
FN  0
Where ki are the spring stiffnesses and Li are the lengths of the springs in their
undeformed state. The ODE system for the velocities of the point masses is then
dvi
1 
 Fi  Fi  
dt mi
M. Dumbser
8/9
Use the second order Runge-Kutta scheme (modified Euler method).
Università degli Studi di Trento
Dipartimento d‘Ingegneria Civile ed Ambientale
Dr.-Ing. Michael Dumbser
Analisi Numerica
Exercise 3
Wite a MATLAB script Orbit.m that integrates the linear ODE system for position
and velocity for a planet, turning around a star.
The ODE system for the position vector of the planet is:

dx 
v
dt
The gravitation force of the star, acting on the planet, is given by:
 

xS  x P
mS mP
FG      
xS  x P xS  x P
2
The ODE system for the velocity vector of the planet is

dv P
1 

FG
dt
mP
M. Dumbser
9/9
Use the second order Runge-Kutta time integration scheme (modified Euler method).