Numerical Solution of Differential Equations
1.2 Runge-Kutta Method
dy
To find numerical solution to the initial value problem f ( x, y ), y (0) y 0 using
dx
Runge-Kutta method we have the following consideration (This method gives more
accurate result compared to Euler’s method):
k1 f ( x n , y n )
k 2 f ( x n h / 2, y n hk1 / 2)
k 3 f ( x n h / 2, y n hk 2 / 2)
k 4 f ( x n 1 , y n hk 3 )
h
y n 1 y n ( k1 2 k 2 2 k 3 k 4 )
6
Exercise 4: Simple first order ODE
Consider the following initial value problem:
dx
t2 t
dt
With the initial condition: x (0) = 0.5 write a code to solve this equation using Runge-
Kutta integration method.
t_in=0; t_fin=2; nsteps=10; dt=(t_fin-t_in)/nsteps;
t(1)=t_in; x(1)= 0.5;
f=inline('t^2+t') ;
for n=1:nsteps
t(n+1)=t(n)+dt ;
k1 = f(t(n));
k2 = f( t(n)+dt/2);
k3 = f( t(n) + dt/2);
k4 = f( t(n) + dt);
x(n+1)=x(n)+(1/6)*dt*(k1 +2*k2 + 2*k3 +k4);
end
x_Exact=t.^3/3+t.^2/2+0.5;
plot (t,x,'k*:',t,x_Exact,'k+-')
xlabel('Time (t)'),ylabel('x(t)')
legend('Runge Kutta integration','analytical integration')
The results of this program are plotted in Figure (4), the error between analytical and
Runge- Kutta integration method is too small.
Computer Programming 125
Second Class \ Lec. 13
Runge Kutta integration
5
analytical integration
4
x(t)
3
0 0.5 1 1.5 2
Time (t)
dx
Figure (4): analytical and Runge-Kutta integration for t 2 t function
dt
Comparison between Euler method and Runge-Kutta Method
Exercise 5:
Solve the following equation using both Eular and Runge-Kutta method and to
dy
approximate the solution of the initial value problem x y, y (0) 1 with step size h =
dx
0.1.
Solution:
Eular Runge Kutta
clear all, clc,format short clear all, clc,format short
x(1)=0; x(1)=0;
y(1)=1; y(1)=1;
h=0.1 h=0.1;
for i=1:5 f=inline('x+y);
x(i+1)=x(i)+h; for i=1:5
dy=x(i)+y(i); x(i+1)=x(i)+h;
y(i+1)=y(1)+h*dy; k1 = f(x(i),y(i));
end k2 = f(x(i)+h/2,y(i)+k1*h/2);
y_exact= -1-x+2*exp(x); k3 = f( x(i) + h/2,y(i)+k2*h/2);
error=y_exact-y k4 = f( x(i) + h,y(i)+k3*h);
table=[x',y',y_exact',error'] y(i+1)=y(i)+(1/6)*h*(k1 +2*k2 + 2*k3 +k4);
end
y_exact= -1-x+2*exp(x);
error=y_exact-y
table=[x',y',y_exact',error']
Computer Programming 126
Second Class \ Lec. 13
error = error =1.0e-05 *
0 0.0103 0.1228 0.2677 0 0.0169 0.0375 0.0621
0.4404 0.6431 0.0915 0.1264
table = table =
0 1.0000 1.0000 0 0 1.0000 1.0000 0
0.1000 1.1000 1.1103 0.0103 0.1000 1.1103 1.1103 0.0000
0.2000 1.1200 1.2428 0.1228 0.2000 1.2428 1.2428 0.0000
0.3000 1.1320 1.3997 0.2677 0.3000 1.3997 1.3997 0.0000
0.4000 1.1432 1.5836 0.4404 0.4000 1.5836 1.5836 0.0000
0.5000 1.1543 1.7974 0.6431 0.5000 1.7974 1.7974 0.0000
Questions: discus the above results and compare between it.
1.3 MATLAB Built-In Routines for solving ODES
MATLAB have several sub programs (Routines) used to solve initial value
problems for ordinary differential equations (ODEs):
ode113:Variable order solution to non-stiff system
ode15s:Variable order, multistep method for solution of stiff system
ode23: Lower order adaptive step size routine for non-stiff systems
ode23s:Lower order adaptive step size routine for stiff systems
ode45: Higher order adaptive step size routine for non-stiff system
ode45: ode45 is a Differential Equation Solver used to solve non-stiff differential
equations — medium order method. This routine uses a variable step Runge-Kutta
Method to solve differential equations numerically.
Syntax
[t,y] = ode45(odefun,tspan,y0)
[t,y] = ode45(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)
Arguments
odefun : A function that evaluates the right-hand side of the differential equations.
tspan : A vector specifying the interval of integration, [t0,tf].
y0 : A vector of initial conditions.
options : Optional integration argument created using the odeset function.
Computer Programming 127
Second Class \ Lec. 13
1.3.1 Solving ODES using 0de45
There are three types of solution to solve ODES using 0de45 :
1-using local function: Write a local function in a separate script file which have
ODES and call this local function by name in main program.
2-using local function and called as Anonymous function: call local function as
Anonymous function by name in the main program.
3-using Anonymous function only for ODES with ode45 in one main program.
1- Solve ODES using 0de45 with local function
Exercise 6:
Solve the following initial value problem:
dx
t2 t
dt
With the initial condition: x (0) = 0.5 write a code to solve this equation using ode45
method.
Solution:
To solve the same example using MATLAB Routine, We just have to write a local
function which returns the rate of change of the vector x.
function dx = example1(t,x)
dx(1,1) = t^2+t;
The above file named example.m must be saved in default MATLAB folder then the
following code must run.
t_in=0;t_fin=2; nsteps=10;
dt=(t_fin-t_in)/nsteps;
Vspan = t_in:dt:t_fin;
x0=0.5;
[t,x] = ode45('example1',Vspan,x0);
plot (t,x,'k*:')
xlabel('Time (t)'),ylabel('x(t)')
legend('Runge Kutta integration')
The numerical solution, computed using ODE45 is given below in Figure (5)
Computer Programming 128
Second Class \ Lec. 13
5 Runge Kutta integration
x(t)
3
0 0.5 1 1.5 2
Time (t)
Figure (5): Solution produced by ODE45.
2- Solve ODES using 0de45 with local function called as Anonymous
function
Exercise 7:
Solve the following initial value problem:
dx
t2 t
dt
With the initial condition: x (0) = 0.5 write a matlab code to solve this equation using
ode45 method with Anonymous function.
Solution: you have the above local function called example, called as Anonymous
function in the main program.
t_in=0;t_fin=2; nsteps=10;
dt=(t_fin-t_in)/nsteps;
Vspan = t_in:dt:t_fin;
x0=0.5;
[t,x] = ode45(@example,Vspan,x0);
plot (t,x,'k*:')
xlabel('Time (t)'),ylabel('x(t)')
legend('Runge Kutta integration')
Computer Programming 129
Second Class \ Lec. 13
Figure (6): Solution produced by ODE45
With local function called as Anonymous function
Questions: what is the difference in matlab code between Exercise 5 and Exercise 6?
Also write if there is a difference in plotting figures Figure (5) and Figure (6)
3- Solve ODES using 0de45 using Anonymous function only in one
main program
Exercise 8:
Solve the following initial value problem:
dx
t2 t
dt
With the initial condition: x (0) = 0.5
Write a matlab code to solve this equation using ode45 method with Anonymous
function only in one main program.
Solution:
t_in=0;t_fin=2; nsteps=10;
dt=(t_fin-t_in)/nsteps;
Vspan = t_in:dt:t_fin;
x0=0.5;
[t,x] = ode45(@(t,y) t^2+t,Vspan,x0);
plot (t,x,'k*:')
xlabel('Time (t)'),ylabel('x(t)')
legend('Runge Kutta integration')
Questions: write your observations about the matlab code use in above Exercise 7 and the
difference in matlab code from Exercise 5 and Exercise 6.
Computer Programming 130
Second Class \ Lec. 13
1.3.2 Solve multiple ODE using ode45
Exercise 9:
Let’s consider a simple example of a model of a plug flow reactor that is described by
a system of ordinary differential equations. A plug flow reactor is operated as shown in
Figure (6) below.
A
k1
B
k2
C
Figure (7) Isothermal plug flow reactor
The plug flow initially has only reactant A, the components A react to form
component B. The mole balance for each component is given by the following differential
equations
dC A
u k1C A
dz
dC
u B k1C A k2CB
dz
dC
u C k 2C B
dz
With the following initial values
CA (z=0) =1 kmol/m3 CB (z=0) =0 CC (z=0) =0 and k1=2 k2=3
If u=0.5 m/s and reactor length z=3 m. Solve the multiple differential equations and
plot the concentration of each species along the reactor length by these three methods:
1- Using ode45 with local function.
2-Using ode45 with calling local function as Anonymous function
3-Using ode45 with Anonymous function only in main program.
Computer Programming 131
Second Class \ Lec. 13
Solution:
1- Solve multiple ODE using ode45 with local function
We’ll start by writing the function defining the right hand side (RHS) of the ODEs.
The following function file ‘Exercise8’ is used to set up the ode solver.
function dC= Exercise8 ( z, C)
u = 0.5;k1=2; k2=3;
dC(1,1) = -k1 *C(1) / u;
dC(2,1) = (k1 *C(1)-k2 *C(2)) / u;
dC(3,1) = k2 *C(2)/ u;
Now we’ll write a main script file to call ode45. CA, CB and CC must be defined
within the same matrix, and so by calling CA as C(1), CB as C(2) and CC as c(3), they are
listed as common to matrix C.
The following main file is created to obtain the solution:
clear all, clc
C0 = [1 0 0];
txspan = [0 3];
[z , C] = ode45('Exercise8', txspan, C0)
plot (z,C(:,1),'k+-',z,C(:,2),'k*:',z,C(:,3),'kd-.')
xlabel ('Length (m)');
ylabel ('Concentrations (kmol/m^3) ');
legend ('A', 'B', 'C');
2- Solve multiple ODE using ode45 with calling local function as
Anonymous function
Solution: you have the above local function ‘Exercise8’ called as anonymous
function in the main program.
clear all, clc
C0 = [1 0 0];
txspan = [0 3];
[z , C] = ode45(@(z,C) Exercise8 ( z, C) ,
txspan, C0)
plot (z,C(:,1),'k+-',z,C(:,2),'k*:',z,C(:,3),'kd-.')
xlabel ('Length (m)');
ylabel ('Concentrations (kmol/m^3) ');
legend ('A', 'B', 'C');
Computer Programming 132
Second Class \ Lec. 13
The produced plot is as in Figure (7)
1
A
B
Concentrations (kmol/m3) 0.8
C
0.6
0.4
0.2
0
0 1 2 3
Length (m)
Figure (7): A, B and C concentrations along plug flow reactor
3- Solve multiple ODES using 0de45 with Anonymous function only
First let’s write differential equations:
dC(1,1) = -k1 *C(1) / u;
dC(2,1) = (k1 *C(1)-k2 *C(2)) / u;
dC(3,1) = k2 *C(2)/ u;
As
dC(1,1) = 2*y(1)/u);
dC(2,1) = (2 *y(1)-3 *y(2))/u;
dC(3,1) = 3*y(2)/u
Solution:
clear all, clc
u = 0.5;k1=2; k2=3;
C0 = [1 0 0];
txspan = [0 3];
[z,y]=ode45(@(z,y) [(2*y(1)/u);(2 *y(1)-3
*y(2))/u;3*y(2)/u],txspan , C0)
plot (z,y(:,1),'k+-',z,y(:,2),'k*:',z,y(:,3),'kd-.')
xlabel ('Length (m)');
ylabel ('Concentrations (kmol/m^3) ');
legend ('A', 'B', 'C');
Computer Programming 133
Second Class \ Lec. 13
Practice Problems
Q 1)
The concentration of salt x in a homemade soap maker is given as a function of
time by
dx
37.5 3.5 x
dt
At the initial time, t 0 , the salt concentration in the tank is 50 g/L.
Using Runge-Kutta 4th order method and a step size of, h 1.5 min , what is the
salt concentration after 3 minutes (compare the numerical integration result with
exact solution result)?
Q 2)
Solve the following initial value problem:
dy/dt =2 t
With the initial condition: y0 = 0 , in the time interval [0, 5],
Write a code to solve the above ordinary differential equation using ode45
method and plot the result with these three methods:
1 -Using ode45 with local function .
2-Using ode45 with calling local function as Anonymous function
3-Using ode45 with Anonymous function only in main program.
Computer Programming 134
Second Class \ Lec. 13
