ODEs:

1
Chapter 20

• ODEs
• Initial-Value Problems

2
Differential Equations

• There are ordinary differential equations -
  functions of one variable
• And there are partial differential equations -
  functions of multiple variables

3
Order of differential equations

• 1st order (falling parachutist)
• 2nd order (mass-spring system with damping)
• etc.

4
Higher-Order ODEs

• Can always turn a higher order ODE into a set of
  1st order ODEs
• Example
• Let then
• So solutions to first-order ODEs are important

5
Linear and Nonlinear ODEs

• Linear No multiplicative mixing of variables, no
  nonlinear functions
• Nonlinear anything else

6
Ordinary Differential Equations

• ODEs show up everywhere in engineering
• Dynamics (Newtons 2nd law)
• Heat conduction (Fouriers law)
• Diffusion (Ficks law)

7

• Example of an ODE

x
k
f(t)
What order is this ODE? If f(t) 0, ODE is
homogenous. If f(t) is not equal to 0, ODE is
non-homogenous.
m
c
kx
Free-body diagram
f(t)
m
8
Solutions for ODEs
The solution for the homogenous ODE.
The solution for the non-homogenous ODE
The arbitrary constants C1 C2 are determined by
the Initial-value or Boundary-value conditions.
9
Initial-Value Boundary-Value Conditions
The I-V Conditions All conditions are
given at the same value of the independent
variable.
The B-V Conditions Conditions are
given at the different values of the independent
variable.
The numerical schemes for solving Initial-value
and boundary-value are different.
10
Ordinary Differential Equations I-V Problems

• Eulers and Heun's methods
• Runge-Kutta methods
• Adaptive Runge-Kutta
• Multistep methods
• Adams-Bashforth-Moulton methods

11
Ordinary Differential Equations

• 1st order Ordinary differential equations (ODEs)
• Initial value problems
• Numerical approximations
• New value old value slope step size

12
Runge-Kutta Methods

• Runge Kutta methods (One Step Methods)
• Idea is that
• New value old value slopestep size
• Slope is generally a function of t, hence y(t)
• Different methods differ in how to estimate ?

13
One-Step Method
All one-step methods can be expressed in this
general form, the only difference being the
manner in which the slope is estimated
14
Initial Conditions
The solution of ODE depends on the initial
condition
Same ODE, but with different initial conditions
t
15
Eulers method(First-order Taylor Series Method)

• Approximate the derivative by finite difference
• Local truncation error
• or

16
Eulers Method

• Eulers (Euler-Cauchy or point-slope) method
• Use the slope at ti to predict yi1

17
Eulers method
Straight line approximation
y0
t0
t1
t2
t3
h
h
h
18
Example Eulers Method

• Analytic solution
• Euler method
• h 0.50

19
Example Eulers Method

• h 0.25

20
Eulers method
y
h 0.25
1
h 0.5
t
0.5
0.75
1.0
0.25
0
21
Eulers Method (modified M-file)
22
Eulers Method
gtgt tt00.01pipi gtgt ye1.5exp(-tt)0.5sin(tt)
-0.5cos(tt) gtgt t,yEulode('example2_f',0
pi,1,0.1pi) step t y
1 0.0000000000 1.0000000000 2
0.3141592654 0.6858407346 3
0.6283185307 0.5674580652 4
0.9424777961 0.5738440394 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 gtgt
Hplot(t,y,'r-o',tt,ye) gtgt set(H,'LineWidth',3,'M
arkerSize',12) gtgt print -djpeg ode01.jpg
23
Eulers Method (h 0.1?)
24
Eulers Method (h 0.05?)
25
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

26
Global Local Errors
Global error
Local error
y
y
yi1
yi1
yi
yi
x
o
xi
o
xi1
xi1
xi2
x
xi
27
Eulers Method

• Eulers method uses Taylor series with only first
  order terms
• The true local truncation error is
• Approximate local truncation error - neglect
  higher order terms (for sufficiently small h)

28
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
• Heuns method
• Midpoint (or improved polygon) method
• Third-order Runge-Kutta methods
• Fourth-order Runge-Kutta methods

29
Heuns Method

• Improvements of Eulers method - Heuns method
• Eulers method - derivative at the beginning of
  interval is applied to the entire interval
• Heuns method uses average derivative for the
  entire interval
• A second-order Runge-Kutta Method

30
Heuns Method

• Heuns method is a predictor-corrector method
• Predictor
• Corrector (may be applied iteratively)

31
Heuns Method

• Iterate the corrector of Heuns method to obtain
  an improved estimate

Predictor
Corrector
32
Heuns Method with Iterative Correctors
33
Heuns Method with Iterative Correctors
te00.02pipi yeexample2_e(xe)
t1,y1Euler('example2_f',0 pi,1,0.1pi)
t2,y2Heun_iter('example2_f',0
pi,1,0.1pi,0) t3,y3Heun_iter('example2_f',
0 pi,1,0.1pi,5) Hplot(te,ye,t1,y1,'r-d',t2,
y2,'g-s',t3,y3,'m-o') set(H,'LineWidth',3,'Mark
erSize',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
34
Heuns Method with Iterative Correctors
Heuns with 5 iterations
Heuns method
Exact
Eulers
t
35
Example Heuns Method

• h 0.5
• First Step

36
Example Heuns Method

• h 0.5
• Second Step

37
Example
t1,y1Eulode('example3',0 5,1,0.5)
t2,y2Heun_iter('example3',0,5,1,0.5,0)
t3,y3Heun_iter('example3',0 5,1,0.5,5)
t1' y1' y2' y3' ye' t yEuler
yHeun yHeun_iter ytrue 0 1.0000
1.0000 1.0000 1.0000 0.5000 1.0000
1.1250 1.1331 1.1289 1.0000 1.2500
1.5523 1.5804 1.5625 1.5000
1.8090 2.4169 2.4859 2.4414 2.0000
2.8178 3.9463 4.0882 4.0000 2.5000
4.4964 6.4619 6.7192 6.5664 3.0000
7.1470 10.3793 10.8046 10.5625
3.5000 11.1570 16.2082 16.8630 16.5039
4.0000 17.0024 24.5532 25.5065 25.0000
4.5000 25.2492 36.1126 37.4407 36.7539
5.0000 36.5552 51.6796 53.4643 52.5625
38
Midpoint Method

• Improved Polygon or Modified Euler Method
• Use the slope at midpoint to represent the
  average slope

39
Runge-Kutta (RK) Methods

• One-Step Method with general form
• Where ? is an increment function which represents
  the weighted-average slope over the interval
• Where as are constants and ks are slopes
  evaluated at selected x locations

40
Runge-Kutta (RK) Methods

• One-Step Method
• General Form of nth-order Runge-Kutta Method
• Where ps and qs are constants
• ks are recurrence relationships

41
Second-Order RK Methods

• Taylor series expansion
• Compare to the second-order Taylor formula
• Three equations for four unknowns (a1, a2, k1,
  q11)

42
Second-Order RK Methods

• Second-order version of Runge-Kutta Methods
• 3 equations for 4 unknowns

43
Second-Order RK Methods

• General Second-order Runge-Kutta methods

k2
a1 k1 a2 k2
k1
Weighted-average slope
xi1 xih
xi
xip1h
44
Heuns Method

• Heuns 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

Use average slope over the interval
45
Midpoint Method

• Another Second-order Runge-Kutta method
• Choose a2 1 ? a1 0, p1 1/2, q11 1/2

k2
k1
xi
xi1/2
xi1
46
Midpoint (2nd-order RK) Method
47
Midpoint (2nd-order RK) Method
gtgt t,y midpoint('example2_f',0
pi,1,0.05pi) 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 10
1.4137166941 0.7833740988 11
1.5707963268 0.8143967658 12
1.7278759595 0.8407772414 13
1.8849555922 0.8596353281 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 gtgt tt
00.01pipi yy example2_e(tt) gtgt H
plot(t,y,'r-o',tt,yy) gtgt set(H,'LineWidth',3,'Mar
kerSize',12)
48
Midpoint (RK2) Method
49
Ralstons Method

• Second-order Runge-Kutta method
• Choose a2 2/3 ? a1 1/3, p1 q11 3/4

k2
k1
k1/3 2k2/3
xi1 tih
ti
ti3h/4
50
Third-Order Runge-Kutta Method

• General form
• Weighted slope

51
Third-Order Runge-Kutta Methods

• General Third-order Runge-Kutta methods

k3
k2
Weighted-average value of three slopes k1 , k2 ,
k3
k1
ti
tip1h
ti1 tih
tip2h
52
Third-order Runge-Kutta Methods
Nystrom Method Nearly Optimum Method
53
3rd-order Runge-Kutta Method

• Reduce to Simpsons 1/3 rule for f f(t)

54
3rd-Order Heun Method
55
Classical 4th-orderRunge-Kutta Method

• One-step method
• Reduce to Simpsons 1/3 rule for f f(t)

56
Classical 4th-orderRunge-Kutta Method
k2
k4
k3
k1
ti
ti h/2
ti h
57
Example Classical 4th-order RK Method

• h 0.5
• First step, t 0.5

58
Example Classical 4th-order RK Method

• Second step, t 1.0

59
Fourth-Order Runge-Kutta Method
60
Fourth-Order Runge-Kutta Method
gtgt tt00.01pipi yyexample2_e(tt) gtgt
t,yRK4('example2_f',0 pi,1,0.05pi) 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 8 1.0995574288
0.7180384347 9 1.2566370614
0.7479364401 10 1.4137166941
0.7804851708 11 1.5707963268
0.8118207434 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 gtgt Hplot(x,y,'r-o',xx,yy) gtgt
set(H,'LineWidth',3,'MarkerSize',12)
61
Fourth-Order Runge-Kutta Method
62
Numerical Accuracy
tt00.01pipi yyexample2_e(tt) t00
y0example2_e(t0) t10.1pi
y1example2_e(t1) t,yaEulode('example2_f',0
pi,y0,0.1pi) t,ybmidpoint('example2_f',0
pi,y0,y1,0.1pi) t,ycHeun_iter('example2_f
',0 pi,y0,0.1pi,5) t,ydRK4('example2_f',
0 pi,y0,0.1pi) Hplot(t,ya,'m-',t,yb,'c-d',t
,yc,'g-s',t,yd,'r-O',tt,yy) set(H,'LineWidth',3
,'MarkerSize',12) t,yaEulode('example2_f'
,0 pi,y0,0.05pi) t,ybmidpoint('example2_f
',0 pi,y0,y1,0.05pi) t,ycHeun_iter('examp
le2_f',0 pi,y0,0.05pi,5) t,ydRK4('example
2_f',0 pi,y0,0.05pi) Hplot(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
h 0.1?
h 0.05?
63
Numerical Accuracy
Euler Midpoint Heun (iterative) RK4
h 0.1?
64
Numerical Accuracy
Euler Midpoint Heun (iterative) RK4
h 0.05?
65
Butchers sixth-order Runge-Kutta Method
66
System of ODEs

• A system of simultaneous ODEs
• n equations with n initial conditions

67
System of ODEs

• Bungee Jumpers velocity and position
• Two simultaneous ODEs

68
Second-Order ODE

• Convert to two first-order ODEs

69
System of Two first-order ODEsEulers Method

• Any method considered earlier can be used
• Eulers method for two ODE-IVPs
• Basic Euler method
• Two ODE-IVPs

70
Hand Calculations Eulers Method

• Solve the following ODE from t 0 to t 1 with
  h 0.5
• Euler method

71
Eulers Method for a System of ODEs
y is a column vector with n variables
72
Euler Method for a System of ODEs
function f example5(t,y) dy1/dt f1 -0.5
y1 dy2/dt f2 4 - 0.1y1 - 0.3y2 let y(1)
y1, y(2) y2 tspan 0 1 initial
conditions y0 4, 6 f1 -0.5y(1) f2 4 -
0.1y(1) - 0.3y(2) f f1, f2'
gtgt t,yEuler_sys('example5',0 1,4 6,0.5)
t y1 y2 y3
... 0.000 4.0000000000 6.0000000000
0.500 3.0000000000 6.9000000000 1.000
2.2500000000 7.7150000000 gtgt
t,yEuler_sys('example5',0 1,4 6,0.2)
t y1 y2 y3
... 0.000 4.0000000000 6.0000000000
0.200 3.6000000000 6.3600000000 0.400
3.2400000000 6.7064000000 0.600
2.9160000000 7.0392160000 0.800
2.6244000000 7.3585430400 1.000
2.3619600000 7.6645424576
(h 0.5)
(h 0.2)
73
Euler Method for Second-Order ODE
Nonlinear Pendulum
function f pendulum(t,y) nonlinear pendulum
d2y/dt2 0.3dy/dt -sin(y) convert to two
first-order ODEs dy1/dt f1 y2 dy2/dt f2
-0.1y2 - sin(y1) let y(1) y1, y(2) y2
tspan 0 15 initial conditions y0 pi/2,
0 f1 y(2) f2 -0.3y(2) - sin(y(1)) f
f1, f2'
74
Eulers Method for Second-Order ODE
Nonlinear Pendulum
t,y1Euler_sys('pendulum',0 15,pi/2
0,15/100) t,y2Euler_sys('pendulum',0
15,pi/2 0,15/200) t,y3Euler_sys('pendulum
',0 15,pi/2 0,15/500) t,y4Euler_sys('pen
dulum',0 15,pi/2 0,15/1000)
Hplot(t1,y1(,1),t2,y2(,1),t3,y3(,1),t4,y4(,1)
)
n 100 n 200 n 500 n 1000
75
Higher Order ODEs

• In general, nth-order ODE

76
System of First-Order ODE-IVPs

• Example
• Convert to three first-order ODE-IVPs

77
Eulers Method for Systems of First-Order ODEs

• Eulers Method
• Example

78
Example Eulers Method

• First step t(0) 0, t(1) 0.5 (h 0.5)
• Second step t(1) 0.5, t(2) 1.0

79
Classical 4th-order Runge-Kutta Method for
Systems of ODE-IVPs
2 equations
Applicable for any number of equations
80
Hand Calculations RK4 Method

• Solve the following ODE from t 0 to t 1 with
  h 0.5
• Classical RK4 method

81
Continued RK4 Method
82
4th-order Runge-Kutta Method for ODEs
Valid for any number of coupled ODEs
83
4th-order Runge-Kutta Method for ODEs
gtgt t,yRK4_sys('example5',0 10,4 6,0.5)
t y1 y2 y3
... 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
84
4th-order Runge-Kutta Method for ODE-IVPs
Nonlinear Pendulum
tspan0 15 y0pi/2 0 t1,y1RK4_sys(pe
ndulum',tspan,y0,15/25) t2,y2RK2_sys(pendul
um',tspan,y0,15/50) t3,y3RK2_sys(pendulum',
tspan,y0,15/100) Hplot(t1,y1(,1),t2,y2(,1),t
3,y3(,1)) set(H,'LineWidth',3.0)
n 25 n 50 n 100
85
Comparison of Numerical Accuracy
Nonlinear Pendulum
tspan0 15 y0pi/2 0 t1,y1Euler_sys(
pendulum',tspan,y0,15/100) t2,y2RK2_sys(pen
dulum',tspan,y0,15/100) t3,y3RK4_sys(pendul
um',tspan,y0,15/100) H1plot(t1,y1(,1),t2,y2(
,1),t3,y3(,1)) hold on H2plot(t1,y1(,2),'b
',t2,y2(,2),'g',t3,y3(,2),'r')
?
?
Euler RK2 RK4
??
??
86
Example More than 2 ODE-IVPs
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 Ay
87
4th-order Runge-Kutta Method for ODE-IVPs
Symbols n 20 Lines n 100
tspan0 2 y01 0 0 0 t1,y1RK4_sys(ex
ample', tspan, y0, 2/20) t2,y2RK4_sys(examp
le', tspan, y0, 2/100) H1plot(t1,y1,'o')
set(H1,'LineWidth',3','MarkerSize',12) hold
on H2plot(t2,y2) set(H2,'LineWidth',3)
88
CVEN 302-501Homework 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