Numerical Methods in Engineering (ENGR20005)
Lecture 22 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L22.1:
Implicit Euler’s method for a system of linear equations
Implicit Euler’s method for a system of equations
So far, we have learnt to solve single ODE using implicit Euler method. What if you are required a system of M + 1 ODEs?
dx0 dt
dx1 dt
dx2 dt
dxM dt
= f0(t, x0, x1, x2, ……, xM)
= f1(t, x0, x1, x2, ……, xM)
= f2(t, x0, x1, x2, ……, xM) ⋮=⋮
= fM(t, x0, x1, x2, ……, xM)
Implicit Eulers method can easily be extended to solve a system of M + 1 equations
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) IMPLICIT EULER
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) IMPLICIT EULER
dx0 dt
dx1 dt
= f0(t, x0, x1) = f1(t, x0, x1)
IMPLICIT EULER
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1)Δt
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) IMPLICIT EULER
dx0 dt
dx1 dt
dx2 dt
= f0(t, x0, x1, x2) = f1(t, x0, x1, x2) = f2(t, x0, x1, x2)
IMPLICIT EULER
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt xn+1,2 = xn,2 + f2(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt
dx0 dt
dx1 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1)Δt
= f0(t, x0, x1) = f1(t, x0, x1)
IMPLICIT EULER
I showed in Lecture 21 that if {f} are linear functions, then the system of equations can be put in matrix form as
dx0 dt
dx1 dt
dx2 dt
= f0(t, x0, x1, x2) = f1(t, x0, x1, x2) = f2(t, x0, x1, x2)
d dt
x0 K11 K12 K13 x0 x1 = K21 K22 K23 x1
x2 K31 K32 K33 x2
d dt
=
x0 x1 x2
K11 K12 K13 K21 K22 K23 K31 K32 K33
x0 x1 x2
{x} = Applying implicit Euler’s formula gives
d dt
[K ]{x}
{xn+1} = {xn} + [K ]{xn+1} Δt [I − [K]Δt]{xn+1} = {xn}
{xn+1} = [I − [K]Δt]−1{xn}
Example L22.1
Use Implicit Euler’s method to find the solution to the following
two ODEs
= 5×0 + 2×1
for 0 ≤ t ≤ 8 with x0(t = 0) = x1(t = 0) = 1.
dx0 =−6×0−3×1 dt
dx1 dt
dx0 =f0(t,x0,x1) dt
dx1 =f1(t,x0,x1) dt
If all fis are linear, then you can write above 2 equations in matrix form
d {x0} [K11 K12] {x0} dt x1 = K21 K22 x1
d {x} = [K]{x} dt
Apply implicit Euler
{xn+1} = {xn} + [K ]{xn+1} Δt [I − [K ]Δt]{xn+1} = {xn}
{xn+1} = [I − [K]Δt]−1{xn}
dx0 =−6×0−3×1 dt
dx1 =5×0+2×1 dt
We can write the above 2 eqns in matrix form
d {x0} = [−6 −3] {x0} dt x1 5 2 x1
For this question
So
[K]=[−6 −3] 52
{xn+1,0} = xn+1,1
[I−[K]Δt]=[1 0]−[−6 −3]Δt 0152
{xn,0} xn,1
[[1 0] − [−6 −3] Δt]−1 0152
Lecture22M01.m
[K] = [−6 −3] 52
[I−[K]Δt]=[1 0]−[−6 −3]Δt 0152
[I−[K]Δt]−1 =[[1 0]−[−6 −3]Δt]−1 0152
{xn+1,0} = [[1 0] − [−6 −3] Δt]−1 {xn,0} xn+1,1 01 52 xn,1
2.5
Δt = 0.2 1.5
1
0.5
0 012345678
t
Implicit Euler
MATLAB ode23()
2
[t,x]=MyEulerLinearSysODEs([tmin tmax],x0,Delta_t);
hold off
plot(t,x(:,2),’bo-‘,’linewidth’,2,’markersize’,10)
hold on
plot(tmat,xmat(:,2),’k-‘,’linewidth’,2)
xlabel(‘t’);
ylabel(‘x_1(t)’);
legend(‘Implicit Euler’,’MATLAB ode23()’)
function [t,x]=MyEulerLinearSysODEs(tspan,x0,Delta_t)
K=[-6 -3;5 2];
temp=inv(eye(2)-K*Delta_t);
x(n+1,:)=(temp*x(n,:)’)’;
close all
clear all
tmin=0.0;
tmax=8.0;
x0=[1 1];
Delta_t=0.2;
[tmat,xmat]=ode23(@f,[0 8],[1 1]);
a=tspan(1);
b=tspan(2);
t=a:Delta_t:b;
x=zeros(length(t),numel(x0));
x(1,:)=x0; %setting initial conditions
for n=1:length(t)-1
end
end
function dxdt=f(t,x)
dxdt=[-6*x(1)-3*x(2) ;5*x(1)+2*x(2)];
end
x1(t)
3
x0(t)
Δt = 1.0
x1(t)
2.5
2
1.5
1
0.5
0
-0.5
-1 012345678
t
3
Δt = 0.5
x0(t)
x1(t)
2.5
2
1.5
1
0.5
0
-0.5
-1 012345678
t
3
2.5
2
1.5
1
0.5
0
-0.5
-1 012345678
t
x0(t)
x1(t)
Δt = 0.2
3
2.5
2
1.5
1
0.5
0
-0.5
-1 012345678
t
Δt = 0.1
x0(t)
x1(t)
Note that the solution computed using implicit Euler’s method is stable for all Δt. This is in contrast with the solution computed using explicit Euler’s method which was unstable for larger values of Δt (see Example L21.5).
5
Explicit Euler Implicit Euler
3
Δt = 1.0
x0(t)
x1(t)
2.5
2
1.5
1
0.5
0
-0.5
-1 012345678
t
1500
1000
500
0
-500
-1000
-1500 012345678
t
Δt = 1.0
x0(t)
x1(t)
Note that the solution computed using implicit Euler’s method is stable for all Δt. This is in contrast with the solution computed using explicit Euler’s method which was unstable for larger values of Δt (see Example L21.5).
x(t)
3
Implicit Euler
Δt = 1.0
x0(t)
x1(t)
2.5
2
1.5
1
0.5
0
-0.5
-1 012345678
t
1500
Explicit Euler
Δt = 1.0
x0(t)
x1(t)
1000
500
0
-500
-1000
-1500 012345678
t
x(t)
d {x0} = [−6 −3] {x0} [K ] = [−6 −3] dtx152x1 52
λ0 = − 1 λ1 = − 3 Δt = 1
λ0Δt = − 1 × 1 = − 1, λ1Δt = − 3 × 1 = − 3
−1 −3 -2 -1
Im t
Im t
−1
Re t −3
Re t
End of Example L22.1
Note that we can easily extend the analysis to get time stepping formula for the Crank-Nicolson method. For a linear system of equations.
d dt
Applying implicit Euler’s formula gives
{x }={x}+[K]{x }Δt n+1 n n+1
[I − [K ]Δt]{xn+1} = {xn} {xn+1} = [I − [K ]Δt]−1{xn}
{x} = [K]{x}
Applying Crank-Nicolson formula gives
{xn+1}={xn}+[K]Δt({xn}+{xn+1}) 2
[I − [K ] Δt ] {xn+1} = [I + [K ] Δt ] {xn} 22
{xn+1} = [I − [K ] Δt ]−1 [I + [K ] Δt ] {xn} 22
L22.2:
Implicit Euler’s method for a system of nonlinear equations
Implicit Euler’s method for a system of equations
So far, we have learnt to solve linear system of ODEs using implicit Euler method. What if you want to use implicit Euler to solve a system of nonlinear M + 1 ODEs?
dx0 dt
dx1 dt
dxM dt
dx2 dt
= f0(t, x0, x1, x2, ……, xM)
= f1(t, x0, x1, x2, ……, xM)
= f2(t, x0, x1, x2, ……, xM) ⋮=⋮
= fM(t, x0, x1, x2, ……, xM)
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) IMPLICIT EULER
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) IMPLICIT EULER
dx0 dt
dx1 dt
= f0(t, x0, x1) = f1(t, x0, x1)
IMPLICIT EULER
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1)Δt
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) IMPLICIT EULER
dx0 dt
dx1 dt
dx2 dt
= f0(t, x0, x1, x2) = f1(t, x0, x1, x2) = f2(t, x0, x1, x2)
IMPLICIT EULER
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt xn+1,2 = xn,2 + f2(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt
dx0 dt
dx1 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1)Δt
= f0(t, x0, x1) = f1(t, x0, x1)
IMPLICIT EULER
If fis are all linear, then you can get an explicit expression for xn+1,i. However, if your fis are nonlinear, then you need to solve a nonlinear
set of equations every time step. I will illustrate this with an example.
Example L22.2: Write a MATLAB program that uses the Implicit Euler method to solve
dx0 dt
dx1 dt
dx2 dt
x0(0) = 0,×1(0) = 1,×2(0) = 1
t ∈ [0,10]
= x1x2
= − x0x2
= − 0.5x0x1
dx0 =f(t,x,x,x) dt 0 0 1 2
dx1 =f(t,x,x,x) dt 1 0 1 2
dx2 =f(t,x,x,x) dt 2 0 1 2
IMPLICIT EULER
xn+1,0 = xn,0 + f0(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt xn+1,1 = xn,1 + f1(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt xn+1,2 = xn,2 + f2(tn+1, xn+1,0, xn+1,1, xn+1,2)Δt
dx0 =x1x2 dt
dx1 =−x0x2 dt
dx2 =−0.5x0x1 dt
IMPLICIT EULER
xn+1,0 = xn,0 + Δtxn+1,1xn+1,2 xn+1,1 = xn,1 − Δtxn+1,0xn+1,2 xn+1,2 = xn,2 − 0.5Δtxn+1,0xn+1,1
xn+1,0 = xn,0 + Δtxn+1,1xn+1,2 xn+1,1 = xn,1 − Δtxn+1,0xn+1,2 xn+1,2 = xn,2 − 0.5Δtxn+1,0xn+1,1
At every time step you need to solve a system of 3 equations for the 3 unknowns, xn+1,0, xn+1,1 and xn+1,2. Note that the values of xn,0, xn,1 and xn,2 are usually known.
We can rearrange the 3 equations above to get
xn+1,0 − xn,0 − Δtxn+1,1xn+1,2 {0}
g= xn+1,1 − xn,1 + Δtxn+1,0xn+1,2 = 0 (22.1)
xn+1,2 − xn,2 + 0.5 Δtxn+1,0 xn+1,1 0
This equation is nonlinear and we need some iterative method to solve it. We already learnt how to solve such nonlinear equations in L7.4 (see example L07.5). I will provide a recap here to jog your memory.
End of first part of Example L22.2
Newton-Raphson method to solve multiple equations
We will now recap Newton-Raphson method to a set of equations
g0(p0,p1,p2,p3,…,pM) = 0
g1(p0,p1,p2,p3,…,pM) = 0
g2(p0,p1,p2,p3,…,pM) = 0
g3(p0,p1,p2,p3,…,pM) = 0 ⋮⋮⋮
gM(p0,p1,p2,p3,…,pM) = 0 gis are nonlinear equations and pis are your unknowns.
g0(p0,p1,p2,p3,…,pM) = 0 g1(p0,p1,p2,p3,…,pM) = 0 g2(p0,p1,p2,p3,…,pM) = 0 g3(p0,p1,p2,p3,…,pM) = 0
⋮⋮⋮ gM(p0,p1,p2,p3,…,pM) = 0
To get the solution to the above set of equations, you need to iterate using the equation below. Each iteration should give you a better guess of pi (pi,new). You need to start with a guess, pi,old.
Note that the equation below is a system of linear algebraic equation and you need to solve in order to obtain pi,new.
∂g0/∂p0 ∂g0/∂p1 ∂g1/∂p0 ∂g1/∂p1 ∂g2/∂p0 ∂g2/∂p1 ∂g3/∂p0 ∂g3/∂p1
… ∂g0/∂pM p0,new −p0,old g0(p0,old,p1,old,….,pM,old)
. . . ∂g1/∂pM p1,new − p1,old g1(p0,old, p1,old, . . . . , pM,old)
… ∂g2/∂pM p2,new−p2,old =− g2(p0,old,p1,old,….,pM,old)
. . . ∂g3/∂pM p3,new − p3,old g3(p0,old, p1,old, . . . . , pM,old) ⋮⋮⋮⋮⋮⋮
∂gM /∂p0 ∂gM /∂p2 . . . ∂gM /∂pM pM,new − pM,old gM(p0,old, p1,old, . . . . , pM,old)
This is called the Jacobian. You know the
value of every element of the Jacobian because it is, evaluated at pi,old
You know this because it is a vector of the functions gi s, evaluated at pi,old
Example L22.3: Write your own MATLAB function to solve the system of equations. g0(p0,p1)=2p0−2p1−e−p0 =0
g1(p0,p1)=−p0+2p1−e−p1 =0
The generic formula to perform iteration is
∂g0 /∂p0 ∂g0 /∂p1 ∂g1/∂p0 ∂g1/∂p1 ∂g2 /∂p0 ∂g2 /∂p1 ∂g3/∂p0 ∂g3/∂p1
… ∂g0/∂pM p0,new −p0,old g0(p0,old,p1,old,….,pM,old)
∂gM /∂p0 ∂gM /∂p2 . . . ∂gM /∂pM pM,new − pM,old For a system consisting of 2 equations, this simplifies to
gM(p0,old, p1,old, . . . . , pM,old) g0(p0,old, p1,old)
[
] { 1,new 1,old} g1(p0, p1) = − p0 + 2p1 − e−p1
∂g0 ∂p1
∂g1 =−1 ∂p0
∂g0/∂p0 ∂g0/∂p1
p0,new − p0,old p−p
{ 1
0,old 1,old }
∂g1/∂p0 ∂g1/∂p1 For this example
So
. . . ∂g1/∂pM p1,new − p1,old g1(p0,old, p1,old, . . . . , pM,old)
… ∂g2/∂pM p2,new−p2,old =− g2(p0,old,p1,old,….,pM,old)
. . . ∂g3/∂pM p3,new − p3,old g3(p0,old, p1,old, . . . . , pM,old) ⋮⋮⋮⋮⋮⋮
=− g0(p0, p1) = 2p0 − 2p1 − e−p0
g(p ,p )
∂g0 =2+e−p0 ∂p0
=−2
= 2 + e−p1
∂g1 ∂p1
g0(p0, p1) = 2p0 − 2p1 − e−p0 g1(p0, p1) = − p0 + 2p1 − e−p1
∂g0 ∂p1
=−2
∂g0 ∂p0
= 2 + e−p0
∂g1 ∂p0
=−1
∂g1 ∂p1
= 2 + e−p1
] { 1,new 1,old}
−1 2 + e−p1,old p1,new − p1,old −p0,old + 2p1,old − e−p1,old
pair of p0 and p1 that satisfies both g0( p0, p1) = 0 g1(p0,p1) = 0
[
p0,new − p0,old p−p
∂g0/∂p0 ∂g0/∂p1
g0(p0,old, p1,old) g(p ,p )
{ 1 0,old 1,old }
So we need to solve the matrix equation above for every iteration for p0,new and p1,new until we get a
=−
[2 + e−p0,old −2 ] {p0,new − p0,old} = − { 2p0,old − 2p1,old − e−p0,old }
∂g1/∂p0 ∂g1/∂p1
temp=jacobian\(-gfunc);
p=temp+p;
gfunc=[2*p(1)-2*p(2)-exp(-p(1)); -p(1)+2*p(2)-exp(-p(2))];
jac=[2+exp(-p(1)) -2;-1 2+exp(-p(2))];
clear all
close all
eps=1.0;
p=ones(2,1); %Guess value for p
while eps>1.0e-14
[gfunc,jacobian]=ExampleFunctionWithJacobian(p);
eps=max(abs(gfunc));
end
options = optimoptions(@fsolve,’Display’,’iter’,’SpecifyObjectiveGradient’,true);
[pmatlab,F,exitflag,output,JAC] = fsolve(@ExampleFunctionWithJacobian,[1;1],options);
%Check to see if same as solution given by MATLAB fsolve()
abs(p-pmatlab)
function [gfunc,jac]=ExampleFunctionWithJacobian(p)
end
[2 + e−p0,old −2 ] {p0,new − p0,old} = − { 2p0,old − 2p1,old − e−p0,old } −1 2 + e−p1,old p1,new − p1,old −p0,old + 2p1,old − e−p1,old
{p0,new} = {temp0} + {p1,old} p1,new temp1 p2,old
{p0,new − p0,old} = − [2 + e−p0,old −2 ]−1 { 2p0,old − p1,old − e−p0,old } = {temp0} p1,new − p1,old −1 2 + e−p1,old −p0,old + 2p1,old − e−p1,old temp1
End of Example L22.3
Going back now to Example L22.2, we need to solve
xn+1,0 − xn,0 − Δtxn+1,1xn+1,2 {0} g(x) = xn+1,1 − xn,1 + Δtxn+1,0xn+1,2 = 0
xn+1,2 − xn,2 + 0.5Δtxn+1,0xn+1,1 Similar to previous examples, and for clarity, we let pi = xn+1,i to get
p0 − xn,0 − Δtp1p2 {0} g(p) = p1 − xn,1 + Δtp0 p2 = 0
p2 − xn,2 + 0.5Δtp0 p1 The Jacobian for this system of equation is given by
∂g0/∂p0 ∂g0/∂p1 ∂g0/∂p2 ∂g1/∂p0 ∂g1/∂p1 ∂g1/∂p2 ∂g2/∂p0 ∂g2/∂p1 ∂g2/∂p2
[J] =
For the system of equation we are required to solve, this becomes
Jmn = ∂gm = ∂pn
1 −p2Δt −p1Δt p2Δt 1 p0Δt
0.5p1Δt 0.5p0Δt 1
The generic formula using Newton -Raphson formula to solve a system of 3 nonlinear algebraic equations is
∂g0 /∂p0 ∂g0 /∂p1 ∂g0 /∂p2 ∂g1/∂p0 ∂g1/∂p1 ∂g1/∂p2 ∂g2/∂p0 ∂g2/∂p1 ∂g2/∂p2
p0,new − p0,old
p1,new − p1,old = − p2,new −p2,old
g0(p0,old, p1,old, p2,old) g1(p0,old, p1,old, p2,old) g2(p0,old,p1,old,p2,old)
Hence, the iterative formula for the Newton Raphson method to solve Eq. (22.1) can be written as
1 −p2,oldΔt p2,oldΔt 1
0.5p1,oldΔt 0.5p0,oldΔt
−p1,oldΔt p0,oldΔt
p0,new − p0,old p1,new − p1,old p2,new − p2,old
= −
−1
p0,old − xn,0 − Δtp1,old p2,old
p1,old − xn,1 + Δtp0,old p2,old p2,old − xn,2 + 0.5Δtp0,old p1,old
p0,old − xn,0 − Δtp1,old p2,old
p0,new − p0,old p1,new − p1,old p2,new − p2,old
1 0.5p1,oldΔt
−p2,oldΔt 1
−p1,oldΔt p0,oldΔt
p1,old − xn,1 + Δtp0,old p2,old p2,old − xn,2 + 0.5Δtp0,old p1,old
1 p2,oldΔt
= −
We need to iterate using equation above at every time step until pi,new converges.
0.5p0,oldΔt
1
Lecture22M02.m
clear all;
close all;
Delta_t=0.1;
t=0:Delta_t:50;
%
%Preallocating Memory
%
x=zeros(length(t),3);
x(1,1)=0.0;x(1,2)=1.0;x(1,3)=1.0;
for n=1:length(t)-1
x(n+1,:)=Findxnplus1(Delta_t,x(n,:));
end
%Check with Matlab solution
[tmat,xmat] = ode23(@Example22p2,[0 50],[0 1 1]);
hold off
plot(t,x(:,1),’k-‘,tmat,xmat(:,1),’ko-‘);
hold on
plot(t,x(:,2),’b-‘,tmat,xmat(:,2),’bo-‘);
legend(‘x_0(t) Implicit Euler’,’x_0(t) MATLAB ode23()’,’x_1(t) Implicit Euler’,’x_1(t) MATLAB
ode23()’);
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop at least once
p=xn; %Guess value for p
while eps>1.0e-13
[gfunc,jacobian]=Example22p2WithJacobian(p,Delta_t,xn);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient than inv(jacobian)*(-gfunc)
p=temp’+p;
end end
function [gfunc,jac]=Example22p2WithJacobian(p, Delta_t,xn)
gfunc=[p(1)-xn(1)-p(2)*p(3)*Delta_t; p(2)-xn(2)+p(1)*p(3)*Delta_t; p(3)-xn(3)+0.5*p(1)*p(2)*Delta_t];
jac=[1 -p(3)*Delta_t -p(2)*Delta_t; p(3)*Delta_t 1 p(1)*Delta_t;0.5*p(2)*Delta_t 0.5*p(1)*Delta_t
1];
end
function dxdt = Example22p2(t,x)
dxdt = [x(2)*x(3); -x(1)*x(3);-0.5*x(1)*x(2)];
end
clear all;
close all;
Delta_t=0.1;
t=0:Delta_t:50;
%
%Preallocating Memory
%
x=zeros(length(t),3);
x(1,1)=0.0;x(1,2)=1.0;x(1,3)=1.0;
for n=1:length(t)-1
x(n+1,:)=Findxnplus1(Delta_t,x(n,:));
end
%Check with Matlab solution
[tmat,xmat] = ode23(@Example22p2,[0 50],[0 1 1]);
hold off
plot(t,x(:,1),’k-.’,tmat,xmat(:,1),’k-‘);
hold on
Use Newton-Raphson method to solve the equation below for x
plot(t,x(:,2),’b-.’,tmat,xmat(:,2),’b-‘); n+1 legend(‘x_0(t) Implicit Euler’,’x_0(t) MATLAB
xn+1,0 − xn,0 − Δtxn+1,1xn+1,2 ode23()’,’x_1(t) Implicit Euler’,’x_1(t) MATLAB
f(x)= x − x + Δtx x = ode23()’); n+1,1 n,1 n+1,0 n+1,2 0
x −x +0.5Δtx x n+1,2 n,2 n+1,0 n+1,1
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’
,18);
{0}
ode23()’);
jac=[1 -p(3)*Delta_t -p(2)*Delta_t; p(3)*Delta_t 1
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’
,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop
at least once
p=xn; %Guess value for p
while eps>1.0e-13
[gfunc,jacobian]=Example22p2WithJacobian(p,Delta_t,x
n);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient
than inv(jacobian)*(-gfunc)
p=temp’+p;
end end
Use Newton-Raphson method to solve the equation below for pi,new function [gfunc,jac]=Example22p2WithJacobian(p,
Del
gfu
xn(
xn(
ta_t,xn)
1 −p Δt −p Δt 2,old 1,old
nc=[p(1)-xn(1)-p(2)*p(3)*D
p2,oldΔt 1 p0,oldΔt 2)+p(1)*p(3)*Delta_t; p(3)
0.5p1,oldΔt 0.5p0,oldΔt 1 3)+0.5*p(1)*p(2)*Delta_t];
p0 − xn,0 − Δt p1 p2 p1 − xn,1 + Δtp0p2 p2 − xn,2 +0.5 Δtp0p1
p0,new − p0,old elta_t; p(2)- -p1,new − p1,old = −
p2,new − p2,old
gfunc=[p(1)-xn(1)-p(2)*p(3)*Delta_t; p(2)-
xn(2)+p(1)*p(3)*Delta_t; p(3)-
xn(3)+0.5*p(1)*p(2)*Delta_t];
Use Newton-Raphson method to solve the equation below for pi,new function dxdt = Example22p2(t,x)
dxdt = [x(2)*x(3); -x(1)*x(3);-0.5*x(1)*x(2)];
end
1 −p2,oldΔt p2,oldΔt 1
0.5p1,oldΔt 0.5p0,oldΔt
−p1,oldΔt p0,oldΔt
1
p0,new − p0,old p1,new − p1,old p2,new − p2,old
= −
p0 −
p1 −
p2 −
xn,0 −
xn,1 + xn,2 +0.5
Δtp1p2 Δtp0 p2 Δtp0 p1
[gfunc,jacobian]=Example22p2WithJacobian(p,Delta_t,x
n);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient
than inv(jacobian)*(-gfunc)
p=temp’+p;
end end
function [gfunc,jac]=Example22p2WithJacobian(p,
Delta_t,xn)
jac=[1 -p(3)*Delta_t -p(2)*Delta_t; p(3)*Delta_t 1
p(1)*Delta_t;0.5*p(2)*Delta_t 0.5*p(1)*Delta_t 1];
end
ode23()’);
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’
,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop
at least once
p=xn; %Guess value for p
while eps>1.0e-13
[gfunc,jacobian]=Example22p2WithJacobian(p,Delta_t,x
n);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient
than inv(jacobian)*(-gfunc)
p=temp’+p;
end
end
Calculate
xn,2 +0.5
function [gfunc,jac]=Example22p2WithJa−1cobian(p,
Δtp1p2 Δtp0 p2 = Δtp0 p1
p0,new − p0,old 1 −p2,oldΔt −p1,oldΔt p0 − Delta_t,xn)
xn,0 − xn,1 +
temp0 temp1
temp
p1,new − p1,old = − p2,oldΔt 1 p0,oldΔt p1 − gfunc=[p(1)-xn(1)-p(2)*p(3)*Delta_t; p(2)-
p2,new − p2,old xn(2)+p(1)*p(3)*Delta_t; p(3)-
0.5p1,oldΔt 0.5p0,oldΔt 1 p2 − xn(3)+0.5*p(1)*p(2)*Delta_t];
2
ode23()’);
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’
,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop
at least once
p=xn; %Guess value for p
while eps>1.0e-13
[gfunc,jacobian]=Example22p2WithJacobian(p,Delta_t,x
n);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient
than inv(jacobian)*(-gfunc)
p=temp’+p;
end end
Calculate
function [gfunc,jac]=Example22p2WithJacobian(p,
Delta_t,xn)
gfunc=[p(1)-xn(1)-p(2)*p(3)*Delta_t; p(2)-
p0,new temp0 p0,old p1,new = temp + p1,old
xn(2)+p(1)*p(3)*Delta_t; p(3)- 1
xn(3)+0.5*p(1)*p(22),n*ewDelta_t]; 2 2,old p temp p
jac=[1 -p(3)*Delta_t -p(2)*Delta_t; p(3)*Delta_t 1
Implicit Euler with Δt = 1.0 Implicit Euler with Δt = 0.1
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
x0 (t) Implicit Euler
x0 (t) MATLAB ode23()
0 5 10 15 20 25 30 35 40 45 50 t
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
x0 (t) Implicit Euler
x0 (t) MATLAB ode23()
0 5 10 15 20 25 30 35 40 45 50 t
x0(t)
x0(t)
End of Example L22.2
L22.3:
Crank-Nicolson method for a system of nonlinear equations
Crank-Nicolson method for a system of equations
So far, we have learnt to solve linear system of ODEs using implicit Euler method. What if you want to use Crank-Nicolson method to solve a system of nonlinear M + 1 ODEs?
dx0 dt
dx1 dt
dxM dt
dx2 dt
= f0(t, x0, x1, x2, ……, xM)
= f1(t, x0, x1, x2, ……, xM)
= f2(t, x0, x1, x2, ……, xM) ⋮=⋮
= fM(t, x0, x1, x2, ……, xM)
dx0 dt
xn+1,0 = xn,0 + Δt (f0(tn, xn,0 + f0(tn+1, xn+1,0)) 2
= f0(t,x0) Crank-Nicolson
dx0 dt
xn+1,0 = xn,0 + Δt (f0(tn, xn,0 + f0(tn+1, xn+1,0)) 2
= f0(t,x0) Crank-Nicolson
dx0 dt
dx1 dt
= f0(t, x0, x1) = f1(t, x0, x1)
Crank-Nicolson
xn+1,0 = xn,0 + Δt (f0(tn, xn,0, xn,1) + f0(tn+1, xn+1,0, xn+1,1)) 2
xn+1,1 = xn,1 + Δt (f1(tn, xn,0, xn,1) + f1(tn+1, xn+1,0, xn+1,1)) 2
dx0 dt
xn+1,0 = xn,0 + f0(tn+1, xn+1,0)Δt
= f0(t,x0) Crank-Nicolson
dx0 dt
dx1 dt
dx2 dt
xn+1,0 = xn,0 + Δt (f0(tn, xn,0, xn,1, xn,2) + f0(tn+1, xn+1,0, xn+1,1, xn+1,2)) 2
xn+1,1 = xn,1 + Δt (f1(tn, xn,0, xn,1, xn,2) + f1(tn+1, xn+1,0, xn+1,1, xn+1,2)) 2
xn+1,2 = xn,2 + Δt (f2(tn, xn,0, xn,1, xn,2) + f2(tn+1, xn+1,0, xn+1,1, xn+1,2)) 2
= f0(t, x0, x1, x2) = f1(t, x0, x1, x2) = f2(t, x0, x1, x2)
Crank-Nicolson
dx0 dt
dx1 dt
xn+1,0 = xn,0 + Δt (f0(tn, xn,0, xn,1) + f0(tn+1, xn+1,0, xn+1,1)) 2
xn+1,1 = xn,1 + Δt (f1(tn, xn,0, xn,1) + f1(tn+1, xn+1,0, xn+1,1)) 2
= f0(t, x0, x1) = f1(t, x0, x1)
Crank-Nicolson
If fis are all linear, then you can get an explicit expression for xn+1,i. However, if your fis are nonlinear, then you need to solve a nonlinear
set of equations every time step. I will illustrate this with an example.
Example L22.4: Write a MATLAB program that uses the Crank-Nicolson method to solve
d2φ dt2
+ φ + φ3 = 0
(22.2)
in the domain t ∈ [0,50] given φ(0) = 1, φ′(0) = 0
Equation (22.2) is a 2nd order ordinary differential equation. We first need to write it as a
pair first order ordinary differential equations because the techniques that have been developed only works for first order ordinary differential equations.
Let
So
Note that the initial conditions become
x0 = φ
x1 = dφ/dt
dx0 dt
= x1
dx1 dt
= − x0 − x3 0
x0(0) = 1 φ′(0) = 0 x1(0) = 0
φ(0) = 1
dx0 dt
dx1 dt
=
= f1(t, x0, x1)
dx0 dt
dx1 dt
= =
x1
f0(t, x0, x1)
− x0 − x3 0
Crank-Nicolson
Crank-Nicolson
xn+1,0 = xn,0 + Δt (xn,1 + xn+1,1) 22
xn+1,0 = xn,0 + Δt (f0(tn, xn,0, xn,1) + f0(tn+1, xn+1,0, xn+1,1))
n+1,1 n,1 Δt(1 n n,0 n,1 1 n+1 n+1,0 n+1,1 ) xn+1,1 =xn,1+ 2 −xn,0−x3 −xn+1,0−x3
x = x + f (t , x , x ) + f (t , x , x ) Δt ( n,0 n+1,0) 2
xn+1,0 = xn,0 + Δt (xn,1 + xn+1,1) 2
xn+1,1 =xn,1+Δt(−xn,0−x3 −xn+1,0−x3 ) 2 n,0 n+1,0
At every time step you need to solve a system of 2 equations for the 2 unknowns, xn+1,0 and xn+1,1. Note that the values of xn,0 and xn,1 are usually known.
x
Using similar notation as in the past, we let pi = xn+1,i to get
()
Δt
xn+1,0−xn,0− 2 xn,1+xn+1,1 = 0
g(x)=
n+1,1
33{} −x +Δt(x +x +x +x ) 0
n,1 2 n,0 n,0 n+1,0 n+1,0
()
g(p)= 0 n,0( n,1 1 )=0
Δt p−x−2×3+p3 {0}
p−x +Δt x +x +p+p 1 n,1 2 n,0 n,0 0 0
()
g(p)= 0 n,0( n,1 1 )=0
Δt p−x−2×3+p3 {0}
p−x +Δt x +x +p+p 1 n,1 2 n,0 n,0 0 0
[(3p +1)Δt/2 1 ]{ }
Δt ( ) p0,old − xn,0 − 2 xn,1 + p1,old
2
p1,new − p1,old
p1,old − xn,1 + Δt (xn,0 + xn,0 + p0,old + p0,old) 2
0,old
3 3
1 −Δt/2
p0,new − p0,old
= −
Lecture22M04.m
clear all; close all;
Delta_t=0.05; t=0:Delta_t:50;
%
%Preallocating Memory %
x=zeros(length(t),2); x(1,1)=1.0;x(1,2)=0.0; for n=1:length(t)-1
x(n+1,:)=Findxnplus1(Delta_t,x(n,:));
end
%Check with Matlab solution
[tmat,xmat] = ode23(@Example22p4,[0 50],[1 0]);
1.5
x0 (t) Crank-Nicolson
x0 (t) MATLAB ode23()
1 x1 (t) Crank-Nicolson x1 (t) MATLAB ode23()
0.5
0
-0.5
-1
xi(t)
hold off
plot(t,x(:,1),’k-.’,tmat,xmat(:,1),’k-‘);
hold on
plot(t,x(:,2),’b-.’,tmat,xmat(:,2),’b-‘);
legend(‘x_0(t) Crank-Nicolson’,’x_0(t) MATLAB ode23()’,’x_1(t) Crank-Nicolson’,’x_1(t) MATLAB ode23()’); xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop at least once
p=xn; %Guess value for p
while eps>1.0e-13 [gfunc,jacobian]=Example22p4WithJacobian(p,Delta_t,xn); eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient than inv(jacobian)*(-gfunc) p=temp’+p;
end end
function [gfunc,jac]=Example22p4WithJacobian(p, Delta_t,xn) gfunc=[p(1)-xn(1)-(Delta_t/2.)*(xn(2)+p(2)); p(2)-xn(2)+(Delta_t/2)*(xn(1)+xn(1).^3+p(1)+p(1).^3)];
-1.5
0 5 10 15 20 25 30 35 40 45 50
t
jac=[1 -Delta_t/2; (3*p(1).^2+1)*Delta_t/2. 1];
end
function dxdt = Example22p4(t,x) dxdt = [x(2); -x(1)-x(1).^3];
end
2 0,old
p 1,new
− p 1,old
0,old n,0
1 −Δt/2 [(3p +1)Δt/2 1 ]{
}
Δt ( n,1 1,old) 3
p0,new−p0,old =−
p −x −2 x +p
3 1,old n,1 Δt ( n,0 n,0 0,old 0,old)
p−x+2x+x+p+p
Solve the following using Crank-Nicolson
dx0 dt
dx1 dt
For 0