Numerical Methods in Engineering (ENGR20005)
Lecture 19 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L19.1: Runge-Kutta Methods
2
Runge-Kutta Methods (Page 129 printed lecture notes)
The Runge-Kutta method is arguably, the most popular method used in the industry. It is fast, accurate and quite stable
In general, Runge-Kutta solution to
dx dt
can be written as
xn+1 = xn + Δtg(tn, xn, Δt)
g = a1k1 + a2k2 + a3k3 + a4k4 + . . . . . + aNkN
= f(t,x)
(19.1)
where
(19.2)
(19.3)
To solve
dx dt
= f(t,x) approximate solution as
+a3k3 +a4k4 +…..+aNkN
(19.2)
(19.3)
(19.4)
xn+1 = xn + Δtg(tn, xn, Δt)
g = a1k1 + a2k2
k1 =f(tn,xn)
k2 = f(tn + p1Δt, xn + q11k1Δt)
k3 = f (tn + p2Δt, xn + q21k1h + q22k2Δt)
k4 = f (tn + p3Δt, xn + q31k1h + q32k2h + q33k3Δt)
⋮⋮ ⋮⋮ ⋮⋮
kN = f(tn + pN−1Δt, xi + qN−1,1k1Δt + qN−1,2k2Δt + ⋯ + qN−1,N−1kN−1Δt)
The more terms you have (i.e. the more kis you retain), the more accurate the solution.
The 2nd order Runge-Kutta method formula can be obtained by setting N = 2 in Eq. (19.3)
xn+1 = xn + (a1k1 + a2k2) Δt k1 =f(tn,xn)
k2 = f(tn + p1Δt, xn + q11k1Δt)
In order to use the formula above, you need to find a1, a2, q11 and p1.
To solve
dx dt
= f(t,x) approximate solution as
xn+1 = xn + (a1k1 + a2k2) Δt
k1 =f(tn,xn)
k2 = f(tn + p1Δt, xn + q11k1Δt)
2nd order Runge-Kutta formula
Taylor series
Strategy to find a1, a2, q11 and p1.
Compare
What we would do is to first analyse the Runge-Kutta solution first and then compare with the Taylor series expansion to come up with values for a1, a2, q11 and p1.
x(tn+1)=x(tn)+Δtdx +Δt2 d2x +Δt3 d3x +…….. 1!dt tn 2! dt2 tn 3! dt3 tn
k1 = f(tn, xn)
k2 = f(tn + p1Δt, xn + q11k1Δt)
(19.5) (19.6)
xn+1 = xn + (a1k1 + a2k2) Δt (19.7) Put Eqs. (19.5) and (19.6) into (19.7) gives
xn+1 = xn + Δt (a1 f(tn, xn) + a2 f(tn + p1Δt, xn + q11k1Δt))
The third term on the RHS is a function of two variables. We can perform Taylor expansion
about the two variables, i.e.
f(t+r,x+s)=f(t,x)+r∂f +s∂f +…. ∂t ∂x
xn+1 = xn + Δt (a1 f(tn, xn) + a2 f(tn + p1Δt, xn + q11k1Δt)) f(t+r,x+s)=f(t,x)+r∂f +s∂f +….
∂t ∂x xn+1 = xn + Δt (a1 f(tn, xn) + a2 f(tn + p1Δt, xn + q11k1Δt))
=xn+Δt a1f(tn,xn)+a2 f(tn,xn)+∂fp1Δt+ ∂fq11k1Δt…..
( ( ∂t ∂x ))
=xn+Δt a1f(tn,xn)+a2 f(tn,xn)+∂fp1Δt+ ∂f f(tn,xn)q11Δt+……
( ( ∂t ∂x ))
= xn + (a1 + a2)Δtf(tn, xn) + a2 p1 ∂f (Δt)2 + a2q11 f(tn, xn) ∂f (Δt)2 + . . . . . ∂t ∂x
(19.8)
2nd order Runge-Kutta formula
Taylor series
Strategy to find a1, a2, q11 and p1
Compare
2nd order Runge-Kutta formula
Taylor series
Strategy to find Compare a1, a2, q11 and p1.
Recall in Lectures 17 and 18, we derived Euler method and Taylor method order 2 by truncating Taylor series at Δt2 & Δt3 respectively
+Δt3 d3x 3! dt3 tn
Euler
x(tn+1)=x(tn)+Δtdx +Δt2 d2x 1!dt tn 2! dt2 tn
x(t )=x(t)+Δt +Δt2 n+1 n 1! t 2!
+……..
(19.9)
(19.10)
dx dt
= f(t, x)
dx dt
d2x dt2
+Δt3 d3x 3! dt3 t
Taylor
Method order
2 nnn
+……..
Initial value problem can be expressed as
Put into Eq. (19.10) and neglecting higher order terms give
x(tn+1) = x(tn) + Δt f(tn, xn) + Δt2 df (tn, xn) 1! 2! dt
t
xn+1 = xn + (a1 + a2)Δtf(tn, xn)
+a2p1∂f(Δt)2 +a2q11f(tn,xn)∂f(Δt)2 +….. ∂t ∂x
x(tn+1) = x(tn) + Δt f(tn, xn) + Δt2 df (tn, xn)
1! 2 (2! dt )
x(tn+1)=x(tn)+Δtf(tn,xn)+Δt ∂f(tn,xn)+ ∂f(tn,xn)dx 1! 2! ∂t ∂x dt
x(tn+1) = x(tn) + Δt f(tn, xn) + Δt2 ( ∂f (tn, xn) + ∂f (tn, xn)f(tn, xn)) 1! 2!∂t ∂x
The equation above is called the Taylor Method order 2. It was shown in Lecture 18 that is it more accurate than Euler method. The main disadvantage is that we have to calculate the derivatives ∂f/∂t and ∂f/∂x in order to use it to compute solution to ODE. The 2nd order Runge-Kutta method allows us to get the accuracy of Eq. (19.11) without having to calculate the derivatives.
Strategy to find Compare a1, a2, q11 and p1
(19.11)
xn+1 =xn+(a1+a2)Δtf(tn,xn)+a2q11f(tn,xn)∂f(Δt)2+a2p1∂f(Δt)2+….. ∂x ∂t
2nd order Runge-Kutta formula
Taylor series
x(tn+1) = x(tn) + 1Δtf(tn, xn) + 1 ∂f (tn, xn)f(tn, xn)Δt2 + 1 ∂f (tn, xn)Δt2 2 ∂x 2 ∂t
2nd order Runge-Kutta formula
Taylor series
Strategy to find a1, a2, q11 and p1.
Compare Eq. (19.8)
with Eq. (19.11)
gives
a2q11 = 12
Compare
xn+1 =xn+(a1+a2)Δtf(tn,xn)+a2q11f(tn,xn)∂f(Δt)2+a2p1∂f(Δt)2+….. ∂x ∂t
x(tn+1) = x(tn) + 1Δtf(tn, xn) + ∂f (tn, xn)f(tn, xn)Δt2 + ∂f (tn, xn)Δt2 ∂x ∂t
1 2
1 2
a1 + a2 = 1 a 2 p 1 = 12
a1 + a2 = 1 a 2 p 1 = 12
a2q11 = 12
We have 4 unknowns, a1, a2, p1 & q11 but only 3 equations, hence there cannot be a unique solution.
One possible solution is to set a2 = 1/2 so that we have a 1 = 12
p1 = 1 q11 = 1
Thus the 2nd order Runge-Kutta method can be written as
k1 = f(tn, xn)
k2 =f(tn+p1Δt,xn+q11k1Δt)
xn+1 =xn+(a1k1+a2k2)Δt
k1 = f(tn, xn)
k2 =f(tn+Δt,xn+Δtk1)
xn+1 =xn+(12k1+ k2)Δt
1 2
Note that there is no unique combinations of a1, a2, p1 & q11. For the equation above, we chose a2 = 1/2. If you choose different value of a2, you will get different formula for 2nd order Runge- Kutta method.
Example L19.1:
Write a MATLAB program that uses 2nd order Runge-Kutta method, solve
=1−x For 0 ≤ t ≤ 8 with x(t = 0) = 0 and
a) Δt = 2 b) Δt = 1 c) Δt = 0.5 d) Δt = 0.1
dx dt
Lecture19M01.m
clear all;
close all;
Delta_t=1.0;
t=0:Delta_t:8.0;
0≤t≤8 %Preallocating memory
x=zeros(size(t));
%Initial condition
x(t = 0) = 0 for n=1:length(t)-1
k1=f(t(n),x(n));
end
x(1)=0.0;
k2=f(t(n)+Delta_t,x(n)+Delta_t*k1);
x(n+1)=x(n)+(Delta_t/2)*(k1+k2);
plot(t,x,’ko-‘,’linewidth’,2,’Markersize’,
10)
hold on
tplot=0:0.01:8.0;
truesolution=@(t)(1-exp(-t))
plot(tplot,truesolution(tplot),’b-‘,’linew
idth’,2)
xlabel(‘t’);
ylabel(‘x’);
legend(‘RK-2’,’True’);
axis([0 8 0 2])
function dxdt=f(t,x)
dxdt=1-x;
end
f(t, x) = 1 − x
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + Δtk1)
x n + 1 = x n + ( 12 k 1 + 12 k 2 ) Δ t
Δt = 1.0 Runge-Kutta 2nd order
Δt = 1.0 Explcit Euler
Δt = 1.0 Taylor 2nd order
Δt = 1.0 Runge-Kutta 2nd order
The Taylor 2nd order
solution and the
Runge-Kutta 2nd order (RK2) solution are very similar.
This is not surprising because we have chosen the
coefficients of RK2 to “mimic” the behaviour of Taylor 2nd order.
a1 + a2 = 1 a 2 p 1 = 12
a2q11 = 12
One possible solution is to set a2 = 1/2 so that we have
a 1 = 12 p1 = 1
k1 =f(tn,xn)
k2 = f(tn + p1Δt, xn + q11k1Δt)
q11 = 1
What happens if we change our coefficients? Let a2 = 1
xn+1 = xn + (a1k1 + a2k2) Δt
End of Example L19.1
To get better accuracy, one can use the 4th order Runge-Kutta (RK-4) method
To solve
= f(t, x)
Use 4th-order Runge-Kutta (page 133 of printed lecture notes):
dx dt
k1 =f(tn,xn)
k2 = f (tn + Δt , xn + Δt k1)
k3=f tn+Δt,xn+Δtk2 (22)
22 k4 = f (tn + Δt, xn + Δtk3)
xn+1 =xn+(16k1+13(k2+k3)+16k4)Δt
The global truncation error for (RK-4) can be shown to be O (Δt4).
4th order Runge-Kutta
k1 =f(tn,xn)
k2 = f(tn + Δt,xn + Δtk1)
k3=f tn+Δt,xn+Δtk2 (22)
k4=f tn+Δt,xn+Δtk3 (22)
xn+1 =xn+(16k1+13(k2+k3)+16k4)Δt
Example L19.2: Modify MATLAB program that uses the 2nd- order Runge-Kutta (RK2) method to instead use the 4th-order Runge-Kutta (RK4) method, to solve
=1−x for 0 ≤ t ≤ 8 with x(t = 0) = 1.
dx dt
2nd order Runge Kutta
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + Δtk1)
x n + 1 = x n + ( 12 k 1 + 12 k 2 ) Δ t
4th order Runge Kutta
k1 =f(tn,xn)
k2 = f (tn + Δt , xn + Δt k1)
k3=f tn+Δt,xn+Δtk2 (22)
22 k4 = f (tn + Δt, xn + Δtk3)
xn+1 =xn+(16k1+13(k2+k3)+16k4)Δt
clear all;
close all;
Delta_t=1.0;
t=0:Delta_t:8.0;
%Preallocating memory
x=zeros(size(t));
%Initial condition
x(1)=0.0;
for n=1:length(t)-1
k1=f(t(n),x(n));
Lecture19M02.m
k2=f(t(n)+Delta_t/2,x(n)+Delta_t*k1/2);
k3=f(t(n)+Delta_t/2,x(n)+Delta_t*k2/2);
k4=f(t(n)+Delta_t,x(n)+Delta_t*k3);
x(n+1)=x(n)+(Delta_t/
6)*(k1+2.0*(k2+k3)+k4);
end
clear all;
close all;
Delta_t=1.0;
t=0:Delta_t:8.0;
%Preallocating memory
x=zeros(size(t));
%Initial condition
x(1)=0.0;
for n=1:length(t)-1
k1=f(t(n),x(n));
Lecture19M01.m
k2=f(t(n)+Delta_t,x(n)+Delta_t*k1);
x(n+1)=x(n)+(Delta_t/2)*(k1+k2);
end
plot(t,x,’ko-‘,’linewidth’,2,’Markersiz
e’,10)
hold on
2nd order Runge Kutta
t=0:Delta_t:8.0;
%Preallocating memory
x=zeros(size(t));
4th order Runge Kutta
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + Δtk1)
1 ) n+1 n 21 22
closeall; (1 Delta_t=1.0;
clear all;
x =x+ k+ k Δt
k1 =f(tn,xn)
k2 = f (tn + Δt , xn + Δt k1)
k3=f tn+Δt,xn+Δtk2 (22)
22 k =f t +Δt,x +Δtk
clear all;
4 (n n 3)
closeall; (1 1
1 ) n+1 n 61 3 2 3 64
Delta_t=1.0; ( )
t=0:Dexlta_=t:x8.0+; k + k +k + k Δt
%Preallocating memory
x=zeros(size(t));
%Initial condition
x(1)=0.0;
for n=1:length(t)-1
end
plot(t,x,’ko-‘,’linewidth’,2,’Markersiz
e’,10)
hold on
tplot=0:0.01:8.0;
truesolution=@(t)(1-exp(-t))
plot(tplot,truesolution(tplot),’b-‘,’li
newidth’,2)
xlabel(‘t’);
ylabel(‘x’);
Lecture19M01.m
Lecture19M02.m
k1=f(t(n),x(n));
k2=f(t(n)+Delta_t,x(n)+Delta_t*k1);
x(n+1)=x(n)+(Delta_t/2)*(k1+k2);
k1=f(t(n),x(n));
k2=f(t(n)+Delta_t/2,x(n)+Delta_t*k1/2);
k3=f(t(n)+Delta_t/2,x(n)+Delta_t*k2/2);
k4=f(t(n)+Delta_t,x(n)+Delta_t*k3);
x(n+1)=x(n)+(Delta_t/
6)*(k1+2.0*(k2+k3)+k4);
legend(‘RK-2’,’True’);
%Initial condition
x(1)=0.0;
for n=1:length(t)-1
end
plot(t,x,’ko-‘,’linewidth’,2,’Markersize’,10)
hold on
tplot=0:0.01:8.0;
truesolution=@(t)(1-exp(-t))
plot(tplot,truesolution(tplot),’b-‘,’linewidt
h’,2)
Lecture19M02.m
k1 =f(tn,xn)
k2 = f (tn + Δt , xn + Δt k1)
k3=f tn+Δt,xn+Δtk2 (22)
22 k4 = f (tn + Δt, xn + Δtk3)
xn+1 =xn+(16k1+13(k2+k3)+16k4)Δt
k1=f(t(n),x(n));
k2=f(t(n)+Delta_t/2,x(n)+Delta_t*k1/2);
k3=f(t(n)+Delta_t/2,x(n)+Delta_t*k2/2);
k4=f(t(n)+Delta_t,x(n)+Delta_t*k3);
x(n+1)=x(n)+(Delta_t/6)*(k1+2.0*(k2+k3)+k4);
Δt = 1.0 Runge-Kutta 4th order
plot(t,x,’ko-‘,’linewidth’,2,’Markersize’,10)
hold on
tplot=0:0.01:8.0;
truesolution=@(t)(1-exp(-t))
plot(tplot,truesolution(tplot),’b-‘,’linewidth’,2)
xlabel(‘t’);
ylabel(‘x’);
legend(‘RK-4′,’True’);
axis([0 8 0 2])
clear all;
close all;
Delta_t=1.0;
t=0:Delta_t:8.0;
%Preallocating memory
x=zeros(size(t));
%Initial condition
x(1)=0.0;
for n=1:length(t)-1
end
function dxdt=f(t,x)
dxdt=1-x;
end
f(t, x) = 1 − x
Δt = 1.0 Explicit Euler
Solve dx=1−x dt
for 0 ≤ t ≤ 8 with x(0) = 1. Analytical
solution:
x(t) = 1 − e−t
Δt = 1.0 Runge-Kutta 2nd order
Δt = 1.0 Runge-Kutta 4th order
End of Example L19.2
L19.2: Stability Analysis: Runge-Kutta Methods
24
Stability analysis of 2nd-order Runge-Kutta (page 137 of printed lecture notes): Apply 2nd-order Runge-Kutta
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + Δtk1)
x n + 1 = x n + ( 12 k 1 + 12 k 2 ) Δ t
to model problem
dx dt
= λx = f(t,x)
where λ is a constant and is a complex number
λRe + iλIm
i = −1 is the imaginary unit.
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + Δtk1)
xn+1 = xn + ( 12 k1 + 12 k2) Δt
(19.11)
= λx = f(t,x)
Substitute f(t, x) = λx into expressions for k1 and k2
dx dt
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + Δtk1)
k1 = λxn
k2 = λ(xn + Δtk1)
k2 = λ(xn + Δt(λxn))
Note that this expression for σ is identical to σ for the
Taylor Method Order 2
Substitute k1 and k2 into Eq. (19.11) gives x = x + Δtλx + 1 Δt2λ2x
n+1 n n
x = (1 + Δtλ + Δt2λ2/2)x
n+1 n n xn+1 = σxn = σ x0
where
σ = 1 + Δtλ + Δt2λ2/2
2 n
xn+1 = σxn Thus, running a program using the RK-2 method will give you
x0 = Given
x1 = σx0 2
x2 = σx1 = σ x0 x3 = σx2 = σ3×0
⋮=⋮=⋮ xn+1 = σxn = x0
σn
For solution to be stable
|σ| ≤ 1
We can get MATLAB to plot contour of
|σ| = 1+λΔt+ λ Δt 2
2 2
= 1
to identify stability region for 2nd order Runge-Kutta method. For 4th order Runge-Kutta method
k1 =f(tn,xn)
k2 = f (tn + Δt , xn + Δt k1)
k3=f tn+Δt,xn+Δtk2 (22)
22 k4 = f (tn + Δt, xn + Δtk3)
xn+1 =xn+(16k1+13(k2+k3)+16k4)Δt It can be shown that the stability region is given by areas where
|σ|= 1+λΔt+λ2Δt2 +λ3Δt3 +λ4Δt4 ≤1 2 6 24
and the stability boundary can be illustrated by getting MATLAB to plot contour of
1+λΔt+λ2Δt2 +λ3Δt3 +λ4Δt4 =1 2 6 24
-2.79
λImΔt -2.0
1+λΔt+λ2Δt2 +λ3Δt3 +λ4Δt4 =1 2 6 24
1+λΔt+λ2Δt2 =1 2
2.83
4th order Runge Kutta
2nd order Runge Kutta Explicit Euler
-2.83
λReΔt
1+λΔt+λ2Δt2 =1 2
•Stability region for RK4 is bigger than RK2
•In order for your code to be stable,
you need to pick t such that it falls within the stability region of your numerical scheme.
• For example, if you want to solve
dx dt
= − 10x
using explicit Euler, the maximum value of Δt you
can use is
tmax,Euler = 2 = 1 = 0.2 10 5
• If you were to use RK4, then the maximum value of Δt you can use is
tmax,RK4 = 2.79 = 0.279 10
1+λΔt+ λ2Δt2 + λ3Δt3 + λ4Δt4 =1 2 6 24
L19.3:
The Implicit Euler Method
31
(see page 125 of “Book”) Our analysis thus far shows that the explicit Euler’s method is not very accurate or stable. To come up with a method that is more stable, consider the Taylor series expanded about t = tn+1
x(tn)=x(tn+1)−Δtdx(tn+1)+Δt2 d2x(tn+1)−Δt3 d3x(tn+1)+….. 1! dt 2! dt2 3! dt3
Ignoring higher order terms give
x(tn) = x(tn+1) − Δt f(tn+1, xn+1) 1!
Hence,
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
The above equation is called the implicit Euler method.
x(tn+1) = x(tn) + Δtf(tn+1, xn+1) Eq. (19.12) is the formula for the implicit Euler method.
Compare now with the formula for the explicit Euler method
x(tn+1) = x(tn) + Δtf(tn, xn)
The only difference is the time level the right hand side, dx/dt = f(t, x) is evaluated at. For the explicit Euler method, the function
For the implicit Euler method, the the function
This will give rise to complications, especially if f(t, x) is a nonlinear function (see later).
(19.12)
f (t, x) is evaluated at time level tn.
f (t, x) is evaluated at time level tn+1.
x(t)
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Illustration of the implicit Euler method
x(tn+1) = x(tn) + Δt dx (tn+1, xn+1) dt
dx (tn+1, xn+1) dt
True value x(tn+1)
xn
Δt
tn
tn+1
x(t)
x(tn+1) = x(tn) + Δt dx (tn+1, xn+1)
Smaller Δt will lead to smaller error
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
xn
tn
dt
True value x(tn+1) Error
Predicted xn+1 Δt
tn+1
x(t)
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Smaller Δt will lead to smaller error x(tn+1) = x(tn) + Δt dx (tn+1, xn+1)
dt
True value of x(tn+1)
dx (tn+1, xn+1) dt
xn
Δt
tn+1
tn
x(t)
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Smaller Δt will lead to smaller error x(tn+1) = x(tn) + Δt dx (tn+1, xn+1)
dt
True value of x(tn+1)
xn
Error Predicted xn+1
t
tn+1
tn
x(t)
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Smaller Δt will lead to smaller error
xn
dx (tn+1, xn+1) dt
True value of x(tn+1) t
tn+1
x(tn+1) = x(tn) + Δt dx (tn+1, xn+1) dt
tn
x(t)
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Smaller Δt will lead to smaller error x(tn+1) = x(tn) + Δt dx (tn+1, xn+1)
dt True value of x(tn+1)
xn
Predicted xn+1 t
Error
tn
tn+1
Example L19.3:
Using Euler’s method, solve
dx dt
For 0 ≤ t ≤ 8 with x(t = 0) = 0 and a) Δt = 2
b) Δt = 1 c) Δt = 0.5 d) Δt = 0.1
=1−x
Start with the implicit Euler’s formula
xn+1 = xn + Δtf(tn+1, xn+1) Replace f (t, x) with the function given in the question
xn+1 = xn + Δt(1 − xn+1)
Rearranging gives
xn+1 =
xn + Δt 1+Δt
Lecture19M03.m
xn+1 =
xn + Δt 1+Δt
clear all;
close all;
Delta_t=2.0;
t=0:Delta_t:8.0;
%Preallocating memory
x=zeros(size(t));
%Initial condition
x(1)=0.0;
for n=1:length(t)-1
x(n+1)=(x(n)+Delta_t)/(1+Delta_t);
end
plot(t,x,’ko-‘,’linewidth’,2,’Markersize’,10
)
hold on
tplot=0:0.01:8.0;
truesolution=@(t)(1-exp(-t))
plot(tplot,truesolution(tplot),’b-‘,’linewid
th’,2)
xlabel(‘t’);
ylabel(‘x’);
legend(‘Implicit Euler’,’True’);
axis([0 8 0 2])
Output
Δt = 2.0
Δt = 2.0
Δt = 1.0
Δt = 0.5
Δt = 0.1
Compare implicit and explicit Euler’s method for Δt = 2.0
Δt = 2.0
Δt = 2.0
Implicit Euler’s method is much more stable than the explicit Euler’s method
Stability analysis of the implicit Euler Method
We will now analyse why the Implicit Euler method is so stable. Again, consider the model problem
= λx Applying the implicit Euler method gives
dx dt
xn+1 = xn + Δtλxn+1 (1 − λΔt)xn+1 = xn 1
xn+1 = 1−λΔtxn xn+1 = σxn
1 1−λΔt
σ=
Observe that
|σ| = =
Hence, the stability region is given by
1 1−λΔt
1
(1 − λReΔt)2 + (λImΔt)2
(1 − λReΔt)2 + (λImΔt)2 ≥ 1 This is any area outside a circle of radius 1 centred at (1,0)
(1 − λReΔt)2 + (λImΔt)2 ≥ 1
Im t
1
2
Note that the implicit Euler method is will be stable as even when λRe is positive.
This does not mean the implicit Euler method is accurate, only that the method is very stable. In general, the implicit Euler method is much more stable than the explicit Euler method.
Re t
Consider the initial value problem dx = x
dt t
with x(t = 0) = 1 and t ∈ [0,8]. The analytical (true) solution is x(t) = e which grows with time. However, if you use implicit Euler method to solve this problem, you could have a solution that decays with time (which is the completely wrong behaviour!). See figure below showing the true solution together with computed solution using the explicit and implicit Euler methods using Δt = 2.
-2
-1
x(t) = et
Im t
Explicit Euler
Implicit Euler
1
λΔt = 1 × 2 =2
Re t