Initial value problems (IVPs)

• IVP: Given a differential equation and initial values

Numerical Methods for Initial Value Problems
dy
= f (t, y), and y(t0 ) = f0 . (1)
dt

• y: Dependent variable, t: independent variable

CoE Systems Science
• If f = f (y) 7→ Autonomous ODE
August 29, 2013
• Not all IVPs can be solved analytically 7→ Numerical methods

• Here we will focus on: First order linear IVPs. IVPs involving higher
order ODEs can be rendered as a system of first order ODEs.

Hiremath, K. R. (IITJ) 1 / 19 August 29, 2013

Euler method (Euler explicit or Euler forward) • If h is sufficiently small, then we can approximate

dy y(t1 ) = y(t0 + h) = y(t0 ) + hy 0 (t0 )

• Given = f (t, y), and y(t0 ) = y0 , and we want to find the solution
dt = y(t0 ) + hf (t0 , y(t0 )), (using given ODE)
for t0 ≤ t ≤ tend .
• Similarly
• Discretise the interval [t0 , tend ] with a constant stepsize h. y(t2 ) = y(t1 + h) = y(t1 ) + hf (t1 , y(t1 ))
y(t3 ) = y(t2 + h) = y(t2 ) + hf (t2 , y(t2 ))
tn = t0 + nh, yn = y(tn ). ..
.
• Knowing y at t0 , how to compute y at t1 = t0 + h?
y(tn+1 ) = y(tn ) + hf (tn , y(tn ))
Taylor series expansion:
This method is called forward Euler method, explicit Euler method, or
0
y(t1 ) = y(t0 + h) = y(t0 ) + hy (t0 ) + Et (h), the standard Euler method.

where • Discretization 7→ Errors

• Is the discretization right (consistent)?
h2 00 As h → 0, discrete problem → continuous problem.
Et (h) = y (ξ), t0 < ξ < (t0 + h) = O(h2 )
2
• Is the scheme stable?
is the local truncation error. Errors don’t grow unbounded.
Hiremath, K. R. (IITJ) 2 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 3 / 19 August 29, 2013
Error analysis of the Euler method: Consistency Error analysis of the Euler method: Stability

• Consistency: A numerical scheme is said to be consistent with the • Stability: Errors propagate 7→ Investigation of unbounded growth of
ODE being discretized, if the local truncation error tends to zero as the errors
mesh size tends to zero.
• Error analysis for general case 7→ can be difficult!
1 From Taylor series expansion 7→ Local truncation error Et (h) = O(h2 ) • To get insight: Simple case
2 Finite precision calculation 7→ Round off errors dy
= −ay, y(0) = y0 , a > 0(constant)
dt
3 Total error/ global error
• Analytic solution: y = y0 e−at (· · · asymptotically approached to 0.)
= Local truncation error propagation + Round off errors = O(h)
• Forward Euler scheme: yn+1 = yn − hayn = (1 − ah)yn
If |1 − ah| > 1 7→ The numerical solution will grow unbounded.
⇒ The stability of the method depends on the stepsize h. ⇐
The forward Euler method is conditionally stable.

• Systematic ways of analysing the stability: Von Neumann method

based on Fourier decomposition of error

Hiremath, K. R. (IITJ) 4 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 5 / 19 August 29, 2013

Example Improvements in the forward Euler method

Solve the IVP • Forward Euler method:

y 0 (t) = −4y, y(0) = 4
y(tn+1 ) = y(tn ) + hf (tn , y(tn ))
over an interval [0, 2] with the Euler forward method with h = 0.5, 0.2, 0.1,
and 0.05. A new value of y(t) = y(tn+1 ) is predicted using
1 the old value of y = y(tn ), and

2 the slope at the old value of t = tn (beginning of the subinterval)

to extrapolate linearly over the stepsize h.
• General scheme:
yn+1 = yn + hφ,
φ: increment function 7→ How to defined it?.

Hiremath, K. R. (IITJ) 6 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 7 / 19 August 29, 2013
Improvements in the forward Euler method One more example. . .
• Why not approximate the constant slope on [tn , tn+1 ] by • Consider
1 the slope at the end-point of the interval? (Backward Euler method) y 0 (t) = 1 + t + y 2 , y(0) = 0
| {z }
yn+1 = yn + hf (tn+1 , yn+1 ) (2) f (t,y)

2 the slope at the mid-point of the interval? The forward Euler scheme:
2
yn+1 = yn + hf (tn+1/2 , yn+1/2 ) (3) yn+1 = yn + h(1 + tn+1 + yn+1 )

3 the average of the slopes at tn and tn+1 ? yn+1 on both the sides 7→ Solve nonlinear algebraic equation.
 
f (tn , yn ) + f (tn+1 , yn+1 ) • This is a typical case when f (t, y) is nonlinear in y. 7→ Implicit scheme:
yn+1 = yn + h (4)
2 Implicit Euler method 7→ Computationally expensive
• Let’s start with the first choice.: Solve the IVP • :-(( Difficulty with the mid-point: yn+1 = yn + hf (tn+1/2 , yn+1/2 )
| {z }
0
y (t) = −4y, y(0) = 4 unknown
over an interval [0, 2] with the backward Euler method with h = 0.5, How to overcome the above difficulties?
and 0.2. Compare these results with that of the forward Euler method. Predictor-Corrector methods
Hiremath, K. R. (IITJ) 8 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 9 / 19 August 29, 2013

The mid-point method Heun’s method (Improved/modified Euler method)

• Motivation: Use the slope information (= f (tn+1/2 , yn+1/2 )) at the • Motivation: Use the slope information at the both end points of the
mid-point tn+1/2 of the interval [tn , tn+1 ]. interval [tn , tn+1 ].
• Difficulty: We don’t know yn+1/2 . • Difficulty: Possibility of getting an implicit scheme!

• Way out: • Way out: Repeat the trick used in the mid-point method :-B
1 Using the standard Euler method, predict the value of yn+1/2 :
• Heun’s method: (LTE=O(h3 ), GE= O(h2 ))
h
yn+1/2 = yn + f (tn , yn ) 0
2 yn+1 = yn + hf (tn , yn ) (Predictor)
 0 )
f (tn , yn ) + f (tn+1 , yn+1
2 Using the above predicted value of yn+1/2 , correct the value of yn+1 :
yn+1 = yn + h (Corrector)
yn+1 = yn + hf (tn+1/2 , yn+1/2 )
2
• Iterated Heun’s method: For m = 1, 2, · · · (+ Termination criterion)
• The mid-point method: (LTE=O(h3 ), GE= O(h2 ))
0
yn+1 = yn + hf (tn , yn ) (Predictor)
h
yn+1/2 = yn + f (tn , yn ) (Predictor)
" #
m−1
2 m f (tn , yn ) + f (tn+1 , yn+1 )
yn+1 = yn + h (Corrector)
yn+1 = yn + hf (, tn+1/2 , yn+1/2 ) (Corrector) 2
Hiremath, K. R. (IITJ) 10 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 11 / 19 August 29, 2013
Example Recap

Solve the IVP • The mid-point method:

0 h
y (t) = −4y, y(0) = 4 yn+1/2 = yn + f (tn , yn )
2
over an interval [0, 2] with the mid-point method and Heun’s method with yn+1 = yn + hf (tn+1/2 , yn+1/2 )
h = 0.5, 0.2, 0.1, and 0.05.. Compare these results with that of the forward
Euler method. • Heun’s method:
0
yn+1 = yn + hf (tn , yn )
f (tn , yn ) + f (tn+1 , ynm−1 )
 
m
yn+1 = yn + h , m = 1, 2, · · ·
2
• General idea:
1 Evaluate the slope function at various points within the time step.

2 Use the “weighted average” of all these slope values to extrapolate the y
value at the next time step.

• Generalize this concept for higher order accurate methods 7→ The

Runge-Kutta (RK) methods
Hiremath, K. R. (IITJ) 12 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 13 / 19 August 29, 2013

Back to basics. . . • Inspect the equation

h
yn+1 = yn + f (tn , yn )
Suppose we want to construct a second order accurate method. 2
• Taylor series: h
h2 + [f (tn , yn ) + hft (tn , yn ) + hfy (tn , yn )f (tn , yn )] + O(h3()8)
yn+1 = y(tn+1 ) = y(tn +h) = y(tn )+hy 0 (tn )+ y 00 (tn )+O(h3 ) (5) 2
2 • Multivariate Taylor series expansion:
• From the given IVP:
f (tn + h, yn + k) = f (tn , yn ) + hft (tn , yi ) + kfy (tn , yn ) + · · ·
y 0 (tn ) = f (tn , yn ) (6)
 
00 d d dy
y (tn ) = f (t, y) + f (t, y) f (tn +h, yn +hf (tn , yn )) = f (tn , yn ) + hft (tn , yn ) + hf (tn , yn )fy (tn , yn )+O(h2 )
dt dy dt (tn ,yn )
• Thus Eq. (8) becomes
= ft (tn , yn ) + fy (tn , yn )y 0 (tn ) (7)
h h
• Using Eq.(7) in Eq.(5) gives yn+1 = yn + f (tn , yn ) + f (tn + h, yn + hf (tn , yn )) + O(h3 )
2 2
h2 
ft (tn , yn ) + fy (tn , yn )y 0 (tn ) + O(h3 ) • Rewriting 7→ Classical second order Runge-Kutta (RK2) method

yn+1 = yn + hf (tn , yn ) +
2  
h 1 1
yn+1 = yn + f (tn , yn ) yn+1 = yn + h k1 + k2
2 2 2
h k1 = f (tn , yn )
+ [f (tn , yn ) + hft (tn , yn ) + hfy (tn , yn )f (tn , yn )] + O(h3 )
2 k2 = f (tn + h, yn + hk1 )
Hiremath, K. R. (IITJ) 14 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 15 / 19 August 29, 2013
The classical second order Runge-Kutta method Generalized second order Runge-Kutta methods

 
1 1 yn+1 = yn + h (b1 k1 + b2 k2 )
yn+1 = yn + h k1 + k2
2 2
k1 = f (tn , yn )
k1 = f (tn , yn )
k2 = f (tn + β2 h, yn + α21 hk1 , )
k2 = f (tn + h, yn + hk1 )

1
1 b1 = b2 = 2 and α21 = β2 = 1 7→ Heun’s method
1 The k1 and k2 are known as the stages of RK method.
1
2 b1 = 0, b2 = 1 and α21 = β2 = 2 7→ The mid-point method
2 They correspond to different estimations for the slope of the solution.
• k1 : slope estimation at the start of the interval tn using yn Conclusion: several versions of RK2 method
• k2 : slope estimation at the end of the interval tn+1 using yn + hk1 Four unknowns and three equations (Ref: Chapra’s book)
This is also true for a general n’th order RK method.
3 This scheme can be interpreted as a predictor-corrector method. 7→
Heun’s method (with m = 1)
Hiremath, K. R. (IITJ) 16 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 17 / 19 August 29, 2013

Generalized m’th order explicit Runge-Kutta methods Classical 4th order Runge-Kutta method

m
X k1 = f (tn , yn )
yn+1 = yn + h bj kj
j=1
1 1
k2 = f (tn + h, yn + k1 h)
k1 = f (tn , yn ) 2 2
1 1
k2 = f (tn + β2 h, yn + α21 hk1 ) k3 = f (tn + h, yn + k2 h)
2 2
.. k4 = f (tn + h, yn + k3 h)
.
m−1 1
X yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )h
km = f (tn + βn h, yn + h αn,j kj ) 6
j=1

Most popular: Classical RK4 method • When f = f (t): Classical RK4 is similar to Simpson’s 1/3 rule

• Local truncation error: O(h5 )

• Global error: O(h4 )

Hiremath, K. R. (IITJ) 18 / 19 August 29, 2013 Hiremath, K. R. (IITJ) 19 / 19 August 29, 2013