Numerical Methods in Engineering (ENGR20005)
Lecture 18 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L18.1:
Higher Order Taylor Method
2
Higher Order Taylor Method
Recall from the last lecture, that we can compute the solution to
= f(t,x) using the explicit Euler formula xn+1 = xn + Δtf(tn, xn).
dx dt
In general, Euler’s method is not very accurate because it is derived by truncating Taylor series at Δt2
Taylor series
x(tn+1) = x(tn) + Δt dx (tn) 1! dt
Euler’s method
Ignore all these terms highlighted in blue
+Δt2 d2x(tn)+ Δt3 d3x(tn)+….. 2! dt2 3! dt3
But for initial value problem, you are usually given
x(tn+1) = x(tn) + Δtf(tn, xn)
dx dt
= f(t,x)
Higher Order Taylor Method
To get a better approximation of the solution let’s consider retaining more terms of the Taylor series.
Taylor series
x(tn+1) = x(tn) + Δt dx (tn) 1! dt
Taylor method order 2
x(tn+1)=x(tn)+Δt 1!
Retain this term
Ignore these terms
For initial value problem, you are usually given
dx(tn) dt
+Δt2 d2x(tn) 2! dt2
x(tn+1) = x(tn) + Δtf(tn, xn) + Δt2 2!
dx = f(t,x)
dt dt
d2x dt2
t=tn
Need to express this in terms of f (t, x).
+Δt2 d2x(tn)+Δt3 d3x(tn)+….. 2! dt2 3! dt3
x(tn+1) = x(tn) + Δtf(tn, xn) + Δt2 d2x 2! dt2
t=tn x(tn+1)=x(tn)+Δtf(tn,xn)+Δt2 d
(dx) dt
2! dt
x(tn+1) = x(tn) + Δtf(tn, xn) + Δt2 d f(t, x(t))
t=tn t=tn
Chain rule
2! dt x(tn+1)=x(tn)+Δtf(tn,xn)+Δt2 (∂f + ∂f dx)
2! ∂t ∂x dt
x(tn+1) = x(tn) + Δtf(tn, xn) + Δt2 ( ∂f (tn, xn) + ∂f f(tn, xn))
xn+1 = xn + Δtf(tn, xn) + Δt2 ( ∂f (tn, xn) + ∂f f(tn, xn)) 2!∂t ∂x
Taylors Method Order 2
t=tn 2!∂t ∂x
5
dx dt
= f(t,x)
xn+1 = xn + Δtf(tn, xn) + Δt2 ( ∂f (tn, xn) + f(tn, xn)) 2! ∂t
∂f ∂x
(18.1) Equation (18.1) is called the Taylor Method order 2. It is more accurate than the Euler’s method because
it has a local truncation error O (Δt3).
Taylor method order 2
Advantage – More accurate
Disadvantage – Need to calculate derivatives of f (t, x)
Example L18.1:
Using Taylors method order 2, solve
dx dt
= 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
Compare your solution with Euler’s method
=1−x Taylors Method Order 2 (Eq. (18.1))
2!
We need to find all terms in the equation above For this question
f(t, x) = 1 − x
dx dt
= f(t,x) xn+1 = xn + Δtf(tn, xn) + Δt2 (
+ f(tn, xn))
Therefore Eq. (18.1) becomes
x(tn + 1) = x(tn) + Δt Δt2 (1 − xn) 2!
∂f ∂t
=0
∂f (tn, xn) ∂t
∂f ∂x
=−1
(1−xn)
−
∂f ∂x
Explicit Euler’s method Taylor method order 2 Example L17.3 (last lecture) Example L18.1
xn+1 = xn + Δt(1 − xn) Lecture17M03.m
Lecture18M01.m
x(tn + 1) = x(tn) + Δt(1 − xn) − Δt2 (1 − xn) 2!
clear all;
close all; Delta_t=2; t=0:Delta_t:8.0;
x(1)=0.0;
for n=1:length(t)-1 x(n+1)=x(n)+Delta_t*(1-x(n))
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-‘,’linewidth’,2) xlabel(‘t’);
ylabel(‘x’);
legend(‘Euler’,’True’);
axis([0 8 0 2])
Only need to change one line of code!
clear all;
close all; Delta_t=2; t=0:Delta_t:8.0;
x(1)=0.0;
for n=1:length(t)-1 x(n+1)=x(n)+Delta_t*(1-x(n))-(Delta_t^2)*(1-x(n))/2;
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-‘,’linewidth’,2) xlabel(‘t’);
ylabel(‘x’);
legend(‘Taylor Method Order 2′,’True’); axis([0 8 0 2])
Lecture18M01.m
clear all;
close all; Delta_t=2; t=0:Delta_t:8.0;
x(1)=0.0;
for n=1:length(t)-1 x(n+1)=x(n)+Delta_t*(1-x(n))-(Delta_t^2)*(1-x(n))/2;
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-‘,’linewidth’,2) xlabel(‘t’);
ylabel(‘x’);
legend(‘Taylor Method Order 2′,’True’); axis([0 8 0 2])
Note: the Taylor Method order 2 for Δt = 2 predicts x = 0 for all values of t!
Output
Δt = 2.0
Δt = 2.0
Δt = 1.0
Δt = 0.5
Δt = 0.1
Explicit Euler
Δt = 1.0
Δt = 0.5
Taylor Method Order 2
Δt = 1.0
Δt = 0.5
For a given Δt, the Taylor Method Order 2 is more accurate than explicit Euler
End of Example L18.1
Example L18.2:
Use Taylors Method Order 2, solve
dx = 1 + x dt t
for t ∈ [1,2] with x(t = 0) = 2 and a)Δt = 0.2
b)Δt = 0.1
c) Δt = 0.5
d) Δt = 0.1
Compare your solution with Euler’s method and the solution given by MATLAB function ode23().
Taylors method order 2
xn+1 = xn + Δtf(tn, xn) + Δt2 ( + f(tn, xn)) 2!
xn+1 = xn + Δtf(tn, xn)
xn+1 = xn + Δt (1 + xn ) tn
dx = 1 + x dt t
= f(t,x)
Explicit Euler
We need to find all terms in the equation above
For this question
f ( t , x ) = 1 + xt Therefore Eq. (18.1) becomes
xn+1 = xn + Δt + Δt2 (− xn + 1 2! tn2 tn
xn+1 = xn + Δt (1 + xn ) + Δt2 ( 1 ) tn 2!tn
)
∂f = − x ∂t t2
∂f (tn, xn) ∂t
∂f ∂x
∂f = 1 ∂x t
(1 + xn ) tn
(1 + xn ) tn
Lecture18M02.m
t ∈ [1,2]
x(0) = 2
xEuler(n+1)=xEuler(n)+Deltat*exampleproblem(tEuler(n),xEuler(n));
xTMO2(n+1)=xTMO2(n)+Deltat*exampleproblem(tEuler(n),xTMO2(n))+((Deltat^2)/2.)*(1./tEuler(n));
xn+1 = xn + Δt (1 + xn ) + Δt2 ( 1 ) tn 2!tn
function dxdt=exampleproblem(t,x)
xn+1 = xn + Δt (1 + xn ) tn
clear all [t,x]=ode23(@exampleproblem,[1 2],2);
Deltat=0.2;
tEuler=1:Deltat:2; xEuler=2*ones(size(tEuler)); xTMO2=2*ones(size(tEuler));
for n=1:length(tEuler)-1 end
hold off
plot(t,x);
hold on plot(tEuler,xEuler,’ko’,’Markerfacecolor’,[0 0 0]); plot(tEuler,xTMO2,’ro’,’Markerfacecolor’,[1 0 0]); xlabel(‘t’,’Fontsize’,14)
ylabel(‘x(t)’,’Fontsize’,14)
dxdt=1+x/t;
end
Explicit Euler
dx = 1 + x dt t
Taylor Method Order 2
Δt = 0.2
Solution by MATLAB ode23()
Taylor Method Order 2 Explicit Euler
•
•
For Δt = 0.2, the solution obtained using Taylor Method Order 2 falls on top of the solution give by MATLAB function ode23().
For Δt = 0.2, the solution obtained using explicit Euler method falls just below the solution obtained using MATLAB function ode23().
L18.2: Stability Analysis: Explicit Euler Method
18
STABILITY OF THE EXPLICIT EULER METHOD
•
•
•
• •
In Lecture 17, we looked at the accuracy of the explicit Euler method. Let’s now look at the stability of the Euler method.
We know from previous lecture, the smaller Δt, the more accurate the solution, i.e. the truncation error go to zero. This feature is known as consistency.
When we make Δt big enough, the solution will “blow up”. Can we predict the value of Δt where the solution will “blow up”? For a numerical method to be useful, we need it to be stable.
Stability and consistency are two quite different things. It is possible for a method can be consistent but not stable.
If a numerical method is both consistent and stable, it is convergent.
•
= λx where λ is a constant and is a complex number
λRe + iλIm i = −1 is the imaginary unit.
To carry out stability analysis, we always consider the model problem
(18.2)
dx dt
Strategy
Compare
dx dt
= λx
Find analytical solution at time level t = nΔt
x(t) = x0eλt
= x0e(λRe+iλIm)t
Find approximated Euler’s solution at time level tn = nΔt.
The analytical solution for Eq. (18.2) can be obtained by writing dx = λdt. Integrating both sides and using the initial condition x(t = 0) = x0 gives x
x(t) = x0eλt
= x0e(λRe+iλIm)t
= x0eλReteiλImt
x(t)=xeλ t(cos(λ t)+isin(λ t)) 0 Re Im Im
λIm
x(t)
x(t)
x(t)
t
x(t)
t
t
t
x(t)
t
λRe
x(t) = x0eλRet(cos(λImt) + i sin(λImt)) 22
x(t) = x0eλRet(cos(λImt) + i sin(λImt))
23
So the real (true) solution is stable as long as λRe ≤ 0
Find analytical solution at time level tn = nΔt
x(t) = x0eλRet(cos(λImt) + i sin(λImt)) Compare
dx dt
= λx
Applying Euler’s formula for the model problem (Eq. (18.2)) gives
xn+1 = xn + Δtf(tn, xn) For the model problem dx/dt = λx, f(t, x) = λx. So
xn+1 = xn + Δtλxn = (1 + λΔt)xn
x = σx n+1 n
where σ = (1 + λΔt) is the amplification factor
(18.3)
Find approximated Euler’s solution at time level tn = nΔt
xn+1 = σxn Thus, running a program using the Euler’s method will give you
x0 = Given
x1 = σx0 2
x2 = σx1 = σ x0 x3 = σx2 = σ3×0
⋮=⋮=⋮ xn+1 = σxn = x0
(18.4)
σn
For solution to be stable
|σ| ≤ 1
= 1+ t
= | |2 =
(1 + t Re + i t Im)
(1 + t Re)2 + ( t Im)2 1
Im t
-2
-1
Re t
x(t) = x0eλRet(cos(λImt) + i sin(λImt))
Find analytical solution at time level tn = nΔt
Find approximated Euler’s solution at time level tn = nΔt
xn+1 = σnx0
σ = (1 + λΔt)
Stability region for analytical solution
-2
-1
27
dx dt
= λx
Compare
Im t
Re t
Example L18.3:
Write a MATLAB program that uses explicit Euler’s method to solve
dx dt
for 0 ≤ t ≤ 2 with x(t = 0) = 1. Using the analysis from the previous slides, what are the values of Δt for the numerical solution to be stable?
= − 8x
We have
-8Δt, Δt large -8Δt
dx dt
= − 8x
So for this problem λ = − 8. So λRe = − 8 and λIm = 0. In order for the numerical solution to be
stable, λΔt need to be within the stable region.
Im t
-2
-1
Re t
2 t 0 2 8 t 0 0 t 1/4
So for solution to
be stable Δt must
be between 0 and 0.25
Im t
-2
-1
Re t
Lecture18M03.m
clear all; close all;
Delta_t=0.5; t=0:Delta_t:2.0;
%Preallocating memory
x=zeros(size(t));
xn+1 = xn + Δt x(n+1)=x(n)+Delta_t*(-8*x(n));
end
plot(t,x,’ko-‘) hold on
tplot=0:0.01:2.0; truesolution=@(t)(exp(-8*t)) plot(tplot,truesolution(tplot),’b-‘,’linewidth’,2) axis([0 2 -2 2])
xlabel(‘t’); ylabel(‘x(t)’); legend(‘Euler’,’True’);
x(1)=1.0;
for n=1:length(t)-1
dx dt
= − 8x
(−8xn)
Δt = 0.5
Δt = 0.25
Δt = 0.1
Δt = 0.01
Euler solution is unstable for Δt = 0.5
Euler solution is neutrally stable for Δt = 0.25
Euler solution is stable for Δt = 0.1 & Δt = 0.01
Im t
-4
-2
-1 Re t
Δt = 0.5
t = 8 ⇥ 0.5 = 4
Im t
-2
Δt = 0.25
-1
Re t
t = 8 ⇥ 0.25 = 2
-2
-1
Re t
Im t
Δt = 0.1
t = 8 ⇥ 0.1 = 0.8
Im t
-2
Δt = 0.01
-1
Re t
t = 8 ⇥ 0.01 = 0.08
End of Example L18.3
Consider now the case where λ is purely imaginary,
= Re + i Im
dx = λx 0 dx = iλImx dt dt
For this case, iλImΔt, is always on the vertical axis of the stability diagram. It is not within the stability region. Hence, in this case, explicit Euler method is always unstable.
Im t
–
–
Re t
Example L18.4:
Write a Matlab program that uses Euler’s method to solve
dx dt
for t ∈ [0,10] with x(t = 0) = 1. Using the analysis from the previous slides, what are the values of Δt for the numerical solution to be stable?
= i2x
dx dt
= i2x x = Aei2t
x = (ARe + iAIm)ei2t
x = (ARe + iAIm)(cos(2t) + i sin(2t))
x = (ARe cos(2t) − AIm sin(2t)) + i(AIm cos(2t) + ARe sin(2t)) But we are given that x = 1 at t = 0
1 + i0 = ARe + iAIm Hence
ARe=1 AIm=0
0 ≤ t ≤ 10
The analytical solution for this problem can be written as
x(t = 0) = 1
dx dt
= i2x
So the analytical (true) solution for this equation is
x(t) = ei2t
= cos(2t) + i sin(2t)
R𝔢(x) = cos(2t) I𝔪(x) = sin(2t)
If we use explicit Euler method to find the approximate the solution
xn+1 = xn + Δtf(tn, xn) For this problem f (tn, xn) = i2xn. So
xn+1 = xn + Δti2xn
0 ≤ t ≤ 10 x = (ARe + iAIm)ei2t
x(t = 0) = 1 ARe = 1
AIm = 0
Lecture18M04.m
clear all; close all;
Delta_t=0.5;
%Preallocating memory
x=zeros(size(t));
for n=1:length(t)-1 end
dx dt
= i2x
0 ≤ t ≤ 10
x(t = 0) = 1 x(n+1)=x(n)+Delta_t*(i*2*x(n));
t=0:Delta_t:10.0;
x(1)=1.0;
xn+1 = xn + Δti2xn
plot(t,imag(x),’ko-‘,’linewidth’,2,’markersize’,10) hold on
tplot=0:0.01:10.0; truesolution=@(t)(sin(2*t)) plot(tplot,truesolution(tplot),’b-‘,’linewidth’,2) axis([0 10 -5 5])
xlabel(‘t’); ylabel(‘x(t)’); legend(‘Euler’,’True’);
Δt = 0.5
I𝔪(x) Δt = 0.5
Δt = 0.2
Δt = 0.1
Δt = 0.02
Im t i
-2 -1
λΔt = i2 × 0.5 = i Re t
Δt = 0.5
Im t
0.4i -2 -1
λΔt = i2 × 0.2 = 0.4i Re t
Δt = 0.2
Im t
-2 -1
λΔt = i2 × 0.1 = 0.2i 0.2i
Re t
Δt = 0.1
Im t
-2 -1
λΔt = i2 × 0.02 = 0.04i 0.04i
Re t
Δt = 0.02
For all values of Δt, λΔt is outside the stable region. Hence the numerical solution obtained using Euler’s method will never be stable for any value of Δt.
End of Example L18.4
L18.3: Stability Analysis: Taylor Method Order 2
50
Recall that the Taylor method order 2 formula is
x(tn+1)=x(tn)+Δtf(tn,xn)+Δt2 (∂f + f(tn,xn)) 2! ∂t
For stability analysis, we use the formula above to obtain solution to the model equation,
For this question
dx =λx=f(t,x) dt
f(t,x) = λx
Therefore formula for Taylor Method Order 2 becomes
xn+1 = xn + Δtλxn + Δt2 (λλxn) 2!
xn+1 = 1+λΔt+ (λΔt) xn ( 2! 2)
xn+1 = σxn
∂f ∂t
=0
∂f ∂x
∂f ∂x
=λ
where σ = 1+λΔt+ (λΔt)
( 2! 2)
So
where
For stability, we need
xn+1 = σxn
σ = 1+λΔt+ (λΔt)
2! 2)
Remember that λ = λRe + iλIm is an imaginary variable, so σ is a complex quantity with real
(
and imaginary parts.
We can use MATLAB to plot contours of
|σ| = |1 + λΔt + (λΔt)2/2| and identify regions where |σ| ≤ 1.
|σ| ≤ 1
52
(-4,4) (4,4)
(-4,-4) (4,-4)
Lecture18M05.m
|σ| = 4 |σ| = 3 |σ| = 2 |σ| = 1
53
relamdt=-4:0.1:4;
imlamdt=-4:0.1:4; [RElamdt,IMlamdt]=meshgrid(relamdt,imlamdt); axis square;
lamdt=RElamdt+i*IMlamdt;
λΔt = λReΔt + iλImΔt |σ| = |1 + λΔt + (λΔt)2/2|
sig=1+lamdt+lamdt.^2/2;
contour_levels=[1 2 3 4]; contour(RElamdt,IMlamdt,abs(sig),contour_levels,’linewidth’,4); xlabel(‘Re \lambda \Delta t’,’Fontsize’,14);
ylabel(‘Im \lambda \Delta t’,’Fontsize’,14);
hold on;
plot([-4 4],[0 0],’k’);
plot([0 0],[-4 4],’k’);
title(‘contour |\sigma|’,’Fontsize’,14);
clear all; close all;
|σ| = 4 |σ| = 3
|σ| = 2 |σ| = 1
Stability region of Taylor Method order 2
54
Model problem dx = λx dt
Stability region analytical solution
x(t) = x0eλRet(cos(λImt) + i sin(λImt))
Stability region Taylor Method Order 2 n
xn+1 = σ x0 where
σ = 1+λΔt+ (λΔt)
( 2! 2)
Stability region explicit Euler n
xn+1 = σ x0 where
σ = 1 + λΔt
55
Model problem dx = λx dt
Stability region analytical solution
x(t) = x0eλRet(cos(λImt) + i sin(λImt))
Stability region explicit Euler n
xn+1 = σ x0 where
σ = 1 + λΔt
Stability region Taylor Method Order 2 xn+1 = σnx0 where
σ = 1+λΔt+ (λΔt)
( 2! 2)
56
Model problem dx = λx dt
Stability region analytical solution
x(t) = x0eλRet(cos(λImt) + i sin(λImt))
Stability region explicit Euler n
xn+1 = σ x0 where
σ = 1 + λΔt
Stability region Taylor Method Order 2 n
xn+1 = σ x0 where
σ = 1+λΔt+ (λΔt)
( 2! 2)
57