CVEN 302, Chapter 12 - Texas A&M University
Download
Report
Transcript CVEN 302, Chapter 12 - Texas A&M University
Chapter 20
ODEs:
Initial-Value Problems
Differential Equations
There are ordinary differential equations functions of one variable
cd 2
dv
g v
dt
m
And there are partial differential equations functions of multiple variables
u
u
2u
u
2
t
x
x
2 2
0
2
y
x
Order of differential equations
1st order (falling parachutist)
cd 2
dv
g v
dt
m
2nd order (mass-spring system with damping)
etc.
d2x
dx
m 2 c
kx 0
dt
dt
Higher-Order ODE’s
Can always turn a higher order ODE into a set of 1st
order ODEs
Example:
Let
d3x
d2x
dx
a 3 b 2 c
d t
dt
dt
dt
dx
d2x
y
, z 2
dt
dt
then
dz 1
dt a d ( t ) bz cy
dy
z
dt
dx
dt y
So solutions to first-order ODEs are important
Linear and Nonlinear ODEs
Linear: No multiplicative mixing of variables,
no nonlinear functions
Nonlinear: anything else
d 2 g
sin 0
2
l
dt
u
u
2u
u
2
t
x
x
Ordinary Differential Equations
ODEs show up everywhere in engineering
Dynamics (Newton’s 2nd law)
dv F
dt m
dT
Heat conduction (Fourier’s law) q k
dx
Diffusion (Fick’s law)
dc
J D
dx
•Example of an
ODE
m
d 2x
dt 2
c
dx
kx f (t )
dt
x
k
What order is this ODE?
m
If f(t) = 0, ODE is homogenous.
If f(t) is not equal to 0,
ODE is non-homogenous.
f(t)
c
Free-body
diagram
kx
m
f(t)
dx
c
dt
d 2x
dx
m
c
kx f (t )
2
dt
dt
Solutions for ODEs
The solution for the homogenous ODE.
x(t ) C1S1 (t ) C2 S 2 (t )
C1 & C2 are two arbitrary constants and
S1 (t ) & S 2 (t ) are two general solutions.
The solution for the non-homogenous ODE
x(t ) C1S1 (t ) C2 S 2 (t ) P (t )
P(t ) is the particular solutions.
The arbitrary constants C1 & C2 are determined by
the Initial-value or Boundary-value conditions.
Initial-Value & Boundary-Value Conditions
The I-V Conditions
All conditions are given
at the same value of the independent variable.
dx
For example at t 0, x 0 &
0.
dt
The two conditions are all given at t 0.
The B-V Conditions
Conditions are given
at the different values of the independent variable.
For example at t 0, x 1 & at t 1, x 0.5.
The numerical schemes for solving Initial-value and
boundary-value are different.
Ordinary Differential Equations
I-V Problems
Euler’s and Heun's methods
Runge-Kutta methods
Adaptive Runge-Kutta
Multistep methods*
Adams-Bashforth-Moulton methods*
Ordinary Differential Equations
1st order Ordinary differential equations (ODEs)
Initial value problems
dy
f ( t , y) ; y( t0 ) y0
dt
Numerical approximations
New value = old value + slope × step size
yi 1 yi h
Runge-Kutta Methods
Runge Kutta methods (One Step Methods)
Idea is that
New value = old value + slope*step size
yi 1 yi h
Slope is generally a function of t, hence y(t)
Different methods differ in how to estimate
One-Step Method
yi 1 yi h
All one-step methods
can be expressed in
this general form, the
only difference being
the manner in which
the slope is estimated
Initial Conditions
The solution of
ODE depends on
the initial condition
t
Same ODE, but
with different
initial conditions
Euler’s method
(First-order Taylor Series Method)
Approximate the derivative by finite difference
dy yi 1 yi
f ( t i , yi ) yi 1 yi hf ( t i , yi )
dt
ti 1 ti
Local truncation error
yi 1
or
yi 1
yi 2
yi( n ) n
yi yih h
h Rn ;
2
n!
Rn
y ( n 1) ( )
( n 1)!
h n 1
f ( t i , yi ) i 2
f ( n1) ( t i , yi ) n
yi f ( t i , yi )h
h
h O ( h n 1 )
2
n!
Euler’s Method
Euler’s (Euler-Cauchy or point-slope) method
Use the slope at ti to predict yi+1
Euler’s method
dy
y f ( t , y );
dt
y(t 0 ) y0
Straight line approximation
y0
t0
h
t1
h
t2
h
t3
Example: Euler’s Method
dy
t y , y(0 ) 1,
dt
0t1
Analytic solution y ( 1 t 2 / 4 ) 2
Euler method
yi 1 yi hfi , f t y
h = 0.50
y(0.5 ) y(0 ) hf (0,1.0 ) 1 (0.5)(0 1 ) 1.0
y( 1.0 ) y(0.5 ) hf (0.5 ,1.0 ) 1.0 (0.5 )( 0.5 1.0 ) 1.25
Example: Euler’s Method
h = 0.25
yi 1 yi hfi , f t y
y(0.25 ) y(0 ) hf (0, 1.0 )
1 (0.25 )( 0 1 ) 1.0
y(0.5 ) y(0.25 ) hf (0.25 , 1.0 )
1.0 (0.25 )( 0.25 1.0 ) 1.0625
y(0.75 ) y(0.5 ) hf (0.5 , 1.0625 )
1.0625 (0.25 )( 0.5 1.0625 ) 1.191347
y( 1.0 ) y(0.75 ) hf (0.75 , 1.191347 )
1.191347 (0.25 )( 0.75 1.191347 ) 1.39600
Euler’s method:
y
dy
y t y ;
dt
y( 0 ) 1
h = 0.25
1
h = 0.5
0
0.25
0.5
0.75
1.0
t
Euler’s Method (modified M-file)
Euler’s Method
>> tt=0:0.01*pi:pi;
>> ye=1.5*exp(-tt)+0.5*sin(tt)-0.5*cos(tt);
>> [t,y]=Eulode('example2_f',[0 pi],1,0.1*pi);
step
t
y
1
0.0000000000
1.0000000000
2
0.3141592654
0.6858407346
3
0.6283185307
0.5674580652
dy
4
0.9424777961
0.5738440394
dt
5
1.2566370614
0.6477258022
6
1.5707963268
0.7430199565
7
1.8849555922
0.8237526182
8
2.1991148575
0.8637463173
9
2.5132741229
0.8465525934
10
2.8274333882
0.7652584356
11
3.1415926536
0.6219259596
>> H=plot(t,y,'r-o',tt,ye);
>> set(H,'LineWidth',3,'MarkerSize',12)
>> print -djpeg ode01.jpg
y sin t ; y(0 ) 1
Euler’s Method (h = 0.1)
dy
y sin t ; y(0 ) 1
dt
Euler’s Method (h = 0.05)
dy
y sin t ; y(0 ) 1
dt
Truncation Errors
There are
Local truncation errors - error from
application at a single step
Propagated truncation errors - previous
errors carried forward
The sum is “global truncation error”
Global & Local Errors
y
y
yi+1
yi
o
Global error
Local error
yi+1
yi
xi
xi+1
x
o
xi
xi+1
xi+2
x
Euler’s Method
Euler’s method uses Taylor series with only
first order terms
The true local truncation error is
f ( t i , yi ) 2
Et
h ... O( hn1 )
2!
Approximate local truncation error - neglect
higher order terms (for sufficiently small h)
f ( t i , yi ) 2
Ea
h O( h2 )
2!
Runge-Kutta Methods
Higher-order Taylor series methods (see
Chapra and Canale, 2002) -- need to compute
the derivatives of f(t,y)
Runge-Kutta Methods -- estimate the slope
without evaluating the exact derivatives
Heun’s method
Midpoint (or improved polygon) method
Third-order Runge-Kutta methods
Fourth-order Runge-Kutta methods
Heun’s Method
Improvements of Euler’s method - Heun’s
method
Euler’s method - derivative at the beginning
of interval is applied to the entire interval
Heun’s method uses average derivative for
the entire interval
A second-order Runge-Kutta Method
Heun’s Method
Heun’s method is a predictor-corrector method
Predictor
y
0
i1
yi f ( t i , yi ) h
Corrector (may be applied iteratively)
yi 1
f ( t i , yi ) f ( t i 1 , y )
yi
h
2
0
i1
Heun’s Method
Iterate the corrector of Heun’s method to obtain
an improved estimate
Predictor
Corrector
Heun’s Method with Iterative Correctors
Heun’s Method with Iterative Correctors
»
»
»
»
»
»
»
te=0:0.02*pi:pi; ye=example2_e(xe);
[t1,y1]=Euler('example2_f',[0 pi],1,0.1*pi);
[t2,y2]=Heun_iter('example2_f',[0 pi],1,0.1*pi,0);
[t3,y3]=Heun_iter('example2_f',[0 pi],1,0.1*pi,5);
H=plot(te,ye,t1,y1,'r-d',t2,y2,'g-s',t3,y3,'m-o');
set(H,'LineWidth',3,'MarkerSize',12)
[te' ye' y1',y2',y3']
t
yEuler
yHeun
yHeun_iter
ytrue
0
1.0000
1.0000
1.0000
1.0000
0.3142
0.6858
0.7837
0.7704
0.7746
0.6283
0.5675
0.7018
0.6830
0.6896
0.9425
0.5738
0.7064
0.6872
0.6951
1.2566
0.6477
0.7559
0.7395
0.7479
1.5708
0.7430
0.8152
0.8036
0.8118
1.8850
0.8238
0.8565
0.8503
0.8578
2.1991
0.8637
0.8592
0.8584
0.8648
2.5133
0.8466
0.8112
0.8149
0.8199
2.8274
0.7653
0.7082
0.7154
0.7188
3.1416
0.6219
0.5540
0.5631
0.5648
Heun’s Method with Iterative Correctors
dy
y sin( t ) ; y(0 ) 1
dt
Heun’s with
5 iterations
Heun’s
method
Exact
Euler’s
t
Example: Heun’s Method
h = 0.5
First Step
Preditor :
dy
t y , y(0 ) 1,
dt
0t1
y10 y0 f ( t0 , y0 ) h 1 (0 1 )( 0.5 ) 1.0
f ( t0 , y0 ) f ( t 1 , y10 )
1st Corrector : y y0
h 1 0.5 (0
2
f ( t0 , y0 ) f ( t 1 , y11 )
2
2nd Corrector : y1 y0
h 1 0.5 (0
2
f ( t0 , y0 ) f ( t 1 , y12 )
3
3rd Corrector : y1 y0
h 1 0.5 (0
2
f ( t0 , y0 ) f ( t 1 , y13 )
4
4th Corrector : y1 y0
h 1 0.5 (0
2
f ( t0 , y0 ) f ( t 1 , y14 )
5
5 th Corrector : y1 y0
h 1 0.5 (0
2
1
1
1 0.5 1.0 )( 0.5 ) 1.1250
1 0.5 1.1250 )( 0.5 ) 1.1326
1 0.5 1.1326 )( 0.5 ) 1.1330
1 0.5 1.1330 )( 0.5 ) 1.1331
1 0.5 1.1331 )( 0.5 ) 1.1331
Example: Heun’s Method
h = 0.5
dy
t y , y(0 ) 1,
dt
0x1
Second Step
Preditor :
y20 y1 f ( t 1 , y1 ) h 1.1331 (0.5 1.1331 )( 0.5 ) 1.3992
f ( t 1 , y1 ) f ( t 2 , y20 )
1st Corrector : y y1
h 1 0.5 (0.5
2
f ( t 1 , y1 ) f ( t 2 , y21 )
2
2nd Corrector : y2 y1
h 1 0.5 (0.5
2
f ( t 1 , y1 ) f ( t 2 , y22 )
3
3rd Corrector : y2 y1
h 1 0.5 (0.5
2
f ( t 1 , y1 ) f ( t 2 , y23 )
4
4th Corrector : y2 y1
h 1 0.5 (0.5
2
f ( t 1 , y1 ) f ( t 2 , y24 )
5
5 th Corrector : y2 y1
h 1 0.5 (0.5
2
1
2
1.1331 1.0 1.3992 )( 0.5 ) 1.5619
1.1331 1.0 1.5619 )( 0.5 ) 1.5786
1.1331 1.0 1.5786 )( 0.5 ) 1.5803
1.1331 1.0 1.5803 )( 0.5 ) 1.5804
1.1331 1.0 1.5804 )( 0.5 ) 1.5804
Example:
dy
t y , y(0 ) 1,
dt
0t5
» [t1,y1]=Eulode('example3',[0 5],1,0.5);
» [t2,y2]=Heun_iter('example3',[0,5],1,0.5,0);
» [t3,y3]=Heun_iter('example3',[0 5],1,0.5,5);
» [t1' y1' y2' y3' ye']
t
0
0.5000
1.0000
1.5000
2.0000
2.5000
3.0000
3.5000
4.0000
4.5000
5.0000
yEuler
1.0000
1.0000
1.2500
1.8090
2.8178
4.4964
7.1470
11.1570
17.0024
25.2492
36.5552
yHeun
1.0000
1.1250
1.5523
2.4169
3.9463
6.4619
10.3793
16.2082
24.5532
36.1126
51.6796
yHeun_iter
1.0000
1.1331
1.5804
2.4859
4.0882
6.7192
10.8046
16.8630
25.5065
37.4407
53.4643
ytrue
1.0000
1.1289
1.5625
2.4414
4.0000
6.5664
10.5625
16.5039
25.0000
36.7539
52.5625
Midpoint Method
Improved Polygon or Modified Euler Method
Use the slope at midpoint to represent the average slope
h
yi 1 / 2 yi f ( t i , yi ) ; yi 1 / 2 f ( t i 1 / 2 , yi 1 / 2 )
2
yi 1 yi f ( t i 1 / 2 , yi 1 / 2 )h
Runge-Kutta (RK) Methods
One-Step Method with general form
yi 1 yi ( t i , yi , h) h
Where is an increment function which represents
the weighted-average slope over the interval
a1 k1 a2 k2 an kn
Where a’s are constants and k’s are slopes evaluated
at selected x locations
Runge-Kutta (RK) Methods
One-Step Method
General Form of nth-order Runge-Kutta Method
yi 1 yi h
a1 k 1 a 2 k 2 a n k n
k1
k2
k3
kn
f ( t i , yi )
f ( t i p1 h, yi q11k1 h)
f ( t i p2 h, yi q21k1 h q22 k 2 h)
f ( t i pn 1 h, yi qn 1, 1 k1 h qn 1, 2 k 2 h qn 1, n 1 kn 1 h)
Where p’s and q’s are constants
k’s are recurrence relationships
Second-Order RK Methods
Taylor series expansion
k1 f ( t , y )
k 2 f ( t p1 h, y q11k1 h) f ( t , y ) ( p1 h) f t ( t , y ) ( q11k1 h) f y ( t , y )
y( t h) y( t ) ha1 f ( t , y ) a2 f ( t p1 h, y q11k1 h)
y( t ) ha1 f a2 ( f p1 hft q11k1 hf y )
Compare to the second-order Taylor formula
y( t h) y( t ) hf ( h2 /2 )( ft ff y )
Three equations for four unknowns (a1, a2, k1, q11)
a1 a2 1
a1 a2 1
2
2
a
p
h
f
(h
/2)f
2 1
a2 p1 1/2
t
t
a q 1/2
2
2
2
a
q
k
h
f
a
q
h
ff
(h
/2)f
f
2 11
y
2 11
y
y
2 11 1
Second-Order RK Methods
Second-order version of Runge-Kutta Methods
yi 1 yi ( a1 k1 a2 k2 )h
k1 f ( t i , yi )
k2 f ( t i p1 h, yi q11k1 h)
3 equations for 4 unknowns
a1 a 2 1
a 2 p1 1 / 2
a q 1 / 2
2 11
Second-Order RK Methods
General Second-order Runge-Kutta methods
k2
k1
a1 k1 + a2 k2
k1 f ( t i , yi ) f1
k2 f ( t i p1 h, yi q11k1 h)
y y ( a k a k )h
i
1 1
2 2
i1
Weighted-average slope
xi
xi+p1h
xi+1
= xi+h
Heun’s Method
Heun’s method with a single corrector (a2 = 1/2):
Choose a2 = 1/2 and solve for the other 3 constants
a1 = 1/2, p1 = 1, q11 = 1
k 1 f ( t i , yi )
k 2 f ( t i h, yi k1 h)
yi 1
1
1
yi k 1 k 2 h
2
2
Use average slope over the interval
Midpoint Method
Another Second-order Runge-Kutta method
Choose a2 = 1 a1 = 0, p1 = 1/2, q11 = 1/2
k2
k 1 f ( t i , yi )
h
1
k2 f ( t i , yi k1 h)
2
2
yi 1 yi k2 h
k1
xi
xi+1/2
xi+1
Midpoint (2nd-order RK) Method
Midpoint (2nd-order RK) Method
>> [t,y] = midpoint('example2_f',[0 pi],1,0.05*pi);
step
t
y
1
0.0000000000
1.0000000000
2
0.1570796327
0.8675816988
3
0.3141592654
0.7767452235
4
0.4712388980
0.7206165079
5
0.6283185307
0.6927855807
6
0.7853981634
0.6872735266
7
0.9424777961
0.6985164603
8
1.0995574288
0.7213629094
9
1.2566370614
0.7510812520
dy
10
1.4137166941
0.7833740988
y sin( t )
11
1.5707963268
0.8143967658
dt
12
1.7278759595
0.8407772414
13
1.8849555922
0.8596353281
y( 0 ) 1
14
2.0420352248
0.8685989204
15
2.1991148575
0.8658156821
16
2.3561944902
0.8499586883
17
2.5132741229
0.8202249154
18
2.6703537556
0.7763257790
19
2.8274333882
0.7184692359
20
2.9845130209
0.6473332788
21
3.1415926536
0.5640309524
>> tt = 0:0.01*pi:pi; yy = example2_e(tt);
>> H = plot(t,y,'r-o',tt,yy);
>> set(H,'LineWidth',3,'MarkerSize',12);
Midpoint (RK2) Method
dy
y sin( t ) ; y(0 ) 1
dt
Ralston’s Method
Second-order Runge-Kutta method
Choose a2 = 2/3 a1 = 1/3, p1 = q11 = 3/4
k2
k 1 f ( t i , yi )
k1
k1/3 + 2k2/3
ti
ti+3h/4
3
3
k 2 f ( t i h, yi k1 h)
4
4
1
2
yi 1 yi ( k 1 k 2 ) h
3
3
xi+1 = ti+h
Third-Order Runge-Kutta Method
General form
k 1 f ( t i , yi )
k 2 f ( t i p1 h, yi q11k1 h)
k 3 f ( t i p2 h, yi q21k1 h q22 k 2 h)
yi 1 yi ( a1 k1 a2 k 2 a3 k 3 ) h
Weighted slope
a1 k1 a2 k2 a3 k3
Third-Order Runge-Kutta Methods
General Third-order Runge-Kutta methods
k3
k2
Weightedaverage
value of
three slopes
k1 , k 2 , k 3
k1
ti
ti+p1h
ti+p2h
ti+1 = ti+h
Third-order Runge-Kutta Methods
k 1 f ( t i , yi )
2
2
k 2 f ( t i h, yi k1 h)
3
3
2
2
k 3 f ( t i 3 h, yi 3 k2 h)
1
yi 1 yi ( 2 k 1 3 k 2 3 k 3 ) h
8
Nystrom Method
k 1 f ( t i , yi )
1
1
k 2 f ( t i h, yi k1 h)
2
2
3
3
k 3 f ( t i 4 h, yi 4 k 2 h)
1
yi 1 yi ( 2 k 1 3 k 2 4 k 3 ) h
9
Nearly Optimum Method
3rd-order
Runge-Kutta Method
Reduce to Simpson’s 1/3 rule for f = f(t)
k 1 f ( t i , yi )
1
1
k 2 f ( t i h, yi k1 h)
2
2
k 3 f ( t i h, yi k1 h 2 k 2 h)
yi 1
1
yi ( k 1 4 k 2 k 3 ) h
6
rd
3 -Order
Heun Method
k 1 f ( t i , yi )
1
1
k 2 f ( t i h, yi k1 h)
3
3
2
2
k 3 f ( t i 3 h, yi 3 k 2 h)
1
yi 1 yi ( k 1 3 k 3 ) h
4
Classical 4th-order
Runge-Kutta Method
One-step method
Reduce to Simpson’s 1/3 rule for f = f(t)
k1
k2
k
3
k
4
yi 1
f ( t i , yi )
1
1
f ( t i h, yi k1 h)
2
2
1
1
f ( t i h, yi k 2 h)
2
2
f ( t i h, yi k 3 h)
1
yi ( k 1 2 k 2 2 k 3 k 4 ) h
6
Classical 4th-order
Runge-Kutta Method
k2
k4
k3
k1
ti
ti + h/2
1
( k1 2 k 2 2 k 3 k 4 )
6
ti + h
Example: Classical
4th-order RK Method
h = 0.5
dy
t y , y(0) 1,
dt
0t1
First step, t = 0.5
k1
k2
k3
k4
f ( t0 , y0 ) 0 1.0 0
f ( t0 0.5 h, y0 0.5 k1 h) f (0.25 , 1.0 ) (0.25 1 ) 0.25
f ( t0 0.5 h, y0 0.5 k2 h) f (0.25 , 1.0625 ) 0.25 1.0625 0.257694
f ( t0 h, y0 k3 h) f (0.5 , 1.128847 ) 0.5 1.128847 0.531236
1
1
( k1 2 k2 2 k3 k4 )h 1 0 2(0.25 ) 2(0.257694 ) 0.531236 (0.5 )
6
6
1.128885 ( t 0.0018%)
y1 y0
Example: Classical
4th-order RK Method
Second step, t = 1.0
k1
k2
k3
k4
f ( t 1 , y1 ) 0.5 1.128885 0.531236
f ( t 1 0.5 h, y1 0.5 k1 h) f (0.75 , 1.261697 ) 0.75 1.261697 0.8424396
f ( t 1 0.5 h, y1 0.5 k2 h) f (0.75 , 1.339495 ) 0.75 1.339495 0.8680242
f ( t 1 h, y1 k3 h) f ( 1.0 , 1.562897 ) 1.0 1.562897 1.2501588
1
( k1 2 k 2 2 k 3 k 4 ) h
6
1
1.128885 0.531236 2(0.8424396 ) 2(0.8680242 ) 1.2501588 (0.5 )
6
1.5624127 ( t 0.0056%)
y1 y0
Fourth-Order Runge-Kutta Method
Fourth-Order Runge-Kutta Method
>> tt=0:0.01*pi:pi; yy=example2_e(tt);
>> [t,y]=RK4('example2_f',[0 pi],1,0.05*pi);
step
t
y
1
0.0000000000
1.0000000000
2
0.1570796327
0.8663284784
3
0.3141592654
0.7745866433
4
0.4712388980
0.7178375776
5
0.6283185307
0.6896194725
6
0.7853981634
0.6839104249
7
0.9424777961
0.6951106492
dy
8
1.0995574288
0.7180384347
y
9
1.2566370614
0.7479364401
dt
10
1.4137166941
0.7804851708
11
1.5707963268
0.8118207434
y( 0 ) 1
12
1.7278759595
0.8385543106
13
1.8849555922
0.8577907953
14
2.0420352248
0.8671448731
15
2.1991148575
0.8647524420
16
2.3561944902
0.8492761286
17
2.5132741229
0.8199036965
18
2.6703537556
0.7763385417
19
2.8274333882
0.7187817811
20
2.9845130209
0.6479057500
21
3.1415926536
0.5648190301
>> H=plot(x,y,'r-o',xx,yy);
>> set(H,'LineWidth',3,'MarkerSize',12);
sin( t )
Fourth-Order Runge-Kutta Method
dy
y sin( t ) ; y(0 ) 1
dt
Numerical Accuracy
» tt=0:0.01*pi:pi; yy=example2_e(tt);
» t0=0; y0=example2_e(t0);
» t1=0.1*pi; y1=example2_e(t1);
dy
y sin( t ) ; y(0 ) 1
dt
» [t,ya]=Eulode('example2_f',[0 pi],y0,0.1*pi);
» [t,yb]=midpoint('example2_f',[0 pi],y0,y1,0.1*pi);
h = 0.1
» [t,yc]=Heun_iter('example2_f',[0 pi],y0,0.1*pi,5);
» [t,yd]=RK4('example2_f',[0 pi],y0,0.1*pi);
» H=plot(t,ya,'m-*',t,yb,'c-d',t,yc,'g-s',t,yd,'r-O',tt,yy);
» set(H,'LineWidth',3,'MarkerSize',12);
»
» [t,ya]=Eulode('example2_f',[0 pi],y0,0.05*pi);
» [t,yb]=midpoint('example2_f',[0 pi],y0,y1,0.05*pi);
» [t,yc]=Heun_iter('example2_f',[0 pi],y0,0.05*pi,5);
h = 0.05
» [t,yd]=RK4('example2_f',[0 pi],y0,0.05*pi);
» H=plot(t,ya,'m-*',t,yb,'c-d',t,yc,'g-s',t,yd,'r-O',tt,yy);
» set(H,'LineWidth',3,'MarkerSize',12);
» print -djpeg075 ode06.jpg
Numerical Accuracy
Euler
Midpoint
Heun (iterative)
RK4
dy
y sin( t ) ; y(0 ) 1
dt
h = 0.1
Numerical Accuracy
Euler
Midpoint
Heun (iterative)
RK4
dy
y sin( t ) ; y(0 ) 1
dt
h = 0.05
Butcher’s sixth-order
Runge-Kutta Method
1
yi 1 yi (7k1 32k 2 12k 4 32k5 7k6 )
90
k1 f (ti , yi )
1
1
k 2 f (ti h, yi k1h )
4
4
1
1
1
k3 f (ti h, yi k1h k 2 h )
4
8
8
1
1
k 4 f ( t i 2 h, yi 2 k 2 h k 3h )
k f ( t 3 h, y 3 k h 9 k h )
i
i
1
4
5
4
16
8
k6 f (ti h, yi 3 k1h 2 k 2h 12 k3h 12 k 4h 8 k5h )
7
7
7
7
7
System of ODEs
A system of simultaneous ODEs
n equations with n initial conditions
dy1
dt
dy2
dt
dyn
dt
f 1 ( t , y1 , y 2 , , y n )
f 2 ( t , y1 , y 2 , , y n )
f n ( t , y1 , y 2 , , y n )
System of ODEs
Bungee Jumper’s velocity and position
Two simultaneous ODEs
dx
dt v
dv g c d v 2
dt
m
x (0 ) v (0 ) 0
gc d
gm
tanh
t
v ( t )
cd
m
gc d
m
x ( t ) c ln cosh m t
d
Second-Order ODE
d 2 y
dy
dt 2 g ( t , y, dt )
y( t ) , dy (t )
0
0
1
0
dt
Convert to two first-order ODEs
dy1 dy
y2
y1 y
dt
dt
let
dy
2
y
dy
d
y
2 dt
2
g ( t , y1 , y 2 )
2
dt
dt
y1 ( t 0 ) 0
I .C . s
y2 ( t 0 ) 1
System of Two first-order ODEs
Euler’s Method
Any method considered earlier can be used
Euler’s method for two ODE-IVPs
Basic Euler method
yi 1 yi hf ( t i , yi )
Two ODE-IVPs
y1, i 1 y1, i hf1 ( t i , y1, i , y2 , i )
y2 , i 1 y2 , i hf 2 ( t i , y1, i , y2 , i )
Hand Calculations: Euler’s Method
Solve the following ODE from t = 0 to t = 1 with h = 0.5
dy1
dt 0.5 y1
y1 (0 ) 4
dy
2 4 0.1 y 0.3 y y2 (0 ) 6
1
2
dt
Euler method
y1, i 1 y1, i f 1 ( t i , y1, i , y2 , i ) h y1, i (0.5 y1, i ) h
y2 , i 1 y2 , i f 2 ( t i , y1, i , y2 , i ) h y2 , i ( 4 0.1 y1, i 0.3 y2 , i ) h
y1 (0.5 ) y1 (0 ) 0.5y1 (0 )h 4 ( 0.5)(4) (0.5) 3
y2 (0.5 ) y2 (0 ) 4 0.1y1 (0 ) 0.3 y2 (0 )h 6 4 0.1( 4 ) 0.3(6 )(0.5 ) 6.9
y1 ( 1.0 ) y1 (0.5 ) 0.5y1 (0.5 )h 3 ( 0.5 )( 3 )(0.5 ) 2.25
y2 ( 1.0 ) y2 (0.5 ) 4 0.1y1 (0.5 ) 0.3 y2 (0.5 )h 6.9 4 0.1( 3 ) 0.3(6.9 )(0.5 ) 7.715
Euler’s Method for a System of ODEs
y is a column vector with n variables
Euler Method for a System of ODEs
function f = example5(t,y)
% dy1/dt = f1 = -0.5 y1
% dy2/dt = f2 = 4 - 0.1*y1 - 0.3*y2
% let y(1) = y1, y(2) = y2
% tspan = [0 1]
% initial conditions y0 = [4, 6]
f1 = -0.5*y(1);
f2 = 4 - 0.1*y(1) - 0.3*y(2);
f = [f1, f2]';
>> [t,y]=Euler_sys('example5',[0 1],[4 6],0.5);
t
y1
y2
y3 ...
0.000
4.0000000000
6.0000000000
0.500
3.0000000000
6.9000000000
(h = 0.5)
1.000
2.2500000000
7.7150000000
>> [t,y]=Euler_sys('example5',[0 1],[4 6],0.2);
t
y1
y2
0.000
4.0000000000
6.0000000000
0.200
3.6000000000
6.3600000000
0.400
3.2400000000
6.7064000000
0.600
2.9160000000
7.0392160000
(h
0.800
2.6244000000
7.3585430400
1.000
2.3619600000
7.6645424576
y3
...
= 0.2)
Euler Method for Second-Order ODE
d2y
dy
0.3
sin y
2
dt
dt
dy1
y
y
1
dt
let
dy
y2 dt
dy2
dt
Nonlinear Pendulum
dy
y2
dt
d2y
dy
0
.
3
sin y 0.3 y2 sin y1
2
dt
dt
function f = pendulum(t,y)
% nonlinear pendulum d^2y/dt^2 + 0.3dy/dt = -sin(y)
% convert to two first-order ODEs
% dy1/dt = f1 = y2
% dy2/dt = f2 = -0.1*y2 - sin(y1)
% let y(1) = y1, y(2) = y2
% tspan = [0 15]
% initial conditions y0 = [pi/2, 0]
f1 = y(2);
f2 = -0.3*y(2) - sin(y(1));
f = [f1, f2]';
Euler’s Method for Second-Order ODE
Nonlinear Pendulum
» [t,y1]=Euler_sys('pendulum',[0 15],[pi/2 0],15/100);
» [t,y2]=Euler_sys('pendulum',[0 15],[pi/2 0],15/200);
» [t,y3]=Euler_sys('pendulum',[0 15],[pi/2 0],15/500);
» [t,y4]=Euler_sys('pendulum',[0 15],[pi/2 0],15/1000);
» H=plot(t1,y1(:,1),t2,y2(:,1),t3,y3(:,1),t4,y4(:,1))
n = 100
n = 200
n = 500
n = 1000
Higher Order ODEs
In general, nth-order ODE
( n)
( n 1 )
)
y f ( t , y, y, y, , y
( n 1 )
y
(
t
)
,
y
(
t
)
,
,
y
( t 0 ) n 1
0
0
1
0
y1
y2
let y3
yn
y
y
y
y ( n 1 )
y1
y
2
y3
yn
y2 ,
y1 ( t 0 ) α 0
y3 ,
y2 ( t 0 ) α 1
y4 ,
y3 ( t 0 ) α 2
f ( t, y1 , y2 , , yn ) , yn ( t0 ) αn 1
System of First-Order ODE-IVPs
Example
y f ( t , y, y, y) t 2 4 ty 3 y 5 y
y(0 ) 2, y(0 ) 4 , y(0 ) 1
Convert to three first-order ODE-IVPs
let y1 y, y2 y dy/dt, y3 y d 2 y / dt 2
y1 f 1 ( t , y1 , y2 , y3 ) y2
y2 f 2 ( t , y1 , y2 , y3 ) y3
2
y
f
(
t
,
y
,
y
,
y
)
t
4 ty1 3 y2 5 y3
3
1
2
3
3
y1 (0 ) 2
y2 ( 0 ) 4
y (0 ) 1
3
In Vector Notation y f ( t , y )
Euler’s Method for Systems of
First-Order ODEs
Euler’s Method
y1 ( i 1) y1 ( i ) f1 ( t ( i ), y1 ( i ), y2 ( i ), y3 ( i )) h
y2 ( i 1) y2 ( i ) f 2 ( t ( i ), y1 ( i ), y2 ( i ), y3 ( i )) h
y ( i 1) y ( i ) f ( t ( i ), y ( i ), y ( i ), y ( i )) h
3
3
1
2
3
3
Example
y f ( t , y, y, y) t 2 4 ty 3 y 5 y
y(0 ) 2, y(0 ) 4 , y(0 ) 1
y1 y2
y2 y3
2
y
t
4 ty1 3 y2 5 y3
3
y1 (0 ) 2
y2 ( 0 ) 4
y (0 ) 1
3
Example: Euler’s Method
First step: t(0) = 0, t(1) = 0.5 (h = 0.5)
y1 ( 1) y1 (0 ) f1 (0 )h y1 (0 ) y2 (0 ) h 2 ( 4 )( 0.5 ) 4.0
y2 ( 1) y2 (0 ) f 2 (0 )h y2 (0 ) y3 (0 )h 4 ( 1)( 0.5 ) 4.5
2
y
(
1
)
y
(
0
)
f
(
0
)
h
y
(
0
)
0
4 (0 ) y1 (0 ) 3 y2 (0 ) 5 y3 (0 )
3
3
3
3
2
1
(
0
)
4 (0 )( 2 ) 3( 4 ) 5 ( 1) (0.5 ) 2.5
Second step: t(1) = 0.5, t(2) = 1.0
y1 ( 2 ) y1 ( 1) f1 ( 1)h y1 ( 1) y2 ( 1)h 4.0 ( 4.5 )( 0.5 ) 6.25
y2 ( 2 ) y2 ( 1) f 2 ( 1)h y2 ( 1) y3 ( 1)h 4.5 ( 2.5 )( 0.5 ) 3.25
2
y
(
2
)
y
(
1
)
f
(
1
)
h
y
(
1
)
(
0
.
5
)
4 (0.5 ) y1 ( 1) 3 y2 ( 1) 5 y3 ( 1) h
3
3
3
3
2
2
.
5
(
0
.
5
)
4 (0.5 )( 4.0 ) 3( 4.5 ) 5 ( 2.5 ) (0.5 ) 11.375
Classical 4th-order Runge-Kutta Method
for Systems of ODE-IVPs
k1,1 f1 (t (i ), y1 (i ), y2 (i ))
k1,2 f 2 (t (i ), y1 (i ), y2 (i ))
2 equations
k 2,1 f1 (t(i ) h/2, y1 (i ) k1,1h/2, y 2 (i ) k 1,2h/2) f1(t * ( i ),y1* ( i ),y*2 (i ))
*
*
*
k 2,2 f 2 (t(i ) h/2, y1 (i ) k1,1h/2, y 2 (i ) k 1,2h/2) f 2 (t ( i ),y1 ( i ),y 2 (i ))
k 3,1 f1 (t(i ) h/2, y1 (i ) k 2,1h/2, y 2 (i ) k 2,2h/2) f1(t * ( i ),y1** ( i ),y**
2 (i ))
*
**
**
k 3,2 f 2 (t(i ) h/2, y1 (i ) k 2,1h/2, y 2 (i ) k 2,2 h/2) f 2 (t (i ),y1 ( i ),y 2 (i ))
k 4,1 f1 (t(i ) h, y1 (i ) k 3,1h, y 2 (i ) k 3,2 h) f1 (t ** (i ),y1*** (i ),y ***
2 ( i ))
**
***
***
k 4,2 f 2 (t(i ) h, y1 (i ) k 3,1h, y 2 (i ) k 3,2h) f 2 (t ( i ),y1 ( i ),y 2 ( i ))
1
y
(
i
1
)
y
(
i
)
( k1, 1 2 k 2 , 1 2 k 3 , 1 k 4 , 1 ) h
1
1
6
y ( i 1) y ( i ) 1 ( k 2 k 2 k k )h
2
1, 2
2, 2
3, 2
4,2
2
6
Applicable for
any number of
equations
Hand Calculations: RK4 Method
Solve the following ODE from t = 0 to t = 1 with h = 0.5
dy1
dt f1 0.5 y1
y1 (0 ) 4
dy
2 f 4 0.1 y 0.3 y y2 (0 ) 6
2
1
2
dt
Classical RK4 method
k1, 1 f 1 (0 , y1 (0 ), y2 (0 )) f 1 (0 , 4 ,6 ) (0.5 )( 4 ) 2
k1, 2 f 2 (0 , y1 (0 ), y2 (0 )) f 2 (0 , 4 ,6 ) 4 (0.1)( 4 ) (0.3 )( 6 ) 1.8
y1* y1 (0 ) k1, 1 h / 2 4 ( 2 )( 0.5 ) / 2 3.5
*
y2 y2 (0 ) k1, 2 h / 2 6 ( 1.8 )( 0.5 ) / 2 6.45
k 2 , 1 f 1 (0.25 , y1* , y2* ) f 1 (0.25 , 3.5 ,6.45 ) (0.5 )( 3.5 ) 1.75
*
*
k 2 , 2 f 2 (0.25 , y1 , y2 ) f 2 (0.25 , 3.5 ,6.45 ) 4 (0.1)( 3.5 ) 0.3(6.45 ) 1.715
Continued: RK4 Method
y1** y1 (0 ) k2 , 1 h / 2 4 ( 1.75 )( 0.5 ) / 2 3.5625
**
y2 y2 (0 ) k2 , 2 h / 2 6 ( 1.715 )( 0.5 ) / 2 6.42875
k3 , 1 f1 (0.25 , y1** , y2** ) f 1 (0.25 , 3.5625 ,6.42875 ) (0.5 )( 3.5625 ) 1.78125
**
**
k3 , 2 f 2 (0.25 , y1 , y2 ) f 2 (0.25 , 3.5625 ,6.42875 ) 4 (0.1)( 3.5625 ) 0.3(6.42875 ) 1.715125
y1*** y1 (0 ) k3 , 1 h 4 ( 1.78125 )( 0.5 ) 3.109375
***
y2 y2 (0 ) k3 , 2 h 6 ( 1.715125 )( 0.5 ) 6.857563
k4 , 1 f1 (0.5 , y1*** , y2*** ) f1 (0.5 , 3.109375 ,6.857563 ) (0.5 )( 3.109375 ) 1.554688
***
***
k4 , 2 f 2 (0.25 , y1 , y2 ) f 2 (0.5 , 3.109375 ,6.857563 ) 4 (0.1)( 3.109375 ) 0.3(6.857563 ) 1.631794
y1 (0.5 ) y1 (0 ) ( 1 / 6 )( k1, 1 2 k2 , 1 2 k3 , 1 k4 , 1 )h
4 ( 1 / 6 )[( 2 ) 2( 1.75 ) 2( 1.78125 ) 1.554688]( 0.5 ) 3.115234
y2 (0.5 ) y2 (0 ) ( 1 / 6 )( k1, 2 2 k2 , 2 2 k3 , 2 k4 , 2 )h
6 ( 1 / 6 )[( 1.8 ) 2( 1.715 ) 2( 1.715125 ) ( 1.631794 )]( 0.5 ) 6.857670
4th-order Runge-Kutta Method for ODEs
Valid for any number of coupled ODEs
4th-order Runge-Kutta Method for ODEs
>> [t,y]=RK4_sys('example5',[0 10],[4 6],0.5);
t
y1
y2
0.000
4.0000000000
6.0000000000
0.500
3.1152343750
6.8576703125
1.000
2.4261713028
7.6321056734
1.500
1.8895230605
8.3268859767
2.000
1.4715767976
8.9468651000
2.500
1.1460766564
9.4976013588
3.000
0.8925743491
9.9849540205
3.500
0.6951445736
10.4148035640
4.000
0.5413845678
10.7928635095
4.500
0.4216349539
11.1245594257
5.000
0.3283729256
11.4149566980
5.500
0.2557396564
11.6687232060
6.000
0.1991722422
11.8901165525
6.500
0.1551170538
12.0829881442
7.000
0.1208064946
12.2507984405
7.500
0.0940851361
12.3966392221
8.000
0.0732743126
12.5232598757
8.500
0.0570666643
12.6330955637
9.000
0.0444440086
12.7282957874
9.500
0.0346133758
12.8107523359
10.000
0.0269571946
12.8821259602
y3
...
4th-order Runge-Kutta Method for ODE-IVPs
Nonlinear Pendulum
»
»
»
»
»
»
tspan=[0 15]; y0=[pi/2 0];
[t1,y1]=RK4_sys(‘pendulum',tspan,y0,15/25);
[t2,y2]=RK2_sys(‘pendulum',tspan,y0,15/50);
[t3,y3]=RK2_sys(‘pendulum',tspan,y0,15/100);
H=plot(t1,y1(:,1),t2,y2(:,1),t3,y3(:,1));
set(H,'LineWidth',3.0);
n = 25
n = 50
n = 100
Comparison of Numerical Accuracy
Nonlinear Pendulum
» tspan=[0 15]; y0=[pi/2 0];
» [t1,y1]=Euler_sys(‘pendulum',tspan,y0,15/100);
» [t2,y2]=RK2_sys(‘pendulum',tspan,y0,15/100);
» [t3,y3]=RK4_sys(‘pendulum',tspan,y0,15/100);
» H1=plot(t1,y1(:,1),t2,y2(:,1),t3,y3(:,1)); hold on;
» H2=plot(t1,y1(:,2),'b:',t2,y2(:,2),'g:',t3,y3(:,2),'r:');
Euler
RK2
RK4
Example: More than 2 ODE-IVPs
dy1
dt
dy2
dt
dy3
dt
dy
4
dt
f 1 36 y1 30 y2 20 y3 10 y4 ; y1 (0 ) 1
f 2 61 y1 50 y2 36 y3 18 y4 ; y2 (0 ) 0
f 3 34 y1 29 y2 25 y3 13 y4 ; y3 (0 ) 0
f 4 10 y1 10 y2 10 y3 6 y4 ;
36 30 20 10
61 50 36 18
A
34 29 25 13
10
10
10
6
t 0 1 0 0 0
y4 ( 0 ) 0
function f = example(t, y)
% solve y' = Ay = f, y0 = [1 0 0 0]'
% four first-order ODE-IVPs
A = [ -36 30 -20 10
-61 50 -36 18
-34 29 -25 13
-10 10 -10
6];
y=[ y(1) y(2) y(3) y(4)]';
f = A*y;
4th-order Runge-Kutta Method for ODE-IVPs
Symbols: n = 20
Lines : n =100
» tspan=[0 2]; y0=[1 0 0 0];
» [t1,y1]=RK4_sys(‘example', tspan, y0, 2/20);
» [t2,y2]=RK4_sys(‘example', tspan, y0, 2/100);
» H1=plot(t1,y1,'o'); set(H1,'LineWidth',3','MarkerSize',12);
» hold on; H2=plot(t2,y2); set(H2,'LineWidth',3);
CVEN 302-501
Homework No. 13
Chapter 20
Prob. 20.1 (40), (Hand Calculation, and use
MATLAB plotting for graph)
20.3 a) , b) and c) (40)(Hand Calculation)
Prob. 20.8 (30) (decomposing into two 1st
ODEs and then using MATLAB Program)
Due on Wed. 11/26/2008 at the beginning of
the period