Numerical Methods in Engineering (ENGR20005)
Lecture 20 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L20.1:
Ordinary Differential Equations:
Nonlinear problems
2
Example L20.1: Use explicit Euler’s method to solve the ODE
dx =x(x−x2)
5 dt −5
in the domain 0 ≤ t ≤ 2 × 10 . You are given that x(0) = 10 . Compare your answer with the answer given by MATLAB function ode23().
This is a nonlinear problem because dx/dt = f(t, x) is a nonlinear function. It does not have an analytical (exact) solution. We will compare our solution to the solution given by the inbuilt MATLAB functions.
Applying the explicit Euler equation gives
xn+1 = xn + Δtf(tn, xn)
xn+1 = xn + Δt
Computing the explicit Euler solution is straight forward using the above equation
x (x −x2) nnn
Lecture20M01.m
x(n+1)=x(n)+f(t(n),x(n))*Delta_t;
xn+1 = xn + Δtf(tn, xn)
f(t, x) = x(x − x2)
x(0) = 10−5
0 ≤ t ≤ 2 × 105
dx dt
= x (x − x2)
dxdt=x*(x-x^2);
clear all; close all;
epsilon=0.00001;
Delta_t=1; t=0:Delta_t:2.0/epsilon;
%Preallocating memory
x=zeros(size(t)); x(1)=epsilon;
for n=1:length(t)-1
end
%
%Computing Matlab solution
%
[tmat,xmat]=ode23(@f,[0 2/epsilon],epsilon);
plot(t,x,’ko-‘,tmat,xmat,’bs-‘) hold on
xlabel(‘t’);
ylabel(‘x’);
legend(‘Explicit Euler’,’MATLAB’); function dxdt=f(t,x)
end
Δt = 1 • The solution stays close to 5
zero until about t = 1 × 10 .
• At t = 1 × 105, there is a sharp change and x quickly increases to 1.
• At t > 1 × 105, the value of x quickly stays at 1.
Some characteristics of the solution
plot(t,x,’ko-‘,tmat,xmat,’bs-‘) 0 hold on
xlabel(‘t’);
ylabel(‘x’);
legend(‘Explicit Euler’,’MATLAB’);
≤ t ≤ 2 × 105
dx dt
= x (x − x2)
dxdt=x*(x-x^2);
function dxdt=f(t,x) end
Example L20.2: Use implicit Euler method to solve the nonlinear ODE
dx
= x (x − x2)
in the domain 0 ≤ t ≤ 2 × 105. You are given that x(0) = 10−5.
dt
The last example (Example L20.1) shows that applying the explicit Euler equation gives
x =x +Δtx (x −x2) n+1 n n n n
At any time level, we know the values of all the terms on the right hand side of this equation. So computing the explicit Euler solution is therefore straight forward using the above equation (see Lecture20M01.m)
Applying the implicit Euler’s formula
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
gives
xn+1 = xn + Δt
f(tn+1, xn+1)
Example L20.2: Use implicit Euler method to solve the nonlinear ODE
dx
= x (x − x2)
in the domain 0 ≤ t ≤ 2 × 105. You are given that x(0) = 10−5.
dt
The last example shows that applying the explicit Euler equation gives
x =x +Δtx (x −x2) n+1 n n n n
At any time level, we know the values of all the terms on the right hand side of this equation. So computing the explicit Euler solution is therefore straight forward using the above equation (see Lecture20M01.m)
Applying the implicit Euler’s formula
x(tn+1) = x(tn) + Δtf(tn+1, xn+1) xn+1 = xn + Δt
gives
xn+1 (xn+1 − x2 ) n+1
xn+1=xn+Δtxn+1(xn+1−x2 ) n+1
Rearranging gives
xn+1−Δt(x2 −x3 )−xn=0 n+1 n+1
For every time step, we typically know xn. We need to find xn+1. So we need to solve a nonlinear equation xn+1−Δt(x2 −x3 )−xn=0
n+1 n+1
every single time step to get xn+1. Details of the steps one needs to take is given on the next slide.
t0=a t1 = a + Δt
t2 = a + 2Δt t3 = a + 3Δt
xx x3 x2
x0 1
x0=Given
x1 − Δt (x2 − x3) − x0 = 0
xn+1−Δt(x2 −x3 )−xn=0 n+1 n+1
Solve to get
Solve to get
11
x1 x2
x2 − Δt (x2 − x3) − x1 = 0 2 2
x3−Δt x2−x3 −x2=0 Solve 3 () togetx
33
t0 =a
t1 =a+Δt
t2 =a+2Δt
t3 =a+3Δt t
( 2 3) known x3−Δt x −x −x2=0
We know x2. We want to find x3
So every time step, we have to solve an equation that looks something like
33
Cubic polynomial in terms of p
2 3 known p−Δt(p −p )−K=0
g(p) = p − Δt(p2 − p3) − K = 0
p is the unknown that you need to find. How do we solve for p? This is a root finding problem.
We need to find p such that g(p) = 0. In this example, I am going to show you how to use Newton- Raphson method but you can use any of the methods taught earlier in the semester.
Let’s now first look at the root finding problem (Example L20.3) and see how we can incorporate it to time step forward the implicit Euler method for solving
dx dt
which is Example L20.2
= x (x − x2)
Example L20.3: Write a MATLAB function that uses the Newton-Raphson method to solve
g(p) = p − 0.1(p2 − p3) − 10 = 0
(p−0.1 (p2−p3)−10)
6 4 2 0
g( p)
−4
−6
−8 −10
−2
The answer should be around p = 4.2
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
p
Newton-Raphson Method
pnew = pold − For this problem
g(p) = p − 0.1(p2 − p3) − 10 dg/dp(p) = 1 − 0.1(2p − 3p2)
g(pold)
dg/dp(pold)
function Lecture20M03()
clear all;
close all;
%initial guess value of p
p=3;
gp=p-0.1*(p^2-p^3)-10;
dgdp=1-0.1*(2*p-3*p^2);
while abs(gp)>1.0e-13
p=p-gp/dgdp;
gp=p-0.1*(p^2-p^3)-10;
dgdp=1-0.1*(2*p-3*p^2);
end
p
g(p) = p − 0.1(p2 − p3) − 10 dg/dp(p) = 1 − 0.1(2p − 3p2
g(pold) dg/dp(pold)
pnew = pold − 6
4
2
g( p) −2
−4
0
−63 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5
p
)
End of Example L20.3…. ..now back to Example L20.2
Recall that we want to solve
dx dt
= x (x − x2) Using the implicit Euler method. We ended up having to solve
xn+1−Δt(x2 −x3 )−xn=0 n+1 n+1
for xn+1 every time step. It is important to note that xn is a known value at every time step.
t0=a t1 = a + Δt
t2 = a + 2Δt t3 = a + 3Δt
xx x3 x2
x0 1
x0=Given
x1 − Δt (x2 − x3) − x0 = 0
xn+1−Δt(x2 −x3 )−xn=0 n+1 n+1
Solve to get
Solve to get
11
x1 x2
x2 − Δt (x2 − x3) − x1 = 0 2 2
x3−Δt x2−x3 −x2=0 Solve 3 () togetx
33
t0 =a
t1 =a+Δt
t2 =a+2Δt
t3 =a+3Δt t
Newton Raphson loop
xn xn+1
Hence, for every time step, we need to solve an equation that looks like
xn+1−Δt(x2 −x3 )−xn=0 (20.1) n+1 n+1
To use the Newton-Raphson formula, we need to find its derivative with respect to the variable which you are trying to solve for, which in this case is xn+1. For the equation above, this is
1 − Δt(2xn+1 − 3×2 ) n+1
Using the notation in Example L20.3, we let p = xn+1 and Eq. (20.1) becomes g(p) = p − Δt(p2 − p3) − xn
dg(p)/dp = 1 − Δt(2p − 3p2)
So given Δt and xn we should be able to use the Newton Raphson formula
pnew = pold − to find a value of p such that g(p) = 0
g(pold)
dg/dp(pold)
Newton-Raphson formula Lecture20M03.m
pnew = pold
Lecture20M02.m
g(pold) dg (pold)
dp
g(p) = p − 0.1(p2 − p3) − 10 dg/dp(p) = 1 − 0.1(2p − 3p2)
g(p) = p − Δt(p2 − p3) − xn dg(p)/dp = 1 − Δt(2p − 3p2)
function p=Findxnplus1( ,xn)
clear all;
close all;
%initial guess value of p
p=xn;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
while abs(gp)>1.0e-15
p=p-gp/dgdp;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
end
Delta_t
function Lecture20M03()
clear all;
close all;
%initial guess value of p
p=3;
gp=p-0.1*(p^2-p^3)-10;
dgdp=1-0.1*(2*p-3*p^2);
while abs(gp)>1.0e-13
p=p-gp/dgdp;
gp=p-0.1*(p^2-p^3)-10;
dgdp=1-0.1*(2*p-3*p^2);
end
function p=Findxnplus1(Delta_t,xn)
clear all;
close all;
%initial guess value of p
p=xn;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
while abs(gp)>1.0e-15
p=p-gp/dgdp;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
end
xn xn+1
clear all;
close all;
epsilon=0.00001;
Delta_t=10;
t=0:Delta_t:2.0/epsilon;
%Preallocating memory
x=zeros(size(t));
x(1)=epsilon;
for n=1:length(t)-1
x(n+1)=Findxnplus1(Delta_t,x(n));
end
plot(t,x);
axis([0 2/epsilon -0.1 1.1])
function p=Findxnplus1(Delta_t,xn)
%initial guess value of p
p=xn;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
while abs(gp)>1.0e-15
p=p-gp/dgdp;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
end end
Newton Raphson loop
clear all;
close all;
epsilon=0.00001;
Delta_t=10;
t=0:Delta_t:2.0/epsilon;
%Preallocating memory
x=zeros(size(t));
x(1)=epsilon;
for n=1:length(t)-1
x(n+1)=Findxnplus1(Delta_t,x(n));
end
plot(t,x);
axis([0 2/epsilon -0.1 1.1])
function p=Findxnplus1(Delta_t,xn)
%initial guess value of p
p=xn;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
while abs(gp)>1.0e-15
p=p-gp/dgdp;
gp=p-Delta_t*(p^2-p^3)-xn;
dgdp=1-Delta_t*(2*p-3*p^2);
end end
g ( p ) = p t ( p 2 p 3 ) x l n2
dg(p)/dp = 1 t(2p 3p )
pnew = pold g(pold) dg/dp(pold)
Implicit Euler solution for
dx dt
= x (x − x2) for Δt = 10.
L20.2:
Some Comments on Methods for Solving Ordinary Differential Equations
24
Methods for solving initial value problems
dx dt
for xn+1 can be calculated explicitly, without further iterations or solving another equation. Implicit methods are methods where value of xn+1 is defined implicitly and further calculations/iterations are
needed to calculate the value of xn+1, especially when f (t, x) is a nonlinear function. Examples of explicit methods are explicit Euler and Runge-Kutta methods. Examples of implicit methods are the
= f(t,x) t ∈ [a,b] x(a) = α
can be classified into implicit or explicit methods. Explicit methods are methods where the value
implicit Euler and the Crank-Nicolson methods (see later).
Methods for solving initial value problems.
= f(t,x)
dx dt
Explicit methods
Implicit methods
Runge-Kutta methods
k1 =f(tn,xn)
k2 = f(tn + Δt, xn + k1Δt)
xn+1 = xn + Δt (k1 + k2) 2
Explicit Euler
xn+1 = xn + f(tn, xn)
Implicit Euler
xn+1 = xn + f(tn+1, xn+1)
Crank-Nicolson
xn+1 = xn + Δt [f(tn, xn) + f(tn+1, xn+1)] 2
When you have a sharp
change in the solution, the implicit
methods tend to perform much better
than the explicit methods. MATLAB
function ode23() is an explicit
method and ode23s() is an implicit
method. Note that there is no oscillation
in the solution predicted by ode23s() but
the solution predicted by ode23() oscillates 5 in the vicinity of the discontinuity at t = 1 × 10 .
dx dt
x(0) = 10−5
0 ≤ t ≤ 2 × 105
= x (x − x2)
Implicit Euler
Matlab (ode23s())
Explicit Euler
MATLAB (ode23())
x
x
Comparing implicit and explicit Euler solution
1.06
1.04
1.02
1
0.98
0.96
1.005
1
0.995
0.99
0.98 1
1.02 1.04
1.06
1.08
1.1 1.12
0.995
1
1.005
1.01
1.015 x 105
Explicit Euler and MATLAB ode23()
Implicit Euler and MATLAB ode23s()
t
x 105
t
• For a given Δt, explicit methods are usually less stable than implicit method • Near a sharp gradient, explicit methods can give oscillatory (wrong) solution
MATLAB function ode23() and ode23s() are explicit and implicit methods respectively
φ
φ
Summar y
The solution using implicit Euler method is much more difficult to compute than the explicit Euler method, especially for nonlinear problems (need to use root finding methods).
Implicit Euler method is much more stable than the explicit Euler method
It is possible to obtain a solution with much larger Δt using the implicit Euler method than the explicit Euler method. For example L20.2, for explicit Euler method, the solution blows up for Δt = 10 whereas the solution using implicit Euler method is stable for Δt = 10.
L20.3: Crank-Nicolson method
29
CRANK-NICOLSON FORMULA
To solve:
= f(t,x) Can use the Crank-Nicolson formula:
dx
dt
xn+1 = xn + Δt [f(tn, xn) + f(tn+1, xn+1)] 2
Note that the gradient f (t, x) is calculated at time levels tn and tn+1 and then the average is taken. Since this f (t, x) is needed at time level tn+1, the Crank-Nicolson method is an implicit method.
Will first illustrate why this works graphically.
x(t)
dx =f(tn,xn)
x(tn+1) x(tn)
Error
Δtdx =Δtf(tn,xn) dt
Explicit Euler Predicted value of x(tn+1)
dt
x(t)
True value of x(tn+1)
xn+1 = xn + Δt f (tn, xn) Explicit Euler’s formula
Δt
tn
tn+1
x(t) x(tn+1)
x(tn)
Explicit Euler x(t) Predicted value of x(tn+1)
Implicit Euler
dx
dt = f(tn, xn)
Error Δtdx =Δtf(tn,xn) dt
x(t)
True value of x(tn+1) x(tn+1)
True value x(tn+1) x(t)
dx =f(tn+1,xn+1) dt
x(tn+1) = x(tn) + Δtf(tn+1, xn+1) Implicit Euler’s formula
tn+1
xn
xn+1 = xn + Δtf(tn, xn) Explicit Euler’s formula
Δt
Δt
tn
tn+1
tn
x(t) x(tn+1)
x(tn)
tn
Explicit Euler Predicted value of x(tn+1)
x(t)
x(tn+1) xn
Implicit Euler
dx
dt = f(tn, xn)
Error
Δtdx =Δtf(tn,xn) dt
x(t)
True value of x(tn+1)
True value x(tn+1) x(t)
dx =f(tn+1,xn+1) dt
Predicted xn+1
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Implicit Euler’s formula
tn+1
xn+1 = xn + Δtf(tn, xn) Explicit Euler’s formula
Δt
Δt
tn+1
tn
x(t) x(tn+1)
x(tn)
tn
Explicit Euler Predicted value of x(tn+1)
x(t)
x(tn+1) xn
Implicit Euler
dx
dt = f(tn, xn)
Error
Δtdx =Δtf(tn,xn) dt
x(t)
True value of x(tn+1)
True value x(tn+1) x(t)
dx =f(tn+1,xn+1) dt
Predicted xn+1
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Implicit Euler’s formula
Error
xn+1 = xn + Δtf(tn, xn) Explicit Euler’s formula
Δt
x(t)
Crank-Nicolson
tn+1
Δt xn+1 = xn + Δ2t [f(tn, xn) + f(tn+1, xn+1)]
tn+1
tn True value x(tn+1)
xn
Δt
tn
tn+1
x(t) x(tn+1)
x(tn)
tn
Explicit Euler Predicted value of x(tn+1)
x(t)
x(tn+1) xn
Implicit Euler
dx
dt = f(tn, xn)
Error
Δtdx =Δtf(tn,xn) dt
x(t)
True value of x(tn+1)
True value x(tn+1) x(t)
dx =f(tn+1,xn+1) dt
Predicted xn+1
x(tn+1) = x(tn) + Δtf(tn+1, xn+1)
Implicit Euler’s formula
Error
Δt
x(t)
Crank-Nicolson
tn+1
Δt xn+1 = xn + Δ2t [f(tn, xn) + f(tn+1, xn+1)]
tn+1
xn+1 = xn + Δtf(tn, xn) Explicit Euler’s formula
tn True value x(tn+1)
xn
Δt
Error Predicted xn+1
tn
tn+1
Example L20.4:
Using Crank-Nicolson method, 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 explicit Euler’s method
Crank Nicolson formula
dx dt
=1−x =f(t,x)
xn+1 = xn + Δt 12 [f(tn, xn) + f(tn+1, xn+1)]
xn+1 =xn+Δt12[1−xn+1−xn+1] This is a linear equation so it is easy to solve for xn+1
(1 + Δt ) xn+1 = (1 − Δt ) xn + Δt 2 ( Δt) 2
x =1−2 xn+Δt
n+1 (1 + Δt ) 2
Lecture20M04.m
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)=((1-Delta_t/2)*x(n)+Delta_t)/(1+Delta_t/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(‘Crank-Nicolson’,’True’);
axis([0 8 0 2])
(1 − Δt ) xn + Δt 1 + Δt
xn+1 =
2
(2)
Compare explicit Euler and method for Δt = 2.0
implicit Euler,
Crank Nicolson
Δt = 2.0
Δt = 2.0
Δt = 2.0
End of Example L20.4
Example L20.4x: Use Crank-Nicolson method to solve the nonlinear ODE
dx
= x (x − x2)
in the domain 0 ≤ t ≤ 2 × 105. You are given that x(0) = 10−5.
dt
xn+1 = xn + Δt [f(tn, xn) + f(tn+1, xn+1)] 2
Applying the Crank-Nicolson equation gives
x =x +Δt/2(x(x −x2)+x (x −x2 ))
n+1 n n n n n+1 n+1 So given Δt and xn we should be able to use the Newton Raphson formula
n+1
We have n+1 terms on both sides of the equations and can’t rearrange easily. Need to use one of our
root finding methods.
we let p = xn+1
dg(p)/dp = Δt/2 ((p − p2) + p(1 − 2p)) − 1
So given Δt and xn we should be able to use the Newton Raphson formula pnew = pold −
g(p)=x +Δt/2(x(x −x2)+p(p−p2))−p nnnn
Now, we should be able to code this up…..
g(pold)
dg/dp(pold)
ACCURACY OF CRANK-NICOLSON
Going to show that the Crank-Nicolson time-marching method is 2nd-order accurate.
xn+1=xn+Δt1 dx +dx +Δt3(…) 2dtn dtn+1
That is, that the local truncation error is proportional to Δt3, or that the global truncation error is proportional to Δt2.
ACCURACY OF CRANK-NICOLSON Taylor-series expand x about t = tn:
xn+1=xn+Δtdx +Δt21d2x +Δt3(…) dt n 2 dt2 n
(20.2)
Then Taylor-series expand x about t = tn+1:
xn = xn+1 − Δt dx + Δt2 1 d2x
+ Δt3(…) +dx +Δt21 d2x −d2x
(20.3)
dt n+1
2 dt2 n+1
dt n+1 2 dt2 n dt2 n+1
Take Eq. (20.2)-(20.3)
xn+1−xn =xn−xn+1+Δt dx dt n
+Δt3(…)
xn+1−xn =xn−xn+1+Δt dx +dx +Δt21 d2x −d2x +Δt3(…)
dt n
Divide the above equation by 2 to give
dt n+1 2 dt2 n dt2 n+1
xn+1=xn+Δt1 dx +dx 2dtn dtn+1
Crank-Nicolson
+Δt2 1 d2x − d2x + Δt3(…) 4 dt2 n dt2 n+1
Error
ACCURACY OF CRANK-NICOLSON
xn+1 = xn + Δt 1 dx
2 dt n
+Δt2 14
Taylor-series expand d2x/dt2 about t = tn: d2x = d2x + Δt d3x
Let’s further analyse thisterm
(20.4)
+ dx
dt n+1
+ Δt3(…)
+ Δt2 d4x + . . . . . . 2!dt4 n
d2x − d2x
dt2 n dt2 n+1
dt2 n+1 dt2 n Substitute back into Eq. (20.4) to get
dt3 n
d2x −d2x =Δtd3x +Δt2 d4x +…… dt2 n+1 dt2 n dt3 n 2! dt4 n
xn+1=xn+Δt1 dx +dx 2dtn dtn+1
+Δt3(…)
xn+1=xn+Δt1 dx +dx 2dtn dtn+1
3
+ Δt (…) + t (…)
Small
Crank-Nicolson formula
So the local truncation error for Crank-Nicolson formula is O (Δt3) and the global truncation error
is O (Δt2).
STABILITY OF CRANK-NICOLSON
Consider the model problem
dx dt
= λx Apply the Crank-Nicolson method to obtain
x n + 1 = x n + Δ t 12 [ λ x n + λ x n + 1 ]
(1 − λΔt/2)xn+1 = (1 + λΔt/2)xn x = 1 + λΔt/2 x
n+1 1−λΔt/2 n
= σxn = σxn
STABILITY OF CRANK-NICOLSON
xn+1 = 1+λΔt/2xn 1 − λΔt/2
= σxn
For stability (solution not to blow up) we need |σ| ≤ 1.
Crank-Nicolson:
That is,
1 + λΔt/2 1 − λΔt/2
σ=
|σ|= 1+λΔt/2 ≤1
1 − λΔt/2 1 + λReΔt/2 + iλImΔt/2
≤1 (1+λReΔt/2)2+(λImΔt/2)2 ≤1
1 − λReΔt/2 − iλImΔt/2 (1 − λReΔt/2)2 + (λImΔt/2)2
(1 + λReΔt/2)2 ≤ (1 − λReΔt/2)2 λReΔt ≤ 0
STABILITY OF CRANK-NICOLSON
For stability (solution not to blow up) we need |σ| ≤ 1. That is,
Im t
Re t
The region of stability (|σ| ≤ 1) for the Crank- Nicolson method is when λReΔt ≤ 0 (left half-plane).
Recall that the stability region for the analytical solution is also when λRe < 0. This is very similar to the stability region of the Crank-Nicolson scheme which is when λReΔt < 0. For comparison, the stability region of the and implicit Euler are also shown on the diagram above. The stability region for is much smaller and the stability region for implicit Euler is much Im t -2 -1 12 Re t explicit Euler explicit Euler bigger than the analytical solution. Let’s now consider the case where λ is purely imaginary, λ = iω where ω is a real number. The analytical solution is x(t) = eiωt. It oscillates with magnitude of oscillation 1. Explicit Euler xn+1 = xn + iωΔtxn xn+1 = (1 + iωΔt)xn x = σx n+1 n σ = (1 + iωΔt) |σ| = 1 + (ωΔt)2 ≥1 Implicit Euler xn+1 = xn + iωΔtxn+1 (1 − iωΔt)xn+1 = xn Crank-Nicolson xn+1 = xn + Δt 12 [iωxn + iωxn+1] (1 − iωΔt/2)xn+1 = (1 + iωΔt/2)xn ( σ = ( 1 ) xn+1 = 1 + iωΔt/2 xn 1−iωΔt/2 1 ) 1 − iωΔt xn dx dt = iωx xn+1 = xn+1 = σxn = σxn σ = 1 + iωΔt/2 1−iωΔt |σ| = 1 1 − iωΔt/2 |σ|= 1+(ωΔt/2)2 1+(ωΔt/2)2 |σ|=1 ≤1 1+(ωΔt)2 So even though the Crank-Nicolson solution has an error associated with it, the behaviour of the approximation is such the amplitude of oscillation is exactly the same as the analytical solution. The solution given by explicit and implicit Euler grows and decays with time respectively. Example L20.5: Write a MATLAB program that uses the Crank-Nicolson method to solve = i2x for t ∈ [0,10] with x(t = 0) = 1. dx dt Analytical solution: Applying the Crank-Nicolson method, we obtain x(t) = ei2t xn+1 = xn + Δt 12 [i2xn + i2xn+1] (1 − iΔt)xn+1 = (1 + iΔt)xn x = 1 + iΔt x n+1 1−iΔt n Lecture20M05.m xn+1 = 1+iΔtxn 1−iΔt clear all; close all; Delta_t=0.1; t=0:Delta_t:10.0; %Preallocating memory x=zeros(size(t)); %Initial condition x(1)=1.0; for n=1:length(t)-1 x(n+1)=x(n)*(1+Delta_t*i)/(1-Delta_t*i); end plot(t,real(x),'ko-','linewidth',2,'Markersize',10) hold on tplot=0:0.01:10.0; truesolution=@(t)(real(exp(i*2*t))) plot(tplot,truesolution(tplot),'b-','linewidth',2) xlabel('t'); ylabel('x'); legend('Crank-Nicolson','True'); axis([0 10 -1 1]) 4 Δt = 0.2 Explicit Euler Implicit Euler Crank-Nicolson true solution 3 2 1 0 -1 -2 -3 -4 0 1 2 3 4 5 6 7 8 9 10 t explicit Euler • Amplitude of solution grows with time • Amplitude of implicit Euler solution decays with time • Amplitude of Crank-Nicolson solution remain constant x End of Example L20.5