Transcript Numerical ODE - Philip H. Dybvig Home Page

Fin500J Mathematical Foundations in Finance
Topic 7: Numerical Methods for Solving
Ordinary Differential Equations
Philip H. Dybvig
Reference: Numerical Methods for Engineers, Chapra and Canale, chapter 25,
2006
Slides designed by Yajun Wang
Fin500J Topic 7
Fall 2010 Olin Business School
1
Numerical Solutions
 Numerical method are used to obtain a graph or a
table of the unknown function
 Most of the Numerical methods used to solve ODE are
based directly (or indirectly) on truncated Taylor series
expansion
 Taylor Series methods
 Runge-Kutta methods
Fin500J Topic 7
Fall 2010 Olin Business School
2
Taylor Series Method
The problem to be solved is a first order ODE
dy ( x)
 f ( x, y ),
dx
y ( x0 )  y0
Estimates of the solution at different base points
y( x0  h), y( x0  2h), y( x0  3h), ....
are computed using truncated Taylor series
expansions
Fin500J Topic 7
Fall 2010 Olin Business School
3
Taylor Series Expansion
T runcatedT aylorSeries Expansion
dy
h2 d 2 y
y ( x0  h)  y ( x0 )  h

dx xyxy0 , 2! dx2
0
hn d n y
 ... 
x  x0 ,
n! dxn
y  y0
x  x0 ,
y  y0
nth order Taylor series method uses nth
order Truncated Taylor series expansion
Fin500J Topic 7
Fall 2010 Olin Business School
4
First Order Taylor Series Method
(Euler Method)
dy
y ( x0  h)  y ( x0 )  h
 o( h 2 )
dx xyxy0 ,
0
Notation:
xn  x0  nh,
yn  y ( xn ),
dy
 f ( xi , yi )
dx xyxyi ,
i
Euler Method
yi 1  yi  h f ( xi , yi )
Fin500J Topic 7
Fall 2010 Olin Business School
5
Euler Method
P roblem:
Given thefirst order ODE
with theinitialcondit ion
y ( x)  f ( x, y )
y0  y ( x0 )
yi  y ( x0  ih)
Determine
for i  1,2,...
Euler Met hod:
y0  y ( x0 )
yi 1  yi  h f ( xi , yi )
Fin500J Topic 7
for i  1,2,...
Fall 2010 Olin Business School
6
Interpretation of Euler Method
y2
y1
y0
x0
Fin500J Topic 7
x1
Fall 2010 Olin Business School
x2
x
7
Interpretation of Euler Method
y1
Slope=f(x0,y0)
hf(x0,y0)
y0
x0
Fin500J Topic 7
y1=y0+hf(x0,y0)
h
x1
Fall 2010 Olin Business School
x2
x
8
Interpretation of Euler Method
y2=y1+hf(x1,y1)
y2
Slope=f(x1,y1)
hf(x1,y1)
Slope=f(x0,y0)
y1
hf(x0,y0)
y0
x0
Fin500J Topic 7
y1=y0+hf(x0,y0)
h
x1
h
Fall 2010 Olin Business School
x2
x
9
High Order Taylor Series methods
dy( x)
Given
 f ( y, x), y ( x0 )  y0
dx
th
n order T aylorSeries method
dy h 2 d 2 y
hn d n y
n 1
yi 1  yi  h 
 ....
 O(h )
2
n
dx 2! dx
n! dx
2
3
n
d y d y
d y
, 3 ,....., n need to be derivedanalatically.
2
dx dx
dx
Fin500J Topic 7
Fall 2010 Olin Business School
10
Runge-Kutta Methods (Motivation)


We seek accurate methods to solve ODE that does
not require calculating high order derivatives.
The approach is to a formula involving unknown
coefficients then determine these coefficients to
match as many terms of the Taylor series expansion
Fin500J Topic 7
Fall 2010 Olin Business School
11
Runge-Kutta Method
Second Order Runge Kutta
K1  h f (t , x)
K 2  h f (t   h, x   K1 )
x(t  h)  x(t )  w1 K1  w2 K 2
P roblem:
Find  ,  , w1 , w2
such thatx(t  h) is as accurat eas possible.
Fin500J Topic 7
Fall 2010 Olin Business School
12
Taylor Series in One Variable
T heT aylorSeries expansionof f(x)
n 1
i
h (i )
f ( x  h)   f ( x )
i  0 i!
Approximation

n
h (n)
f (x)
n!
Error
where x is between x and x  h
Fin500J Topic 7
Fall 2010 Olin Business School
13
Taylor Series in One Variable
another look
Define
i
 d 
i d f ( x)
(i )
i
 f ( x) h
 h  f ( x)  h
i
dx
 dx 
T heT aylorSeries expansionof f(x)
i
n 1
i
n
1 d 
1 d 
f ( x  h)    h  f ( x )   h  f ( x )
n!  dx 
i  0 i!  dx 
x is between x and x  h
Fin500J Topic 7
Fall 2010 Olin Business School
14
Definitions
Define
i

f
  
i
h
f
(
x
,
y
)

h


x i
 x 
i
0
 
 
 h
 f ( x, y )  f ( x, y )
k
y 
 x
1
 
 
 f ( x, y )
 f ( x, y )
 h
 f ( x, y )  h
k
k
y 
x
y
 x
2
2
2
 
 
f ( x, y )
 2 f ( x, y )
f ( x, y )
2 
2 
 h

k
f ( x, y )  h
 2h k
k
2
2


x

y

x

x

y

y


Fin500J Topic 7
Fall 2010 Olin Business School
15
Taylor Series Expansion
f(x, y)  ( x  1)(x  y  2)
2
Parialderivatives evaluatedat (0,0)
0
 
 
 h  k  f ( x, y )
4
y 
 x
( 0.0 )
1
 
 
 h  k  f ( x, y )
 h fx  k f y
y 
 x
( 0.0 )
Fin500J Topic 7
Fall 2010 Olin Business School
( 0,0)
16
Taylor Series in Two Variables
T heT aylorSeries expansionof f(x, y)
i
n
1 
 
1 
 
f ( x  h, y  k )    h  k  f ( x, y )   h  k  f ( x , y )
y 
n!  x
y 
i  0 i!  x
n 1
approxim ation
error
( x , y ) is on theline joiningbetween ( x, y ) and ( x  h, y  k )
y+k
y
x
Fin500J Topic 7
x+h
Fall 2010 Olin Business School
17
Runge-Kutta Method
Second Order Runge Kutta
K1  h f (t , x)
K 2  h f (t   h, x   K1 )
x(t  h)  x(t )  w1 K1  w2 K 2
P roblem:
Find  ,  , w1 , w2
such thatx(t  h) is as accurat eas possible.
Fin500J Topic 7
Fall 2010 Olin Business School
18
Runge-Kutta Method
P roblem: Find  ,  , w1 , w2
to matchas many termsas possible.
h2
h3
x(t  h)  x(t )  hx' (t )  x' ' (t )  x' ' ' (t )  ...
2
6
x(t  h)  x(t )  w1h f (t , x)  w2 h f (t   h, x   h f (t , x))
1

 
f (t   h, x   h f )  f   h f t   h f f x    h   hf  f (t , x )
2
t
x 
x(t  h)  x(t )  w1  w2 h f (t , x)   w2 h 2 f t   w2 h 2 f f x  O(h 3 )
2
Fin500J Topic 7
Fall 2010 Olin Business School
19
Runge-Kutta Method
x(t  h)  x(t ) 
hx' (t ) 
h2
x' ' (t )
2
h3

x' ' ' (t )  ...
6
x(t  h)  x(t )  w1  w2 h f (t , x)   w2 h 2 f t   w2 h 2 f f x  O(h3 )
 w1  w2  1,  w2  0.5,
One possible solution
w1  0.5, w2  0.5,   1,
Fin500J Topic 7
Fall 2010 Olin Business School
 w2  0.5
 1
20
Runge-Kutta Method
Second Order Runge Kutta
Alt ernat ive Form
K1  h f (t , x)
F1  f (t , x)
K 2  h f (t  h, x  K1 )
F2  f (t  h, x  hF1 )
1
x(t  h)  x(t )  K1  K 2 
2
h
x(t  h)  x(t )  F1  F2 
2
Fin500J Topic 7
Fall 2010 Olin Business School
21
Runge-Kutta Method
Alternative Formulas
 w1  w2  1,  w2  0.5,
 w2  0.5
another solution
Pick  any non - zero number
  ,
1
w1  1 
,
2
1
w2 
2
Second Order Runge Kutta Formulas(select  0)
K1  h f (t , x)
K 2  h f (t   h, x   K1 )
1 
1

x(t  h)  x(t )  1 
F2
 F1 
2
 2 
Fin500J Topic 7
Fall 2010 Olin Business School
22
Runge-Kutta Methods
T hirdOrder Runge Kutta(RK3)
K1  f ( xi , yi )
1
1
K 2  f ( xi  h, yi  K1h)
2
2
1
K 3  f ( xi  h, yi  K1h  2 K 2 h)
2
1
y ( x  h )  y ( x )   K1  4 K 2  K 3 
6
Fin500J Topic 7
Fall 2010 Olin Business School
23
Runge-Kutta Methods
Fourt hOrder Runge Kut t a (RK4)
K1  f ( xi , yi )
1
1
K 2  f ( xi  h, yi  K1h)
2
2
1
1
K 3  f ( xi  h, yi  K 2 h)
2
2
K 4  f ( xi  h, yi  K 3 h)
h
yi 1  yi  K1  2 K 2  2 K 3  K 4 
6
Fin500J Topic 7
Fall 2010 Olin Business School
24
Runge-Kutta Methods
Higher order Runge-Kutta methods are available
Higher order methods are more accurate but
require more calculations.
Fourth order is a good choice. It offers good
accuracy with reasonable calculation effort
Fin500J Topic 7
Fall 2010 Olin Business School
25
Example 1
Second Order Runge-Kutta Method
P roblem:
dy
 1  y 2  x3 ,
y (1)  4
dx
Use RK 2 to find y (1.01), y (1.02)
h  0.01
f ( x, y )  1  y 2  x 3
x0  1, y0  4
Step 1 :
K1  f ( x0 , y0 )  (1  y0  x0 )  18.0
2
3
K 2  f ( x0  h, y0  K1h)  (1  ( y0  0.18) 2  ( x0  .01)3 )  16.92
h
0.01
y1  y0  K1  K 2   4 
(18  16.92)  3.8254
2
2
Fin500J Topic 7
Fall 2010 Olin Business School
26
Example 1
Second Order Runge-Kutta Method
P roblem:
dy
 1  y 2  x3 ,
y (1)  4
dx
Use RK 2 to find y (1.01), y (1.02)
h  0.01
f ( x, y)  1  y 2  x 3
x1  1.01, y1  3.8254
Step 2 :
K1  f ( x1 , y1 )  (1  y1  x1 )  16.66
2
3
K 2  f ( x1  h, y1  K1h)  (1  ( y1  0.1666) 2  ( x1  .01)3 )  15.45
y2  y1 
Fin500J Topic 7
h
K1  K 2   3.8254 0.01 (16.66  15.45)  3.6648
2
2
Fall 2010 Olin Business School
27
Example 1
Summary of the solution
P roblem:
dy
 1  y 2  x3 ,
y (1)  4
dx
Use RK 2 to find y (1.01), y (1.02)
Summary of the solution
i
xi
yi
0
1
2
Fin500J Topic 7
1.00
1.01
1.02
 4.0000
 3.8254
 3.6648
Fall 2010 Olin Business School
28
Solution after 100 steps
Fin500J Topic 7
Fall 2010 Olin Business School
29
Numerically Solving ODE in Matlab
 Matlab has a few different ODE solvers, Matlab
recommends ode45 is used as a first solver for a
problem
 ode45 uses simultaneously fourth and fifth order
Runge-Kutta formula (Dormand–Prince)
Fin500J Topic 7
Fall 2010 Olin Business School
30
Numerically Solving ODE in Matlab (Example 1)
P roblem:
dy
 1  y 2  x3 ,
dx
y (1)  4, over t heinterval[1,2].
 Step 1: Create a M-file for dy/dx as firstode.m
function yprime=firstode(x,y);
yprime=1+y^2+x^3;
 Step 2: At a Matlab command window
>>[x,y]=ode45(@firstode,[1,2],-4);
>> [x,y]
Matlab returns two column vectors, the first with values of x and
the second with value of y.
Fin500J Topic 7
Fall 2010 Olin Business School
31
Numerically Solving ODE in Matlab (Example 1)
>> plot(x,y,'+')
5
4
3
2
1
0
-1
-2
-3
-4
Fin500J Topic 7
1
1.1
1.2
1.3
1.4
1.5
Fall 2010 Olin Business School
1.6
1.7
1.8
1.9
2
32
Solving a system of first order ODEs

Methods discussed earlier such as Euler, RungeKutta,…are used to solve first order ordinary
differential equations

The same formulas will be used to solve a system of
first order ODEs. In this case, the differential
equation is a vector equation and the dependent
variable is a vector variable.
Fin500J Topic 7
Fall 2010 Olin Business School
33
Euler method for solving a system of first
order ODEs
Recall Euler method for solving first order ODE.
dy( x)
Given
 f ( y, x), y (a)  ya
dx
Euler Method :
y (a  h)  y (a)  h f ( y (a), a)
y (a  2h)  y (a  h)  h f ( y (a  h), a  h)
y (a  3h)  y (a  2h)  h f ( y (a  2h), a  2h)
Fin500J Topic 7
Fall 2010 Olin Business School
34
Solving a system of n first order ODEs
using Euler method



Exactly the same formula
is used but the scalar
variables and functions
are replaced by vector
variables and vector
values functions.
Y is a vector of length n
F(Y,x) is vector valued
function
 y1 ( x) 
 y ( x)
Y ( x)   2 
Y is n 1 vector
 ... 


y
(
x
)
 n 
 d y1 
 dx   f1 (Y , x) 
d y  

d Y ( x)  2   f 2 (Y , x)
 dx 
 F (Y , x)


 ... 
dx
 ...  

 d yn   f n (Y , x)
 dx 
Y (a  h)  Y (a )  h F (Y (a ), a )
Y (a  2h)  Y (a  h)  h F (Y (a  h), a  h)
Y (a  3h)  Y (a  2h)  h F (Y (a  2h), a  2h)
Fin500J Topic 7
Fall 2010 Olin Business School
35
Example :
Euler method for solving a system of first order ODEs
 y '1 ( x)   y2 
 y1 (0)   1
 y ' ( x)  1  y   F (Y , x), Y (0)   y (0)   1 
1
 2  
 2   
Two steps of Euler Method with h  0.1
STEP1 :
Y (0  h)  Y (0)  h F (Y (0),0)
 y1 (0.1)   y1 (0) 
 y2 (0)    1  0.1   0.9
 y (0.1)   y (0)  0.11  y (0)  1  0.1(1  1)   1.2 
 

1
 2
  2 

 
STEP 2 :
Y (0  2h)  Y (h)  h F (Y (h), h)
 y1 (0.2)   y1 (0.1) 
 y2 (0.1)    0.9  0.12   0.78
 y (0.2)   y (0.1)  0.11  y (0.1)  1.2  .1(1  0.9)   1.39 
 

1
 2
  2


 
Fin500J Topic 7
Fall 2010 Olin Business School
36
Example :
RK2 method for solving a system of first order ODEs
 y '1 ( x)   y2 
 y1 (0)   1


F
(
Y
,
x
),
Y
(
0
)

 y ' ( x) 1  y 
 y (0)   1 
1
 2  
 2   
Two steps of secondorder Runge  Kut ta Met hod with h  0.1
STEP1 :
 y2 (0)   0.1
K1  h F (Y (0),0)  0.1
 

1  y1 (0) 0.2
 y (0)  0.2  0.12
K 2  h F (Y (0)  K1,0  h)  0.1 2



1

(
y
(
0
)

0
.
1
)
0
.
19

1

 
Y (0  h)  Y (0)  0.5( K1  K 2)
 y1 (0.1)   1 1   0.1 0.12   0.89
 y (0.1)   1   2  0.2  0.19    1.195 
 

 2
  
  
Fin500J Topic 7
Fall 2010 Olin Business School
37
Example :
RK2 method for solving a system of first order ODEs
 y '1 ( x)   y2 
 y1 (0)   1
 y ' ( x)  1  y   F (Y , x), Y (0)   y (0)   1 
1
 2  
 2   
STEP 2 :
 y2 (0.1)  0.1195
K1  h F (Y (0.1),0.1)  0.1



1

y
(
0
.
1
)
0
.
1890

1

 
 y2 (0.1)  0.189  0.1384
K 2  h F (Y (0.1)  K1,0.1  h)  0.1



1

(
y
(
0
.
1
)

0
.
1195
)
0
.
1771

1

 
Y (0.1  h)  Y (0.1)  0.5( K1  K 2)
 y1 (0.2)   0.89 1  0.1195 0.1384   0.7611
 y (0.2)   1.195   2  0.1890   0.1771    1.3780 

 
 

 2
 

Fin500J Topic 7
Fall 2010 Olin Business School
38
The general approach to solve high
order ODE
high order ODE
y ' '3 y '6 y  1
y ' (0)  1; y (0)  4
Second order ODE
Fin500J Topic 7
convert
System of first order ODE
solve
convert
z2
 z1 '  

 z '  1  3z  6 z ,
2
1
 2 
 4
Z (0)   
1 
solve
Two first order ODEs
Fall 2010 Olin Business School
39
Conversion Procedure
high order ODE
1.
convert
System of first order ODE
solve
Select of dependent variables
One way is to take the original dependent variable and its
derivatives up to one degree less than the highest order
derivative.
2. Write the Differential Equations in terms of the new
variables. The equations comes from the way the new
variables are defined or from the original equation.
3. Express the equations in matrix form
Fin500J Topic 7
Fall 2010 Olin Business School
40
Example of converting High order ODE to
first order ODEs
Convert y ' '3 y '6 y  1,
y ' (0)  1; y (0)  4
t o a systemof first order ODE
1.Select a new set of variable
(secondorder ODE  We need two variables)
z1  y
z2  y'
Fin500J Topic 7
One degree less than the
highest order derivative
Fall 2010 Olin Business School
41
Example of converting High order ODE to
first order ODEs
old
new Initial
name name cond.
Equation
y
z1
4
z '1  z 2
y'
z2
1
z '2  1  3 z 2  6 z1
z2
 z '1  

 4
 z '   1  3 z  6 z , Z (0)  1 
 
2
1
 2 
Fin500J Topic 7
Fall 2010 Olin Business School
42
Numerically Solving high order ODE in Matlab
Example:
y' '   xy  e x y'3 sin 2 x,
y(0)  2, y' (0)  8, over theinterval[0,4].
 Step 1: First convert the second order equation to an
equivalent system of first order ODEs, let z1=y, z2=y’:
z2
 z1 '  

 
,

x
 z 2 '  xz1  e z 2  3 sin 2 x 
 2
Z (0)   
8 
Fin500J Topic 7
Fall 2010 Olin Business School
43
Numerically Solving high order ODE in Matlab
Example:
y' '   xy  e x y'3 sin 2 x,
y(0)  2, y' (0)  8, over theinterval[0,1].
 Step 2: Create the following M-file and save it as F.m
function zprime=F(x,z)
zprime=zeros(2,1); %since output must be a column vector
zprime(1)=z(2);
zprime(2)=-x*z(1)+exp(x)*z(2)+3*sin(2*x);
 Step 3: At Matlab prompt
>> [x,z]=ode45(@F,[0,1],[2,8])
Since z1(x)=y, to print out the solution y
>> [x,z(:,1)]
Fin500J Topic 7
Fall 2010 Olin Business School
44
Numerically Solving ODE in Matlab (Example 1)
20
To plot y against x
>> plot(x, z(:,1))
18
16
14
Because the
vector z has
first
component
z1=y
12
10
8
6
4
2
Fin500J Topic 7
0
0.1
0.2
0.3
0.4
0.5
Fall 2010 Olin Business School
0.6
0.7
0.8
0.9
1
45