Numerical Methods for IVP ODEs
CS/SE 4X03
Ned Nedialkov McMaster University March 19, 2021
Outline
The problem ODE examples ODEs
Euler’s method Backward Euler Stability
The problem Examples ODEs Euler’s method The problem
Backward Euler
Stability
• Given
compute y(t) on [a, b]
• y′ ≡ y′(t) ≡ dy dt
y′ =f(t,y),
y(a)=c
• This is an Initial Value Problem (IVP) in Ordinary Differential Equations (ODEs)
• We approximate y(t) at points ti in [a, b] using a numerical method
Copyright © 2021 N. Nedialkov. All rights reserved. 3/23
The problem Examples ODEs Euler’s method Backward Euler Stability ODE examples
y′ = −y + t • Solutionisy(t)=t−1+αe−t
• Giveny(0)=c,e.g. c=5,
y(0) = −1 + α = c = 5,
y(t) = t − 1 + 6e−t is the solution with this initial condition
α = 6
Copyright © 2021 N. Nedialkov. All rights reserved.
4/23
The problem Examples ODEs Euler’s method Backward Euler Stability
θ′′ =−gsinθ, θ′′ = d2θ(t) dt2
• Motion of a pendulum
◦ ball of mass 1 attached to the end of a rigid, massless rod of
length r = 1
◦ g ≈ 9.81 is gravity ◦ t is time
• This is a second-order ODE
• Towriteasafirst-orderODE,sety1 =θ,y2 =θ′ =y1′:
y 1′ = y 2
y2′ =−gsin(y1)
• Needed initial conditions are y1(0) and y2(0)
Copyright © 2021 N. Nedialkov. All rights reserved. 5/23
The problem Examples ODEs Euler’s method Backward Euler ODEs
Stability
System of n first-order equations in n variables
y′ =f(t,y), f :R×Rn →Rn
Nonlinear: if f is nonlinear in y, linear otherwise
Copyright © 2021 N. Nedialkov. All rights reserved.
6/23
The problem Examples ODEs Euler’s method Backward Euler Stability Autonomous ODE
y′ = f(y), y(a) = c
is an autonomous ODE, does not depend on time explicitly
y′ =f(t,y), y(a)=c
To convert a non-autonomous ODE to an autonomous set x = t
is non-autonomous and then
x′ = 1
y′ =f(x,y), x(a)=a, y(a)=c
Setz=(z1,z2)T =(x,y)T. Thenz′ =f(z): z 1′ = 1
z 2′ = f ( z )
Copyright © 2021 N. Nedialkov. All rights reserved. 7/23
The problem Examples ODEs Euler’s method Backward Euler High-order ODEs
Stability
y(n) = f t,y,y′,…,y(n−1) can be converted to first-order by setting
y1 = y
y 2 = y ′ = y 1′ y3 = y′′ = y2′
.
yn = y(n−1) = yn′ −1 Copyright © 2021 N. Nedialkov. All rights reserved.
8/23
The problem
Examples
ODEs Euler’s method Backward Euler
Stability
Then
y 1′ = y 2 y 2′ = y 3
.
yn′ = f(t,y1,y2,…,yn) To solve, we need initial values for
y1(a), y2(a), . . . , yn(a)
Copyright © 2021 N. Nedialkov. All rights reserved.
9/23
The problem Examples ODEs Euler’s method Backward Euler Stability Euler’s method
• Leth=(b−a)/N,N>1isaninteger • h is stepsize
• Lett0 =a,ti =a+ih,i=0,1,…,N • From y′(ti) = f(ti,y(ti)), we write
′ h2′′ y(ti+1) = y(ti) + hy (ti) + 2 y (ξi),
ξi between ti andti+1 =y(ti)+hf(ti,y(ti))+ 2y (ξi)
h2 ′′ ≈ y(ti) + hf(ti, y(ti))
Copyright © 2021 N. Nedialkov. All rights reserved. 10/23
The problem Examples ODEs Euler’s method Backward Euler Stability • Euler’s method:
y0 = c
yi+1 =yi+hf(ti,yi), i=0,1,…,N−1
• Example: Euler’s method on y′ = −y + t, y(0) = y0 = 5, with h = 0.1:
yi+1 =yi +hf(ti,yi)=yi +h(−yi +ti)
y1 =y0 +h(−y0 +t0)=5+0.1(−5+0)=4.5
y2 =y1 +h(−y1 +t1)=4.5+0.1(−4.5+0.1)=4.06
y3 =y2 +h(−y2 +t2)=4.06+0.1(−4.06+0.2)=3.674
• Exactsolutionisy(t)=t−1+6e−t
• The corresponding exact values are y(0.1) ≈ 4.5290,
y(0.2) ≈ 4.1124, y(0.3) ≈ 3.7449
Copyright © 2021 N. Nedialkov. All rights reserved. 11/23
The problem Examples ODEs Euler’s method Backward Euler Stability Example: Forward Euler on y′ = −y + t
5 4.5 4 3.5 3 2.5 2 1.5
y’ = -y+t’, y(0) = 5, h = 0.1
y(t) = t-1+6*exp(-t)
initial condition forward Euler
0 0.5 1 1.5 2 2.5 3
Copyright © 2021 N. Nedialkov. All rights reserved.
12/23
The problem Examples ODEs Euler’s method Backward Euler Stability Backward Euler
• We can write
y(ti) = y(ti+1) − hy (ti+1) + 2 y (ηi)
′ h2 ′′ ≈ y(ti+1) − hf(ti+1, y(ti+1))
y(ti+1) ≈ y(ti) + hf(ti+1, y(ti+1)) • Backward Euler
yi+1 =yi +hf(ti+1,yi+1)
• This is an implicit method; forward Euler is explicit
Copyright © 2021 N. Nedialkov. All rights reserved.
13/23
The problem Examples ODEs Euler’s method Backward Euler Stability
• Example: Backward Euler method on y′ = −y + t, y(0)=y0 =5,withh=0.1:
yi+1 =yi +hf(ti+1,yi+1) = yi + h(−yi+1 + ti+1)
• We need to solve for yi+1:
yi+1 = yi − hyi+1 + hti+1
yi+1 + hyi+1 = yi + hti+1 yi+1 = yi + hti+1
1+h
Copyright © 2021 N. Nedialkov. All rights reserved.
14/23
The problem Examples ODEs Euler’s method Backward Euler Stability
• We compute
y1 = y0 +ht1 = 5+0.1·0.1 ≈4.5545 1+h 1+0.1
y2 = y1 +ht2 ≈ 4.5545+0.1·0.2 ≈4.1586 1+h 1+0.1
y3 = y2 +ht3 ≈ 4.1586+0.1·0.3 ≈3.7987 1+h 1+0.1
• The corresponding exact values are y(0.1) ≈ 4.5290, y(0.2) ≈ 4.1124, y(0.3) ≈ 3.7449
• Hereitwaseasytosolveforyi+1: f(t,y)=−y+tislinearin y
• In general, it is non-linear: apply Newton’s method
Copyright © 2021 N. Nedialkov. All rights reserved. 15/23
The problem Examples ODEs Euler’s method Backward Euler Stability Example: FEandBEony′ =−y+t
5 4.5 4 3.5 3 2.5 2 1.5
y’ = -y+t’, y(0) = 5, h = 0.1
y(t) = t-1+6*exp(-t)
initial condition forward Euler backward Euler
0 0.5 1 1.5 2 2.5 3
Copyright © 2021 N. Nedialkov. All rights reserved.
16/23
The problem Examples ODEs Euler’s method Backward Euler Stability Stability
Forward Euler
• Consider y′ = λy, y(0) = y0
• The exact solution is y(t) = eλty0
• Forward Euler with constant stepsize h is
yi+1 =yi +hf(ti,yi)=yi +hλyi = (1 + hλ)yi
= (1 + hλ)2yi−1 .
= (1 + hλ)i+1y0
• If λ < 0, y(t) is decaying. Since |y(ti+1)| < |y(ti)|, we want
|yi+1| ≤ |yi|
Copyright © 2021 N. Nedialkov. All rights reserved. 17/23
The problem Examples ODEs Euler’s method Backward Euler Stability Stability cont.
• For the method to be numerically stable, we require |yi+1| = |1 + hλ| · |yi| ≤ |yi|
• Thatis|1+hλ|≤1,or
−1 ≤ 1 + hλ ≤ 1
−2 ≤ hλ ≤ 0 h≤2
|λ|
• If |λ| is large, we can have a severe restriction on the stepsize
Ife.g. y′ =−106y,h≤2·10−6
Copyright © 2021 N. Nedialkov. All rights reserved. 18/23
The problem Examples ODEs Euler’s method Backward Euler
Stability cont.
Example 1.
• Consider y′ = −10y, y(0) = y0
• Euler’s method is
yi+1 =yi +hλyi =(1−10h)yi
• For stability h ≤ 0.2 • If e.g. h = 0.21 then
y1 =(1−10·0.21)y0 =−1.1y0 y2 = −1.1y1 = 1.21y0
y3 = −1.1y2 = −1.331y0
.
yi = (−1.1)iy0
Copyright © 2021 N. Nedialkov. All rights reserved.
Stability
19/23
The problem Examples ODEs Euler’s method
Stability
Example 1. cont. 11
Backward Euler
Stability
= -10, h = 0.05
= -10, h = 0.15
t ey
0
Forward Euler
0.8
0.6
0.4
0.2
0
0.5
0
-0.5
0 0.5 1 1.5 2 2.5 3
= -10, h = 0.21
t ey
0
Forward Euler
0 0.5 1 1.5 2 2.5 3
= -10, h = 0.2
14
2
0
-2
-4
t ey
0
Forward Euler
t ey
0
Forward Euler
0.5
0
-0.5
-1 -6 0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3
Copyright © 2021 N. Nedialkov. All rights reserved.
20/23
The problem Examples ODEs Euler’s method Backward Euler Stability Stability
Backward Euler
• Consider the backward Euler on y′ = λy, where λ < 0 yi+1 = yi + hλyi+1
yi+1= 1 yi 1−hλ
|yi+1| = 1 |yi| |1−hλ|
≤ |yi| for any h > 0
Copyright © 2021 N. Nedialkov. All rights reserved.
21/23
The problem Examples ODEs Euler’s method Backward Euler Stability Stability
Example 2.
• y′ = −10y
• Backward Euler is
yi+1= 1 yi 1+10h
• Stableforanyh>0
• Backward Euler is absolutely (for any h > 0) stable
Copyright © 2021 N. Nedialkov. All rights reserved.
22/23
The problem Examples ODEs Euler’s method
Stability
Example 2. cont. 11
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
00 0 0.5 1 1.5 2 2.5 3
= -10, h = 0.2
11
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
00 0 0.5 1 1.5 2 2.5 3
Backward Euler
Stability
= -10, h = 0.05
= -10, h = 0.15
t ey
0
Backward Euler
t ey
0
Backward Euler
0 0.5 1 1.5 2 2.5 3
= -10, h = 0.21
t ey
0
Backward Euler
t ey
0
Backward Euler
Copyright © 2021 N. Nedialkov. All rights reserved.
23/23
0 0.5 1 1.5 2 2.5 3