ODEs: Adaptive Methods and Stiff Systems

1
Chapter 21

• ODEs Adaptive Methods and Stiff Systems

2
Adaptive Runge Kutta Method
Use small step size in high gradient region
(abrupt change) automatic step size adjustment
3
Adaptive Runge Kutta Method

• First approach Step halving
• Estimate local truncation error using two
  different step sizes
• Solve each step twice, once as a full step and
  then as two half steps
• Second approach Embedded RK methods
• (also called RK-Fehlberg methods)
• Estimate local truncation error between two
  predictions using two different-order RK methods

4
Step Halving Method

• Step Halving (Adaptive RK) method
• Compute the solutions at each step twice using
  4th-order classic RK method
• Once as a full step h and independently as two
  half steps (h/2)
• y1 one full step y2 two half steps

Error estimate Correction 5th-order
5
Adaptive 4th-order RK Method

• One full-step with h
• Two half-steps with h/2
• Therefore,

Two half-steps
6
Embedded Runge-Kutta Method

• MATLAB function ODE23
• BS23 algorithm (Bogacki and Shampine, 1989
  Shampine, 1994)
• Use 3rd-order 4th-order RK methods
  simultaneously to solve the ODE and estimate the
  error for step-size adjustment
• Error estimate (Note k1 is the same as k4 from
  previous step)

7
Embedded RK Method ODE23

• Uses only three function evaluations (k1, k2, k3)
• After each step, the error is checked to
  determine whether it is within desired tolerance.
  If it is, yi1 is accepted and k4 becomes k1 for
  the next time step
• If the error is too large, the step is repeated
  with a reduced step sizes until the estimate
  error is acceptable
• RelTol relative tolerance (default 10?3)
• AbsTol relative tolerance (default 10?6)

8
Adaptive RK Method ode23

• Example 21.2 Use ode23 to solve the following
  ODE from t 0 to 4

function yp ex21_2(t, y) Example 21.2 in the
text book yp 10exp(-(t-2)(t-2)/(20.0752))
- 0.6y
9
Example 20.2 ode23
gtgt options odeset('RelTol',1.e-4) gtgt
ode23('ex21_2', 0 4, 0.5, options)
gtgt options odeset('RelTol',1.e-3) gtgt
ode23('ex21_2', 0 4, 0.5, options)
(a) RelTol 10?3
(b) RelTol 10?4
10
Other MATLAB Functions

• MATLAB Function ode45
• Dormand and Prince (1990)
• Solve fourth- and fifth-order RK formulas
  simultaneously
• Make error estimates for step-size adjustment
• MATLAB Function ode113
• Adams-Bashforth-Moulton solver (order 1-12)
• Predictor-corrector method

11
Runge-Kutta Fehlberg Method

• Fourth-order
• Fifth-order

These coefficients were developed by Cash and
Karp (1990). Also called Cash-Karp RK Method
12
Runge-Kutta Fehlberg Method

• Identical coefficients (k1, k2, k3, k4, k5, k6)
  for both the fourth and fifth order
  Runge-Kutta-Fehlberg methods
• Save CPU time
• Error estimate use two RK methods of different
  order to estimate the local truncation error

13
Runge-Kutta Fehlberg Method

• First, calculate yi1 using 4th-order Runge-Kutta
  Felberg method ? (y1)4th
• Then, calculate yi1 using 5th-order Runge-Kutta
  Felberg method ? (y2)5th
• Calculate Error Ea (y2)5th - (y1)4th
• Adjust step size according to error estimate

14
Step Size Control

• Use step-halving or Runge-Kutta Fehlberg to
  estimate the local truncation error
• Adjust the step size according to error estimate
• Increase the step size if the error is too small
  and decrease it if the error is too large
• For fourth-order schemes

15
Step Size Adjustment
For nth-order RK Method

• For step size increases (RK4, n 4)
• For step size decreases, h is implicit in ?new

Reduce hnew also reduce ?new
16
Example21.2
Forcing function
Runge-Kutta-Fehlberg method with adaptive step
size control
t
small step size around t 2
solution
t
17
Predator-Prey Equation

• A simple predator-prey relationship is described
  by the Lotka-Volterra model, which we write in
  terms of a fox population f(t), with birth rate
  bf and death rate df and a geese population g(t)
  with birth rate bg and death rate dg

18
Predator-Prey Equation

• Example Given bf 0.3, df 0.8, bg 1.2, dg
  0.6, find the populations of predators and preys
  as a function of time (t 0 to 20) using ode45.

Predator Prey
function yp predprey(t,y) Fox (Predator)
population y1(t), birth rate bf, death rate df
Geese (prey) poupulation y2(t), birth rate bg,
death rate dg bf 0.3 df 0.8 bg 1.2 dg
0.6 yp y(1)(bfy(2)-df) y(2)(bg-dgy(1))
19
Adaptive RK method ode45
gtgt tspan0 20 y0 1, 2 gtgt t,y
ode45('predprey',tspan,y0) gtgt out t y out
0 1.00000000000000
2.00000000000000 0.08372954771699
0.98466335925883 2.10387606351447
0.16745909543397 0.97217906418030
2.21469691968449 0.25118864315096
0.96261463455261 2.33264806768044
0.33491819086794 0.95606000043195
2.45787517860514 0.62870201024789
0.95940830414765 2.95257534365817
0.92248582962784 1.00915715534212
3.53554512325941 1.21626964900779
1.11944001342999 4.17743733245175
1.51005346838774 1.31420023093313
4.80288705430266 1.70029730481714
1.49877617603179 5.14418160032910
1.89054114124654 1.73932623590652
5.37426899115905 2.08078497767594
2.03724630610217 5.44284022577897
2.27102881410534 2.38123550265509
5.31619117469119 2.46127265053474
2.74867929767025 4.98122010458805
2.65151648696414 3.09373931102751
4.48531246964154 2.84176032339354
3.37201992118282 3.89965751560944
3.03200415982294 3.55414367377695
3.29772307544844 3.21873982743099
3.62508046596480 2.75627545562992
3.40547549503903 3.59592495497541
2.29786856084044 3.59221116264708
3.48499741488926 1.93369109247072
3.77894683025512 3.31638942318525
1.65488750920928 3.94714350340851
3.13529220765148 1.46223282653043
4.11534017656189 2.93901369774356
1.31685686503153 4.28353684971528
2.73779832900552 1.20974621108398
4.45173352286867 2.53877488937042
1.13408335792525 4.65875081398499
2.30372716989349 1.07647855767312
4.86576810510130 2.08503830477357
1.05104469525303 5.07278539621762
1.88590370110986 1.05298367705029 ...
... ...
19.13258672508553 1.39621015109817
1.21506308763510 19.36078504002272
1.26894875046392 1.33168491347240
19.58898335495991 1.16392668931057
1.48284270918163 19.69173751621993
1.12358612912238 1.56325743063018
19.79449167747995 1.08747228241032
1.65192298666833 19.89724583873998
1.05554028549720 1.74928308667197
20.00000000000000 1.02776980338246
1.85579178601841
20
Predator-Prey Model

• Time history of predator (fox) and prey (geese)
  populations

21
Predator-Prey Model
State-space plot
22
Multistep Methods

• Runge-Kutta methods
• -- one-step method
• -- use intermediate values between ti and
  ti1
• -- several evaluations of slope per step
• Multistep methods
• -- use values at ti , ti-1 , ti-2 etc
• -- only one evaluation of derivative per step

23
One-Step and Multistep Methods
One-step
Multistep
24
Multi-step methods use several previous points
(yi-1 , yi-2 ,) in addition to the present point
yi

• - explicit (b0 0) implicit methods.
• - of the previous steps.
• Non-self start Huen method
• Adams-Bashforth Methods (explicit methods b0
  0)
• Admas-Moulton Methods (implicit methods)
• Predictor-Corrector Methods

25
Heuns Method (Review)

• Start with second-order method
• Look at Heuns method (Self-Starting)
• Euler predictor - O(h2)
• Trapezoid corrector - O(h3)

26
Non-Self-Starting Heuns Method

• To improve predictor, use 3rd-order Euler method
• along with the same old corrector
• iterated until converged

Note are the final results of the
corrector iterations at the previous time step
27
Non-Self Starting Heuns Method
(a) 3rd-order predictor (b) 3rd-order
corrector
Non-self starting, need yi-1 and yi (two initial
values)
28
Heuns Method
(a) 2nd-order predictor (b) 3rd-order
corrector
Self-starting, need yi only
29
Truncation Errors

• Non-self-starting Heun method
• Predictor truncation error
• Corrector truncation error

30
Modifiers

• Corrector Modifier
• Predictor Modifier (not used in textbook)

Error in textbook
31
Non-Self Starting Heuns Method

• So the sequence is
• Predict
• Adjust prediction
• Correct
• Converged?
• If not, correct again

32
Hand Calculations

• Non-self-starting Heuns method
• Need two initial conditions
• Example 21.3 (p523- 524)

Exact solution y 3.076923e 0.8t 1.076923 e
0.5 t
33
Hand Calculations First Step
Predictor Corrector

• First Step (i 0) need two initial conditions
• Predictor

34
Hand Calculations First Step

• Corrector

35
Hand Calculations Second Step

• Second Step (i 1)
• Predictor

36
Hand Calculations First Step

• Corrector

37
Stiffness

• A stiff system involves rapidly changing
  components (usually die away quickly) together
  with slowly ones
• Long-time solution is dominated by slow varying
  components

38
Numerical Stability

• Amplification or decay of numerical errors
• A numerical method is stable if error incurred at
  one stage of the process do not tend to magnify
  at later stages
• Ill-conditioned differential equation
• -- numerical errors will be magnified
  regardless
• of numerical method
• Stiff differential equation
• -- require extremely small step size to
  achieve
• accurate results

39
Stability

• Example problem

40
Euler Explicit Method

• Stability criterion
• Region of absolute stability

Amplification factor
41
Euler Implicit Method
Forward, Backward, or Centered Scheme?

• Unconditionally stable !

42
Stability

• Explicit Euler method
• Second-order Adams-Moulton
• Implicit methods are in general more stable than
  explicit methods.

If a 1000, then h ? 0.002 to ensure stability
Unconditionally stable
43
Example 21.5
Euler Explicit Euler Implicit
44
MATLAB Functions for Stiff Systems

• ode15s based on Gear backward differentiation
• for stiff problems of low to
  medium accuracy
• ode23s based on modified Rosenbrock formula
• ode23t trapezoidal rule with a free interpolant
• for moderately stiff problems
  without
• numerical damping
• ode23tb implicit Runge-Kutta formula
• for more explanation please see page 529

45
Stiff ODE van del Pol equation

• Van der Pol equation for electronic circuit
• Convert to two first-order ODEs

46
Stiff ODE van del Pol equation

• M-file for van del Pol equation

function yp vanderpol(t,y,mu) van der Pol
equation for electronic circuit d2y1/dt2 -
mu(1-y12) dy1/dt y1 0 initial
conditions, y1(0) dy1/dt 1 convert to two
first-order ODEs The solution becomes
progressively stiffer as mu gets large yp
y(2) mu(1-y(1)2)y(2)-y(1)
47
Nonstiff ODE van del Pol equation
gtgt t,yode45(_at_vanderpol,0 20,1 1,, 1) gtgt
Hplot(t,y(,1),t,y(,2),'m--') gtgt
legend('y1','y2',2)

• ? 1

ode45
48
Stiff ODE van del Pol equation
gtgt t,yode23s(_at_vanderpol,0 6000,1 1,,
1000) gtgt Hplot(t,y(,1)) gtgt set(H,'LineWidth',3
)

• ? 1000

ode23s
49
Bungee Jumper Coupled ODEs

• Vertical Dynamics of a jumper connected to a
  stationary platform with a bungee cord (Ex21.4)
• Air resistance depending on whether the cord is
  slack or stretched use sign(v) for drag force
• Spring constant k (N/m) and damping coefficient ?
  (Ns/m)
• Need to solve two simultaneous ODEs for x and v

k ? 0 if x ? L
50
Non-stiff ODEs Bungee Jumper

• Use ode45 for non-stiff ODEs to solve for the
  distance and velocity of a bungee jumper

function dydt bungee(t,y,L,cd,m,k,gamma) g
9.81 cord 0 if y(1) gt L determine if the
cord exerts a force cord k/m (y(1)-L)
gamma/m y(2) end dydt y(2) g - sign(y(2))
cd/my(2)2 - cord
51
gtgt t,y ode45(_at_bungee,0 50,0
0,,30,0.25,68.1,40,8) gtgt h1
plot(t,-y(,1),t,y(,2),'m--') gtgt h2 legend('x
(m)','v(m/s)') gtgt set(h1,'LineWidth',3)
set(h2,'FontSize',12)
Bungee jumper velocity
Bungee jumper distance
52
CVEN 302-501Homework No. 14

• Chapter 21
• Prob. 21.3 (30) (Hand Calculation) 21.5 (30)
• Prob. 21.8 (40) (using MATLAB program except for
  part a))
• Due 12/02/08 before 500 pm