Introduction to Scientific
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
Chapter 9 The Initial Value
Problem
? Basic concepts
? The Runge-Kutta Methods
? The Adams Methods
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.40
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Solutions to y'(t) = -5 y(t)
tty 5)(' ??
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.40
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Five Steps of Euler Method (y'=-5y,y(0)=1)
),(1 nnnnn ytfhyy ???
Euler Method
?
?
?
?
?
???
?
?
?
?
?
??
??
???
??
?
1
0
N0
1
0
10N00
||||
:1f o r t h e n i n t e r s e c t,
:0}:))(,{(
i e s t r a j e c t o r t h eof n o n e a n d ]t,[ a l lf o r
0
),(
If
.)(
bye r r o r n t r u n c a t i ol o c a l t h ea n d
)(b y e r r o r g l o b a l t h eD e f in e
, a n dg i v e n a r e )y,(),.,,,y,( w h e r e
)()),(
I V P t h es o l v e s
t h a te x i s t s ( t )f u n c t io n a,:0f o r t h a t,A s s u m e 9 T h e o r e m
k
kn
n
y
nnnn
nnn
NN
nn
n
L T Eg
Nn
nnNtttyt
tt
y
ytf
f
ytyL T E
ytyg
ttttt
ytytf ( t,yy ' ( t )
yNn
?
-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.40
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
y' = -5.0y,y(0) = 1
* = exact solution
o = computed solution
Euler Solution (h=0.1) of y'=-5y,y(0) = 1
Showtrunc.m
Nt(th /)0m a x ??
2||
2
2
hML T E
n ?
hMtthNML T Eyty
N
n
nN 2
1
0ma x
2
2ma x 22|||)(| ?
?
?????
to lMN tt ?? 2
2
0m a x
2
)(
n = ceil(((tmax-t0)^2*M2)/(2*tol))+1;
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50
2
4
6
8 x 10
-4 Fixed h Euler Error for y'=-y,y(0) = 1
tol = 0.010,n = 1251
Showfixedeuler.m
tol=.01;
[tvals,yvals] =
FixedEuler('f1',1,0,5,1,tol);
It is better to determine h adaptively so that longer step sizes can
be taken in regions where the solution is smooth
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
Euler in 3-Digit Arithmetic
h = 1/140
h = 1/160
h = 1/180
Note that the error gets worse as h get smaller
because the step size are in the neighborhood of
unit roundoff.
Stability
)(10)(' tyty ??
nn yhy )101(1 ???
1|101| ?? h 5/1?h
Backward Euler Method
),( 111 ??? ?? nnnnn ytfhyy
n
n
nnnn yhyhyy ?? ???? ?? 1
1
11
yy ??'
Implicit method
?
?
?
?
?
?
?
?
?
?
?
)(
)(
)(
1
tz
tz
ty
d
?system
The Runge-Kutta Methods
211
12
1
),(
),(
bkakyy
kyhthfk
ythfk
nn
nn
nn
???
???
?
?
??
12
12
1
?
?
??
?
?
b
b
ba )()( 311 hOyty nn ?? ??
2/)(
),(
),(
211
12
1
kkyy
kyhthfk
ythfk
nn
nn
nn
???
???
?
?
RK2:
? ?
)22(
6
1
,
2
1
,
2
2
1
,
2
),(
43211
34
23
12
1
kkkkyy
kyhthfk
ky
h
thfk
ky
h
thfk
ythfk
nn
nn
nn
nn
nn
?????
???
?
?
?
?
?
?
???
?
?
?
?
?
?
???
?
?
Fourth order Runge-Kutta Methods
0 100 200 300 400 500 600
10-14
10-12
10-10
10-8
10-6
10-4
10-2
100
Runge-Kutta on y'(t) = -y(t),y(0) = 1,0<=t<=5,h = 5/n
n (Number of Steps)
Ma
xim
um
ab
so
lut
e
err
or
k=1
k=2
k=3
k=4
k=5
Do not conclude from the example that high-order methods are
necessarily more accurate,If the higher derivatives of the
solution are badly behaved,then it may well be the case that a
lower-order method gives more accurate results,One must also
be mindful of the number of f-evaluations that are required to
purchase a given level of accuracy,The situation is analogous to
what we found in the quadrature unit,Of course,the best
situation is for the IVP software to handle the selection of method
and step.
The Matlab IVP Solvers ode23 and ode45
ode23 is based on a pair of second- and third-order Runge-
Kutta methods,Ode45 based on fourth- and fifth-order RK,
2)0(,0)0(,
))()(((
)(
)(
0)0(,4)0(,
))()(((
)(
)(
2/322
2/322
??
?
??
??
?
??
yy
tytx
ty
ty
xx
tytx
tx
tx
???
???
Two-body problem
2)0())()(/()()(
0)0()()(
0)0())()(/()()(
4)0()()(
4
2/32
3
2
134
343
2
2/32
3
2
112
121
????
??
????
??
ututututu
ututu
ututututu
ututu
?
?
?
?
Kepler
-1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Kepler Problem,ode23 with Default Tolerances
Number of Steps = 53
Ode23 for Kepler
-1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
Kepler Problem,ode23 with Specified Output Times
Spline Fit
ode23 Output
Ode23 and spline
-1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Kepler Problem,ode23 with RelTol = 10-6 and AbsTol = 10-8
Number of Steps = 517
Ode23
-1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Kepler Problem,ode45 with RelTol = 10-6 and AbsTol = 10-8
Number of Steps = 208
Ode45
The Adams method
?? ?? ????? 11 )()()(')()( 1 nnnn ttnttnn dttftydttytyty
? ? ?? ?? 1 )(11 nntt knn dttpyy
Kth order Adams-Bashforth(AB) method
One-point Newton-cotes ruler
The first-order method is the Euler method.
),(1 nnnnn ytfhyy ???
)()( 1
1
1
11 ?
?
?
?? ?
???
n
n
nn
nk tth
ffftp
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
??
?
?
?
?
???
??
1
11
1
1
1
1
1
2
2
)())(,(
11
n
n
n
n
n
nnn
t
t
n
n
nn
n
t
t
f
h
h
f
h
hhh
dttt
h
ff
fdttytf
n
n
n
n
)3(
2 11 ??
??? nnnn ffhyyhhh nn ?? ? 1
Second order Adams method
)3(
2 11 ??
??? nnnn ffhyy
)51623(
12 211 ???
???? nnnnn fffhyy
nnn hfyy ??? 1
)9375955(
24 3211 ????
????? nnnnnn ffffhyy
)2511274261627741901(720 43211 ????? ?????? nnnnnnn fffffhyy
? ? ??? ??? 1 ))())(,(()( 111 nntt knnn dttptytfyty
Implement
Multistep method
How do we perform the first step when there is no,history”?
For a kth order Adams framework we could use a kth
order RK method.
0 100 200 300 400 500 60010
-12
10-10
10-8
10-6
10-4
10-2
Adams-Bashforth on y'(t) = -y(t),y(0) = 1,0<=t<=5,h = 5/n
n (Number of Steps)
Ma
xim
um
ab
so
lut
e e
rro
r
k=1
k=2
k=3
k=4
k=5
The Adams-Moulton method
? ? ?? ?? 1 )(11 nntt knn dttpyy
The kth-order Adams-Moulton method uses a degree k-1
interpolant of the points (tn+1-j,fn+1-j),j=0:k-1
The first-order method is the backward Euler method.
)),((
2 111 nnn
n
nn fytf
hyy ???
???
Second-order method
),( 111 ??? ?? nnnn ythfyy
)519),(9(
24 21111 ?????
????? nnnnnnn fffytfhyy
)),((
2 111 nnn
n
nn fytf
hyy ???
???
)8),(5(
12 1111 ????
???? nnnnnnn ffytfhyy
)191 0 62 6 46 4 6),(2 5 1(7 2 0 321111 ?????? ?????? nnnnnnnn ffffytfhyy
The Predictor-Corrector Idea
To couple an Adams-Bashforth method with an Adams-
Moulton method of the same order,The idea is to predict
yn+1 using an Adams-Bashforth method and then to correct
its value using the corresponding Adams-Moulton method.
)(
2 1
)(
1 ?? ??? nnn
P
n ff
hyy
)),((
2
)(
11
)(
1 n
P
nnn
C
n fytf
hyy ???
???
0 100 200 300 400 500 60010
-14
10-12
10-10
10-8
10-6
10-4
10-2
Adams PC on y'(t) = -y(t),y(0) = 1,0<=t<=5,h = 5/n
n (Number of Steps)
Ma
xim
um
ab
so
lut
e
err
or
k=1
k=2
k=3
k=4
k=5
Stepsize Control
The idea behind error estimation in adaptive quadrature is to
compute the integral in question in two ways and then accept
or reject the better estimate based on the observed
discrepancies,The predictor-corrector framework presents us
with a similar opportunity,The Quality of y(c) can be estimated
from |y (c) - y(p) |,If the error is too large,we can reduce the
step,If the error it is too small,then we can length the step,
Properly handled,we can use this mechanism to integrate the
IVP across the interval of interest with steps that are as long as
possible given a prescribed error tolerance,In this way we can
compute the required solution,more or less minimizing the
number of f evaluations,The Matlab ode23 and ode45 are RK
based and do just that,
s t e pr r e c t p r e d i c t / c oa n o t h e r t r y a n d H a l v e
)),((
2
:R e p e a t
:f o r m u l a A M 2 t h eof
n a p p l i c a t i o r e p e a t e d t h r o u g h R e f i n e
:esp o s i b i l i t i t oa r e t h e r en o t,If
i o n t oa p p r o x i m a to u r as a c c e p t i n g w o r t h a n d
g o o df a i r l y p r o b a b l y is t h e n s m a l l,is || If
111
11
1
)(
1
)(
1
h
fytf
h
yy
yy
y
yy
nnnnn
( C )
nn
( C )
n
C
n
P
n
?
???
?
?
???
???
??
?
??
],[)(
12
1
)(
],[)(
12
5
)(
22
)3(3)(
11
11
)3(3)(
11
httyhyty
httyhyty
nn
C
nn
nn
P
nn
????
????
??
??
??
??
],[)(
2
1|| )3(3)(
1
)(
1 httyhyy nn
P
n
C
n ???? ?? ??
||
6
1|)(| )(
1
)(
11
)(
1
P
n
C
nn
C
n yytyy ???? ???
hyy
hyy
h
yy
fytf
h
yy
ff
h
yy
( C )
nn
( C )
nn
P
n
C
n
n
P
nnn
C
n
nnn
P
n
i n c r e a s e a n d s e t
t h e ns m a l l t o ois if e l s e
k e e p a n d s e t
r i g h t t h e na b o u t is if e l s e
a g a i n t r y a n d r e d u c e
t h e n b i g t o o is If
||
6
1
)),((
2
)(
2
11
11
)(
1
)(
1
)(
11
)(
1
1
)(
1
??
??
??
???
??
?
?
??
???
???
?
?
?
?
??
?
?
m ax
1
L T E
n
n
n 0m a x
|| tt hLT E n ?? ?
??
??
?
?
?
m a xm a x
1 0m a x1
L T E
n
n
n
n
n
n tt
h ??
0m a x10
1
tt
h
?
? ??