Numerical Methods in Engineering (ENGR20005)
Lecture 21 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L21.1:
Solving ODEs with MATLAB inbuilt functions
• There are many inbuilt MATLAB functions to solve ODEs: Initial Value Problems. Some of them are
ode45() ode23() ode23s() ode34s() ode15s()
’s’ is for stiff equation
• We will now learn how to use these functions to solve ODEs and then learn a bit of theory that provide insight on the need and differences between the different solvers.
Example L21.1: (taken from “Applied numerical methods with MATLAB” -S. Chapra) Use the ode23() function and write a MATLAB program to solve the following ordinary differential equations
dx dt
= 10e”(t”2)2/[2(0.075)2] ” 0.6x
for t ! [0,4] and the initial value is given as x(0) = 0.5.
t ! [0,4]
x(0) = 0.5
Lecture21M01.m
[t,x]=ode23(@f,[0 4],0.5); plot(t,x,’ko-‘) xlabel(‘t’,’fontsize’,16); ylabel(‘x(t)’,’fontsize’,16);
function dxdt=f(t,x) dxdt=10*exp(-(t-2)^2/(2*0.075^2))-0.6*x; end
Output
dx = 10e”(t”2)2/[2(0.075)2] ” 0.6x dt
large #t
small #t
• Note that the step size, #t, is not constant.
• ode23() adjusts the step size, #t, of the numerical solution in
order to keep the estimated error under a certain specified
value.
• The step size are smaller in regions where the gradient is large.
How does MATLAB determine #t?
In MATLAB, you can specify the bound of the estimated error,
E $ max(RelTol % !x!, AbsTol)
RelTol and AbsTol are the Relative Error and Absolute Error which the user can specify. By default the values of
RelTol=10″3
AbsTol=10″6
In some applications, especially if you are dealing with micro/nano scale problems (i.e. if the values of x in your solution are very small, i.e. of the order of 10″5 or 10″6) you will need to use smaller values for RelTol or AbsTol in order to get acceptable solution. To change the value of RelTol (or AbsTol), use the odeset() function. The following MATLAB code carries out simulations with the value of RelTol=10″5.
Setting the relative tolerance of the error to 10″5.
Output
2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0 0.5 1 1.5 2 2.5 3 3.5 4 t
Lecture21M01b.m
options=odeset(‘RelTol’,1.0e-5);
[t,x]=ode23(@f,[0 4],0.5,options); plot(t,x,’ko-‘) xlabel(‘t’,’fontsize’,16); ylabel(‘x(t)’,’fontsize’,16);
function dxdt=f(t,x) dxdt=10*exp(-(t-2)^2/(2*0.075^2))-0.6*x; end
x(t)
RelTol=10″3 RelTol=10″5
Note that the solution obtained using the smaller value for RelTol is more accurate but has smaller time steps, #t. Conversely, the larger the value of RelTol or AbsTol, the larger the time steps.
2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0 0.5 1 1.5 2 2.5 3 3.5 4 t
x(t)
End of Example L21.1
Systems of Ordinary Differential Equations
The functions to solve ODEs in MATLAB can also be used to solve a system of ODEs
where
d {x} = {f(t, {x})} dt
{x} =
x0
x1
x2
x3
f0(t, {x})
f1(t, {x}) {f}= f2(t,{x}) f3(t, {x})
&&
xM fM(t, {x})
Example L21.2: (taken from “Applied numerical methods with MATLAB” -S. Chapra) Use the ode23() function and write a MATLAB program to solve the following set of ordinary differential equations
dx0 =1.2×0″0.6x0x1 dt
dx1 =”0.8×1+0.3x0x1 dt
for 0 $ t $ 20 and the initial value is given as x (0) = 2 and x (0) = 1. . 0 x1(0) = 1
1
In this problem, we need to find how both x0(t) and x1(t) changes with time, t. Note that this problem is “coupled” because dx0/dt and dx1/dt depends on both x0(t) and x1(t).
0 $ t $ 20
Lecture21M02.m
0 $ t $ 20
dxdt=[1.2*x(1)-0.6*x(1)*x(2);
-0.8*x(2)+0.3*x(1)*x(2)];
options=odeset(‘RelTol’,1.0e-3);
[t,x]=ode23(@f,[0 20],[2 1],options);
plot(t,x(:,1),’k–‘,t,x(:,2),’r-‘)
xlabel(‘t’,’fontsize’,16);
ylabel(‘x’,’fontsize’,16);
legend(‘x_0(t)’,’x_1(t)’);
axis([0 20 0 6]);
function dxdt=f(t,x)
end
Output
x0(0) = 2 x1(0) = 1
Data is stored in the following manner
x0(t) x1(t) t=0
dx0 =1.2×0″0.6x0x1 dt
dx1 =”0.8×1+0.3x0x1 dt
6
5
4
3
2
1
0
0 2 4 6 8 10 12 14 16 18 20
x0(t)
x1(t)
21
x1,0 ‘ x0(t1)
t = #t t = 2#t t = 3#t t = 4#t
x2,0 ‘ x0(t2)
x3,0 ‘ x0(t3)
x4,0 ‘ x0(t4)
x1,1 ‘ x1(t1)
x2,1 ‘ x1(t2)
x3,1 ‘ x1(t3)
x4,1 ‘ x1(t4)
t
x
Note that in MATLAB, data for the dependent vector is stored in matrix form as follows
….
…. …..
…….
& ….
x0(t0)
x0(t1)
x0(t2)
x0(t3) &
x0(tNt)
xM(t0)
xM(t1)
xM(t2)
xM(t3) &
xM(tNt)
x1(t0)
x1(t1)
x1(t2)
x1(t3) &
x1(tNt)
x2(t0)
x2(t1)
x2(t2)
x2(t3) &
x2(tNt)
[x] = Timeseries for x0
[x] is a (Nt + 1) % (M + 1) matrix. In the example above, Nt + 1 is the number of time steps and M + 1 is the number of dependent variables in the system
Timeseriesforx1 Timeseriesforx2
Timeseries for xM
x (t ) 0 0
x (t ) 0 1
x (t ) 0 2
x (t ) 0 3
&& x (t )
0 N tt
[x] =
x1(t0) x2(t0)
x1(t1) x2(t1)
x1(t2) x2(t2)
. . . .
. . . . . . . . .
xM(t0) xM(t1) xM(t2)
The MATLAB command x(:,1) will give you this column vector
x1(t3) x2(t3) ……. xM(t3)
&&&& x1(tNt) x2(tNt) . . . . xM(tNt)
x0(t0) x1(t0)
x0(t1) x1(t1)
. . . .
. . . . . . . . .
xM(t0) xM(t1) xM(t2)
x (t ) 2 0
x (t ) 2 1
x (t ) 2 2
x (t ) 2 3
&& x (t )
2 N tt
x0(t2) x1(t2)
[x] =
x0(tNt) x1(tNt) . . . . xM(tNt)
The MATLAB command x(:,3) will give you this column vector
x0(t3) x1(t3) &&&&
……. xM(t3)
x0(t0) x1(t0)
x0(t1) x1(t1)
x2(t0) x2(t1) x2(t2) x2(t3)
. . . .
. . . .
xM(t0) xM(t1) xM(t2) xM(t3)
x0(t2) x1(t2)
. . . . . . . . . . . .
[x] =
x0(tNt) x1(tNt) x2(tNt) . . . . xM(tNt)
The MATLAB command x(2,:) will give you this row vector
x0(t3) x1(t3) &&&&&
L21.2:
Systems of
Ordinary Differential Equations: Explicit EULER
18
In general, a system of M + 1 ODEs can be written as dx0 =f0(t,x0,x1,x2,……,xM)
dx1 =f1(t,x0,x1,x2,……,xM) dt
dx2 =f2(t,x0,x1,x2,……,xM) dt
dt
&=&
dxM =fM(t,x0,x1,x2,……,xM)
dt
Explicit Euler’s method can easily be extended to solve a system of M + 1 equations
dx0 dt
x =x +f(t,x )#t n+1,0 n,0 0 n n,0
= f0(t,x0) EXPLICIT EULER
xn,0
is the approximate value of x0(t) at time level t = tn, i.e. x0(tn).
In general, xn,i
is the approximate value of xi(t) at time level t = tn, i.e. xi(tn).
dx0 dt
x =x +f(t,x )#t n+1,0 n,0 0 n n,0
= f0(t,x0) EXPLICIT EULER
dx0 =f0(t,x0) dt
EXPLICIT EULER
xn+1,0 = xn,0 + f0(tn, xn,0)#t
dx0 =f0(t,x0,x1) dt
dx1 =f1(t,x0,x1) dt
EXPLICIT EULER
xn+1,0 = xn,0 + f0(tn, xn,0, xn,1)#t xn+1,1 = xn,1 + f1(tn, xn,0, xn,1)#t
dx0 =f0(t,x0) dt
EXPLICIT EULER
x0,n+1 = x0,n + f0(tn, x0,n)#t dx0 =f0(t,x0,x1)
dx1 =f1(t,x0,x1) dt
EXPLICIT EULER
xn+1,0 = xn,0 + f0(tn, xn,0, xn,1)#t xn+1,1 = xn,1 + f1(tn, xn,0, xn,1)#t
dx0 =f0(t,x0,x1,x2) dt
dx
1 =f1(t,x0,x1,x2)
dt
dx2 = f2(t, x0, x1, x2) dt
EXPLICIT EULER
xn+1,0 = xn,0 + f0(tn, xn,0, xn,1, xn,2)#t xn+1,1 = xn,1 + f1(tn, xn,0, xn,1, xn,2)#t xn+1,2 = xn,2 + f2(tn, xn,0, xn,1, xn,2)#t
dt
EXAMPLE L21.3:
Write a MATLAB program that uses explicit Euler’s method to solve the following system of two ODEs
dx0 =”6×0″3×1 dt
= 5×0 + 2×1
for 0 $ t $ 8 with x0(t = 0) = x1(t = 0) = 1.
dx1 dt
So we want to calculate
x
x0(t1)
x0(t0) x1(t1)
x0(t2) x1(t2)
x1(t0)
x0(t3) x0(t4) x1(t3) x1(t4)
x1(t5)
x0(t)
x0(t5)
x1(t)
t #t t #t t #t t #t t #t t t 012345
We will organise the data in the following manner (similar to MATLAB, see L21.1)
Columns of matrix indicate variable
x0(t0) x1(t0)
x0(t1) x1(t1)
x0(t2) x1(t2)
x0(t3) x1(t3)
x0(t4) x1(t4)
x0(t5) x1(t5)
Rows of matrix indicate time level
Lecture21M03.m
clear all; close all;
Delta_t=0.2; t=0:Delta_t:8; x=zeros(length(t),2);
for n=1:length(t)-1
end
[tmat,xmat]=ode23(@f,[0 8],[1 1]);
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)’);
See how this loop works
x(1,:)=[1.0 1.0];
x0(0) = x1(0) = 1
k1=f(t(n),x(n,:))’; x(n+1,:)=x(n,:)+Delta_t*k1;
Columns of matrix indicate variable
101 (t) (t)
10
1
function dxdt=f(t,x) dxdt=[-6*x(1)-3*x(2) ;5*x(1)+2*x(2)]; end
dx0 =”6×0″3×1
Rows of matrix indicate time level
dt 04 0
x (t )
dx1 =5×0+2×1 dt
1
2
x (t ) 0 1
x (t )
x (t ) 01 1
x (t ) 01 2
0
02
x (t ) 0 3
x (t )
0
13
x (t ) 0
14
x (t )
x (t ) 0
0
05
15
clear all; close all;
Delta_t=0.2; t=0:Delta_t:8; x=zeros(length(t),2);
x(1,:)=[1.0 1.0];
for n=1:length(t)-1 k1=f(t(n),x(n,:))’;
x(n+1,:)=x(n,:)+Delta_t*k1;
end
[tmat,xmat]=ode23(@f,[0 8],[1 1]);
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)’);
function dxdt=f(t,x) dxdt=[-6*x(1)-3*x(2) ;5*x(1)+2*x(2)]; end
Columns of matrix indicate variable
)
01 (t ) 2
1 (t 10
1
1
Rows of matrix indicate time level
x0(t1) x0(t2) x0(t3) x0(t4) x0(t5)
x1(t1) x1(t2) x1(t3) x1(t4) x1(t5)
clear all; close all;
Delta_t=0.2; t=0:Delta_t:8; x=zeros(length(t),2);
x(1,:)=[1.0 1.0];
for n=1:length(t)-1 k1=f(t(n),x(n,:))’;
x(n+1,:)=x(n,:)+Delta_t*k1;
end
[tmat,xmat]=ode23(@f,[0 8],[1 1]);
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)’);
function dxdt=f(t,x) dxdt=[-6*x(1)-3*x(2) ;5*x(1)+2*x(2)]; end
3
2.5
2
1.5
1
0.5
0 012345678
Explicit Euler
#t = 0.2 MATLAB ode23()
t
x1(t)
End of Example L21.3
L21.3:
Systems of
Ordinary Differential Equations: Runge-Kutta Methods
32
In general, a system of M + 1 ODEs can be written as dx0 =f0(t,x0,x1,x2,……,xM)
dx1 =f1(t,x0,x1,x2,……,xM) dt
dx2 =f2(t,x0,x1,x2,……,xM) dt
&=&
dxM =fM(t,x0,x1,x2,……,xM)
Runge-Kutta method can easily be extended to solve a system of M + 1 equations
dt
dt
dx0 dt
xn+1,0 = xn,0 + #t ( 12 k1 + 12 k2)
k1 = f(tn, xn,0)
k2 = f(tn + #t, xn,0 + #tk1)
= f0(t, x0) RK-2
dx0 =f(t,x) dt 0 0
RK-2
xn+1,0=xn,0+#t(12k1+12k2)
k1 = f (tn, xn,0)
k2 = f (tn + #t, xn,0 + #tk1)
dx0 =f0(t,x0,x1) dt
dx1 =f1(t,x0,x1) dt
RK-2
k1 associated with x0
xn+1,0 = xn,0 +#t
(
12k11 + 12k21
) k2 associated with x0
k2 associated with
xn+1,1=xn,1+#t(12k12+12k22) x1
k =f(t,x ,x ) k1associatedwith
11 0 n n,0 n,1 x1 k12 = f1(tn, xn,0, xn,1)
k21 = f0 (tn + #t, xn,0 + #tk11, xn,1 + #tk12) k22 = f1 (tn + #t, xn,0 + #tk11, xn,1 + #tk12)
dx0 =f0(t,x0) dt
RK-2
xn+1,0=xn,0+#t(12k1+12k2) k1 = f(tn,xn,0)
dx0
dt = f0(t, x0, x1, x2)
d x1
dt =f1(t,x0,x1,x2)
dx2 = f (t,x ,x ,x )
k1 associated with x0
k2 associated with x0
k2 associated with x1
k2 = f (tn + #t, xn,0 + #tk1) d x0
dt xn+1,0=xn,0+#t(1k11+1k21)
2 0 1
2
dt =f0(t,x0,x1) dx1 =f(t,x,x)
RK-2
dt 1 0 1 RK-2
xn+1,1 = xn,1 +#t 12k12 + 12k22 (2 2 )
k1 associated with x
xn+1,0 = xn,0 +#t 2k11 + 2k21 (1 1 )
x =x+#t(1k+1k) n+1,1 n,1 2 12 2 22
k11 = f0(tn, xn,0, xn,1) k12 = f1(tn, xn,0, xn,1)
k21 = f0 (tn + #t, xn,0 + #tk11, xn,1 + #tk12) k22 = f1 (tn + #t, xn,0 + #tk11, xn,1 + #tk12)
(11) k 2 associated with x =x +#t k + k 1
k21 = f0 (tn + k22 = f1 (tn + k23 = f2 (tn +
n+1,2 n,2 2 13 2 23 x2
k1 associated with
k11 = f0(tn, xn,0, xn,1, xn,2) x2 k12 = f1(tn, xn,0, xn,1, xn,2)
k13 = f2(tn, xn,0, xn,1, xn,2)
#t, xn,0 + #tk11, xn,1 + #tk12, xn,2 + #tk13)
#t, xn,0 + #tk11, xn,1 + #tk12, xn,2 + #tk13) #t, xn,0 + #tk11, xn,1 + #tk12, xn,2 + #tk13)
EXAMPLE L21.4:
Write a MATLAB program that uses 2nd order Runge-Kutta method to solve the following system of two ODEs
for 0 $ t $ 8 with x (t = 0) = x (t = 0) = 1. 0 1
dx0 =”6×0″3×1 dt
dx1 =5×0+2×1 dt
Lecture21M04.m
dx0 =f0(t,x0,x1) dt
dx1 =f1(t,x0,x1) dt
clear all; close all;
Delta_t=0.2;
0$t$8 x=zeros(length(t),2);
t=0:Delta_t:8;
x(1,:)=[1.0 1.0];
for n=1:length(t)-1
end
x0(t =0)=x1(t =0)=1
k1=f(t(n),x(n,:))’;
k2=f(t(n)+Delta_t,x(n,:)+Delta_t*k1)’;
x(n+1,:)=x(n,:)+Delta_t*(k1/2.0+k2/2.0);
[tmat,xmat]=ode23(@f,[0 8],[1 1]);
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)’);
function dxdt=f(t,x) dxdt=[-6*x(1)-3*x(2) ;5*x(1)+2*x(2)]; end
k11 = f0(tn, xn,0, xn,1) k12 = f1(tn, xn,0, xn,1)
k21 = f0 (tn + #t, xn,0 + #tk11, xn,1 + #tk12) k22 = f1 (tn + #t, xn,0 + #tk11, xn,1 + #tk12)
xn+1,0 = xn,0 +#t(12k11 + 12k21) xn+1,1 = xn,1 +#t(12k12 + 12k22)
dx0 =”6×0″3×1 dt
dx1 =5×0+2×1 dt
Lecture21M04.m
clear all; close all;
Delta_t=0.2; t=0:Delta_t:8; x=zeros(length(t),2); x(1,:)=[1.0 1.0];
for n=1:length(t)-1
k1=f(t(n),x(n,:))’; k2=f(t(n)+Delta_t,x(n,:)+Delta_t*k1)’; x(n+1,:)=x(n,:)+Delta_t*(k1/2.0+k2/2.0);
end
[tmat,xmat]=ode23(@f,[0 8],[1 1]);
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)’);
function dxdt=f(t,x) dxdt=[-6*x(1)-3*x(2) ;5*x(1)+2*x(2)]; end
2.5
2
1.5
1
0.5
MATLAB ode23()
RK-2
#t = 0.2
0 012345678
t
x1(t)
L21.4:
Stability of System of Equations
40
Stability of a system of ODEs
We will now perform stability analysis for a system of linear M + 1 ODEs.
dx0 =f0(t,x0,x1,x2,……,xM) dt
dx1 =f1(t,x0,x1,x2,……,xM) dt
dx2 =f2(t,x0,x1,x2,……,xM) dt
&=&
dxM =fM(t,x0,x1,x2,……,xM)
If ALL the fis are linear, then we can write the above set of equations in matrix form as
dt
dx0 = f0(t,x0,x1,x2,……,xM)
dt d
dx1 =f1(t,x0,x1,x2,……,xM) dt
dx2 = f2(t,x0,x1,x2,……,xM) dt
&=&
dxM =fM(t,x0,x1,x2,……,xM)
For example
dx0 =”6×0″3×1 =f0(t,x0,x1) dt
dx1 =5×0+2×1=f1(t,x0,x1) dt
So in this example
dt
=
x0 x1
x& M
[K]
{x}
dt
d {x} = dt
{x} = {x0}
d {x0} = [“6 “3] {x0} dtx1 52×1
[K] = [“6 “3] x1 52
K0,0 K0,1 . . K0,M
K1,0 K1,1 . . K1,M
K2,0 K2,1 . . K2,M …..
KM,0 KM,2 . . KM,M
x0 x1
x& M
d {x} = [K]{x} dt
We will assume that [K] is diagonalizable, meaning that [K] can be written as
[K] = [(][)][(]”1
where [(] is the matrix of eigenvectors and [)] is a diagonal matrix of the eigenvalues of [K].
d {x} = [K]{x} [K] = [)][(]”1
dt !0,0 !0,1 . . !0,M
!1,0 !1,1 . . !1,M !2,0 !2,1 . . !2,M
[(] =
…..
!M,0 !M,1 . . !M,M
“0 0 . . 0
0 “1 . . 0
[)] =
0 0 . . “M
!i,j is the eigenvector corresponding to the eigenvalue “i
0 0 “2 . 0 …..
[(]
d {x} = [K]{x} [K] = [(][)][(]”1
dt
d {x} = [K] {x} dt d {x} = [(][)][(]”1{x}
d dt
dt [(]”1{x} = [)][(]”1{x}
If we now let
[(]”1{x} = {#} The equation above can be written as
d {#} = [)]{#} dt
Diagonal matrix
(21.1)
d #1 = “1#1 00dt
d#0 =”0#0 dt
“0 0 . . 0
0 “1 . . 0
0 0 “2 . 0 …..
0 0 . . “M
##
d #1 #1 d#2 =”2#2
==dt
dt #& #& d#3=”3#3
MMdt
&=&
d #M = “M #M dt
d {#} = [)]{#} dt
#0 “00..0 #0 d #1 = 0 “1 . . 0 = #1
dt& 00″2.0 & #M ….. #M
0 0 . . “M
• Since # is a linear combination of xj (see Eq. (21.1)), each of the
#js must be stable in order for the numerical solution for xj to be stable.
• We have already analysed the stability of numerical solution for equation that looks very much like Eq. (21.2).
• From lectures 18, 19 and 20 it was shown that in order for the numerical solution to be stable “j#t must lie within the stability region
of the numerical schemes. Hence, for a system of linear ODEs, ALL “j#t must be within the stability region of the numerical scheme.
d#0 = “0#0 dt
d#1 = “1#1 dt
d#2 = “2#2 dt
d#3 = “3#3 dt
&=&
d#M =”M#M dt
(21.2)
Example L21.5:
For the set of ODEs shown below, use the stability analysis to find the maximum value of #t that you can use with the explicit Euler and RK-2 methods to find numerical solution to the following set of ordinary differential equation
dx0 =”6×0″3×1 dt
dx1 =5×0+2×1 dt
for 0 $ t $ 8 with x0(t = 0) = x1(t = 0) = 1.
We can write the above 2 eqns in matrix form
d {x0} = [“6 “3] {x0} dtx1 52×1
So for this example,
d {x0} = [“6 “3] {x0} dtx1 52×1
[K] = [“6 “3] 52
The eigenvalues and eigenvectors of [K ] are {#0,0} {“3/5}
“0=”1 #1,0 = 1
“1 =”3 {#0,1}={“1}
#1,1 1
For stability, ALL “#t must lie within the stability region.
“0 = ” 1 “1 = ” 3
#t = 2
“0#t = ” 1 % 2 = ” 2,
“1#t = ” 3 % 2 = ” 6 -2
-6 -2
Im
Still well outside the stability region of the explicit Euler method. So the numerical solution will be unstable.
stability region of the explicit Euler method
-1
Re t
2500 t 2000 1500 1000 500 0 -500 -1000 -1500 -2000 -2500
012345678 t
x0(t)
x1(t)
x(t)
“0 = ” 1 “1 = ” 3
#t = 1
“1#t = ” 1 % 1 = ” 1, “2#t = ” 3 % 1 = ” 3
Im
-3
-1
-2
-1
Re t
1500
x0(t)
x1(t)
t 1000 500 0 -500 -1000 -1500
012345678 t
x(t)
“0 = ” 1 “1 = ” 3
#t = 0.5
“1#t = ” 1 % 0.5 = ” 0.5, “2#t = ” 3 % 0.5 = ” 1.5
Im t
-2 -1.5 -1
-0.5
Re t
5 4 3 2 1 0
-1 -2 -3 -4
x (t) 0
x (t) 1
012345678 t
x(t)
“0 = ” 1 “1 = ” 3
#t = 0.1
“1#t = ” 1 % 0.1 = ” 0.1, “2#t = ” 3 % 0.1 = ” 0.3
Im
-2
-1
Re t
t 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5
012345678 t
x0(t)
x1(t)
x(t)
End of Example L21.5
Example 21.6:
Analyse the stability of if you use Euler’s method to solve the 2nd order ordinary differential equations
d2$ + d$ + 4$ = 0 (21.3) dt2 dt
for 0 $ t $ 10 with $(t = 0) = 1 and d$/dt(t = 0) = 0.
If we let
x0 = $
x1 = d$ dt
(21.4) (21.5)
d2$ + d$ + 4$ = 0 dt2 dt
Then Eq. (21.3) becomes
dx1 dt
(21.3)
x0 = $ x1 = d$
(21.4) (21.5)
dt = ” 4×0 ” x1
(21.6)
Eqs. (21.5) and (21.6) can be written in matrix form as
1 ] {x0}
d {x0} = [ 0
dt x1 “4 “1 x1
We should now be able to solve the equation above using explicit Euler and RK-2 with initial conditions x0(0) = 1.0 and x1(0) = 0. The stability of the solution would be dependent on the eigenvalues (multiply by #t) of the matrix [K ] shown below
“0 = ” 0.5000 + 1.9365i “1 = ” 0.5000 ” 1.9365i
[K]
If you use explicit Euler’s method to solve
Eq. (21.3), then
it should be stable for
#t = 0.1
Stability region of explicit Euler’s method.
“0#t #t = 1
“0#t #t = 0.1
“Im#t
2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2
-2.5
-5 -4.5 -4
“1#t #t = 0.1
“1#t #t = 0.5
“1#t #t = 1
“0#t #t = 0.5
-3.5 -3
-2.5 -2
-1.5 -1 -0.5 0
“Re#t
End of Example 21.6
L21.5:
System of Nonlinear Equations
58
Example L21.7: Write a MATLAB program that uses the explicit Euler and RK2 methods to solve
t ! [0,50]
dx0 =x1x2 dt
dx1 =”x0x2 dt
dx2 =”0.5x0x1 dt
x0(0) = 0,×1(0) = 1,×2(0) = 1
Lecture21M07.m
dx0 =f0(t,x0,x1,x2) dt
dx1 =f1(t,x0,x1,x2) dt
dx2 =f2(t,x0,x1,x2) dt
clear all;
close all;
t ! [0,50] xRK2=zeros(length(t),3);
xEuler=zeros(length(t),3);
for n=1:length(t)-1
Delta_t=0.05;
t=0:Delta_t:50;
xRK2(1,:)=[0.0 1.0 1.0];
xEuler(1,:)=[0.0 1.0 1.0];
k1=f(t(n),xRK2(n,:))’;
k2=f(t(n)+Delta_t,xRK2(n,:)+Delta_t*k1)’;
xRK2(n+1,:)=xRK2(n,:)+Delta_t*(k1/2.0+k2/2.0);
xEuler(n+1,:)=xEuler(n,:)
+Delta_t*f(t(n),xEuler(n,:))’;
end
[tmat,xmat]=ode23(@f,[0 50],[0 1 1]);
hold off
plot(t,xEuler(:,2),’b-‘,’linewidth’,2,’markersize’,10)
hold on
plot(t,xRK2(:,2),’r-‘,’linewidth’,2,’markersize’,10)
plot(tmat,xmat(:,2),’k-‘,’linewidth’,2)
xlabel(‘t’);ylabel(‘x_1(t)’);
legend(‘Explicit Euler’,’RK-2′,’MATLAB ode23()’)
x0(0) = 0,×1(0) = 1,×2(0) = 1
function dxdt=f(t,x)
dxdt=[x(2)*x(3); -x(1)*x(3) ;-0.5*x(1)*x(2)];
end
Explicit Euler
xn+1,0 = xn,0 + f0(tn, xn,0, xn,1, xn,2)#t xn+1,1 = xn,1 + f1(tn, xn,0, xn,1, xn,2)#t xn+1,2 = xn,2 + f2(tn, xn,0, xn,1, xn,2)#t
dx0 =x1x2 dt
= f0(t, x0, x1, x2) dx1 =”x0x2 =f1(t,x0,x1,x2)
dx2 =”0.5x0x1 =f2(t,x0,x1,x2) dt
dt
RK-2
k21 = f0 (tn + #t , xn,0 + #t k11, xn,1 + #t k12, xn,2 + #t k13) k22 = f1 (tn + #t , xn,0 + #t k11, xn,1 + #t k12, xn,2 + #t k13) k23 = f2 (tn + #t , xn,0 + #t k11, xn,1 + #t k12, xn,2 + #t k13)
xn+1,0= xn,0+#t(12k11+ 12k21) xn+1,1= xn,1+#t(12k12+ 12k22) xn+1,2= xn,2+#t(12k13+ 12k23)
k11 = f0(tn, xn,0, xn,1, xn,2) k12 = f1(tn, xn,0, xn,1, xn,2) k13 = f2(tn, xn,0, xn,1, xn,2)
[tmat,xmat]=ode23(@f,[0 50],[0 1 1]);
hold off
plot(t,xEuler(:,2),’b-‘,’linewidth’,2,’markersize’,10)
hold on
plot(t,xRK2(:,2),’r-‘,’linewidth’,2,’markersize’,10)
plot(tmat,xmat(:,2),’k-‘,’linewidth’,2)
xlabel(‘t’);ylabel(‘x_1(t)’);
legend(‘Explicit Euler’,’RK-2′,’MATLAB ode23()’)
function dxdt=f(t,x)
dxdt=[x(2)*x(3); -x(1)*x(3) ;-0.5*x(1)*x(2)];
end
3
2
1
0
-1
-2
d x dt
0 = x1x2
dx1 =”x0x2 =f1(t,x0,x1,x2)
dx2 =”0.5x0x1 =f2(t,x0,x1,x2) dt
Explicit Euler
RK-2
MATLAB ode23()
= f0(t, x0, x1, x2) dt
x1(t)
-3
0 5 10 15 20 25 30 35 40 45 50
t
End of Example 21.7