B. Tech. Chemical Engineering IV Sem.
(2019-2020)
CHC2910: Computer Programming Lab
Study Materials and Home Assignments
1. Euler’s Method for ODEs
This method is the simplest one used to solve ordinary differential equations (ODEs) with a
given initial value.
Euler’s Method is derived from the Taylor’s series
1
𝑦(𝑥 + ℎ) = 𝑦(𝑥) + ℎ ∗ 𝑦 ′ (𝑥) + 2! ℎ2 ∗ 𝑦 ′′ (𝑥)+ . .. (1)
Neglecting higher order terms we get
𝑦(𝑥 + ℎ) = 𝑦(𝑥) + ℎ ∗ 𝑦 ′ (𝑥)
Here 𝑦 ′ =dy/dx and is evaluated at the previously known values of x & y. Therefore
𝑦𝑖+1 = 𝑦𝑖 + ℎ ∗ 𝑦′(𝑥𝑖 )
h indicates step size
ℎ = 𝑡𝑖+1 − 𝑡𝑖
Now, Consider a differential equation dy/dt = f(t, y) = 2y + t with initial condition y(0)=1
Then using Euler’s method, successive approximation of this equation can be given by:
y(t+h) = y(t) + h * f(t(t), y(t))
𝑦𝑖+1 = 𝑦𝑖 + ℎ ∗ 𝑓(𝑡𝑖 , 𝑦𝑖 )
h indicates step size. Choosing smaller values of h leads to more accurate results but more
computation time.
From above problem, given 𝑦0 = 𝑦(0) = 1 at 𝑡0 = 𝑡(0) = 0
If h = 0.02, then value of y at t=0.02 can be calculated as:
𝑦1 = 𝑦0 + 0.02 ∗ 𝑓(𝑡0 , 𝑦0 )
𝑦1 = 𝑦(0.02) = 1 + 0.02 ∗ (2 ∗ 1 + 0) = 1.04
Note: Note: 𝑦1 or y(0.02) is the value of y at 𝑡1 i.e. t=0.02
Similarly 𝑦2 = 𝑦(0.04) at 𝑡2 i.e t=0.04 can be calculated as
B. Tech. Chemical Engineering IV Sem. (2019-2020)
CHC2910: Computer Programming Lab
𝑦2 = 𝑦1 + 0.02 ∗ 𝑓(𝑡1 , 𝑦1 )
𝑦2 = 𝑦(0.04) = 1.04 + 0.02 ∗ (2 ∗ 1.04 + 0.02) = 1.082
And so on.
Similarly you can calculate the value of y at any time t using a small time step.
2. Runge-Kutta 2nd order method for ODEs
Runge Kutta 2nd order method is a very popular method for solving ODEs and can be derived
from Taylor’s series (S.K. Gupta, 1999). This method is more accurate than Euler’s method
but requires more computation.
For this we consider first 3 terms of Taylor’s series while in Euler’s method we consider only
two terms.
1 2
𝑦(𝑥 + ℎ) = 𝑦(𝑥) + ℎ ∗ 𝑦 ′ (𝑥) + ℎ ∗ 𝑦′′(𝑥)
2!
1
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 )ℎ + 𝑓′(𝑥𝑖 , 𝑦𝑖 )ℎ2
2!
For dy/dx = f(x,y) , y(0) = 𝑦0 at x=0
Runge Kutta method for 2nd order can be written as
1 1
𝑦𝑖+1 = 𝑦𝑖 + ( 𝑘1 + 𝑘2 ) ℎ
2 2
where,
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 ) & 𝑘2 = 𝑓(𝑥𝑖 + ℎ, 𝑦𝑖 + 𝑘1 ℎ)
This form allows one to take advantage of the 2nd order method without having to calculate
f′(x, y).
For the previous problem dy/dt = f(t, y) = 2y + t , given 𝑦0 = 𝑦(0) = 1 at 𝑡0 = 𝑡(0) = 0
If h = 0.02, then value of y at t=0.02 can be calculated as:
1 1
𝑦1 = 𝑦0 + 0.02 ∗ ( 𝑘1 + 𝑘2 )
2 2
𝑘1 = 𝑓(𝑡0 , 𝑦0 ) = 2 ∗ 1 + 0 = 2
𝑘2 = 𝑓(𝑡0 + ℎ, 𝑦0 + 𝑘1 ℎ) = 2 ∗ (1 + 2 ∗ 0.02) + (0 + 0.02) = 2.1
1 1
𝑦1 = 1 + 0.02 ∗ ( ∗ 2 + ∗ 2.1) = 1.041
2 2
B. Tech. Chemical Engineering IV Sem. (2019-2020)
CHC2910: Computer Programming Lab
Similarly you can calculate the value of y at any time t using a small time step.
Now, suppose you have two first order simultaneous equations:
𝑑𝑦
𝑓1 = = 2𝑥 + 2𝑦𝑡
𝑑𝑡
𝑑𝑥
𝑓2 = = 2𝑦 + 2𝑥𝑡
𝑑𝑡
Both equations will have to be solve simultaneously as
1 1
𝑥𝑖+1 = 𝑥𝑖 + ( 𝑘11 + 𝑘12 ) ℎ
2 2
1 1
𝑦𝑖+1 = 𝑦𝑖 + ( 𝑘21 + 𝑘22 ) ℎ
2 2
where,
𝑘11 = 𝑓1 (𝑡𝑖 , 𝑥𝑖 , 𝑦𝑖 )
𝑘12 = 𝑓1 (𝑡𝑖 , 𝑥𝑖 + ℎ, 𝑦𝑖 + 𝑘11 ℎ)
𝑘21 = 𝑓2 (𝑡𝑖 , 𝑥𝑖 , 𝑦𝑖 )
𝑘22 = 𝑓2 (𝑡𝑖 , 𝑥𝑖 + ℎ, 𝑦𝑖 + 𝑘11 ℎ)
B. Tech. Chemical Engineering IV Sem. (2019-2020)
CHC2910: Computer Programming Lab
Assignment #3:
(i) Solution of coupled ODEs
A reactant A decomposes according to the series reaction:
𝑘1 𝑘2
𝐴 → 𝐵 → 𝐶
where k1 = 3s-1 and k2 = 1s-1
The rate expressions follow simple first order kinetics as follows:
𝑑𝐶𝐴
= −𝑘1 𝐶𝐴
𝑑𝑡
𝑑𝐶𝐵
= 𝑘1 𝐶𝐴 − 𝑘2 𝐶𝐵
𝑑𝑡
𝑑𝐶𝐶
= 𝑘2 𝐶𝐵
𝑑𝑡
with CA0 = 1 gmol/lit and CB0 = CC0 = 0.
Write a C program to solve the above set of ordinary differential equations using Runge-
Kutta and Euler methods. Use this program to calculate the concentration of A, B, and C after
t = 10 s. Also sketch the concentration profile of A, B, and C upto 10 s.
