CHC204: Computational Tools for Chemical Engineers Lab
Winter Semester, 2024 - 2025
Instructors: Dr. Priyanka D Pantula & Dr. Bidhan Chandra
Department of Chemical Engineering
Indian Institute of Technology (ISM) Dhanbad
(pantula@[Link] & bidhan@[Link])
Department of Chemical Engineering
Module 1: Solve a system of Linear & Nonlinear equations
Linear Algebraic equations Nonlinear Algebraic equations
Department of Chemical Engineering
Representing System of Linear equations in Matrix form
Department of Chemical Engineering
Reviewing the Rank of a Matrix
Rank - The number of linearly independent columns in a matrix.
Solution Exists -
Department of Chemical Engineering
Rank of a matrix in MATLAB
• Keyword in MATLAB – rank
A = [1 1;2 3]; Matrix representation A = [2 3;4 6]; Matrix representation
b = [4;11]; Column Matrix b = [11;20]; Column Matrix
augA = [A b]; Augmented Matrix augA = [A b]; Augmented Matrix
rank(A) rank(A)
rank(augA) rank(augA)
Department of Chemical Engineering
Homogeneous Case
Department of Chemical Engineering
Review
Department of Chemical Engineering
Solve the following system of Linear equations in MATLAB
Step-wise Procedure
2x+ y + z = 2 1. Declare the variables or unknowns.
-x + y - z = 3 2. Declare the system of equations.
x + 2y + 3z = -10 3. Convert the equations into matrix format Ax = b.
4. Use linsolve to solve Ax = b for the vector of unknowns x.
syms x y z A= b=
[ 2, 1, 1] 2
eqn1 = 2*x + y + z == 2; [-1, 1, -1] 3
eqn2 = -x + y - z == 3; [ 1, 2, 3] -10
eqn3 = x + 2*y + 3*z == -10; X=
3
[A,b] = equationsToMatrix([eqn1, eqn2, eqn3], [x, y, z]) 1
X = linsolve(A,B) -5
Department of Chemical Engineering
Solve the following system of Linear equations in MATLAB
• Method 2 - Use solve instead of linsolve if you have the equations in the form of
expressions and not a matrix of coefficients.
2x+ y + z = 2 • Method 3 - Use \ symbol (mldivide)
-x + y - z = 3
x + 2y + 3z = -10 syms x y z
eqn1 = 2*x + y + z == 2;
syms x y z eqn2 = -x + y - z == 3;
eqn1 = 2*x + y + z == 2; sol = eqn3 = x + 2*y + 3*z == -10;
eqn2 = -x + y - z == 3; x: 3 [A,b] = equationsToMatrix([eqn1,
eqn3 = x + 2*y + 3*z == -10; y: 1 eqn2, eqn3], [x, y, z])
sol = solve([eqn1, eqn2, eqn3], [x, y, z]) z: -5 sol = A\b
• Some of the numerical methods used for solving system of Linear equations –
Gauss Elimination; Gauss Jordan; Gauss Seidel
Department of Chemical Engineering
Solve the following system of Linear equations in MATLAB
2x1+x2 = 1
x1 + 2x2 + x3 = 2 x1 = 3, x2 = -5, x3 = 9
x2 + x3 = 4
x1 – x2 + 3x3 = -8
2x1 – 2x2 + 3x3 = -20 x1 = -23/3, x2 = 13/3, x3 = 4/3
x1 + x2 +x3 = -2
Department of Chemical Engineering
Gauss Elimination Technique
a11x1 + a12x2 + a13x3 = b1 a11 a12 a13 x1 a14
a21x1 + a22x2 + a23x3 = b2 0 a22 a23 x2 = a24
a31x1 + a32x2 + a33x3 = b3 0 0 a33 x3 a34
a11 a12 a13 x1 b1
a21 a22 a23 x2 = b2
a31 a32 a33 x3 b3 a33 ∗ x3 = a34
a34
x3 =
a33
a11 a12 a13 a14
a21 a22 a23 a24 x2 * a22 + x3 * a23 = a24
Back a24 − [x3 ∗ a23]
a31 a32 a33 a34 Substitution x2 =
a22
x1 * a11 + x2 *a12 + x3 * a13 = a14
a11 a12 a13 a14 a14 − [x3 ∗ a13 + x2 ∗ a12]
0 a22 a23 a24 x1 =
a11
0 0 a33 a34
Department of Chemical Engineering
System of Linear Equations – Chemical Engineering
Department of Chemical Engineering
Application of system of linear equations: Balancing Chemical Reactions
Department of Chemical Engineering
Application of system of linear equations: Balancing Chemical Reactions
• 3 equations (homogeneous), 4 unknowns infinitely
many solutions
• Given: x2 = 5. Now unknowns are only three.
On solving Ax = b
3 −1 0 0 0 x1 x1 = 1
A= 8 0 −2 b= 0 = 0 X = x3 x3 = 3
0 −2 −1 −2 ∗ 5 −10 x4 x4 = 4
Department of Chemical Engineering
Application of system of linear equations: Reactors
Ex: A chemical reaction takes place in a series of four continuous stirred tank reactors arranged as shown.
• The chemical reaction
is a first-order
irreversible reaction of
the type
• The conditions of temperature in
each reactor are such that the
value of the rate constant 𝑘𝑖 is
different in each reactor. Also, the
volume of each reactor 𝑉𝑖 is
different.
• The values of 𝑘𝑖 and 𝑉𝑖 are given in the Table
• Assumptions regarding this system:
1. The system is at steady state.
2. The reactions are in the liquid phase.
3. There is no change in volume or density of the liquid.
4. The rate of disappearance of component A in each reactor is given by
−𝑟𝑖 = 𝑉𝑖 𝑘𝑖 𝑐𝐴 𝑖 mol/hr
Department of Chemical Engineering
Application of system of linear equations: Reactors
Solution:
1) Since the system is at steady state, the material balance for each reactor is
Input = output + disappearance by reaction
This balance applied to each of the four reactors yields the following set of equations:
The above is a set of four simultaneous linear algebraic equations with four unknowns.
Department of Chemical Engineering
Application of system of linear equations: Reactors
Current Example:
1) Set up the material balance equation for each of the four reactors. What type of
equations do you have in this set of material balances?
2) Solve the above set of equations in MATLAB and find the exit concentration from
each reactor.
Department of Chemical Engineering
Application of system of linear equations: Reactors
Solution:
2) Solve the four simultaneous linear algebraic equations in MATLAB -
clc
clear all
A = [1100 0 0 0;1000 -1400 100 0;0 1100 -1240 100;0 0 1100
-1250];
b = [1000;0;0;0];
x = A\b
Output (Exit concentrations from each reactor) -
x =
0.9091
0.6969
0.6654
0.5856
Department of Chemical Engineering
Module 1: System of Nonlinear Algebraic Equations
Consider a system expressed as follows -
Example -
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Keyword – fsolve
Syntax – fsolve(fun,x0)
x0 – Initial or Starting point for the algorithm
Requires starting point.
Returns a vector value x.
Implements Levenberg – Marquardt algorithm in the background.
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Example: Solve two nonlinear equations in two variables
Step wise Procedure –
1. Convert the equations to the form .
2. Write a MATLAB script that takes x as input and gives F(x) as output.
3. Using fsolve, solve the nonlinear equations (created in Step 2) by calling giving an initial
or starting point.
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Example: Solve two nonlinear equations in two variables
Solution -
function F = root2d(x)
F(1) = exp(-exp(-(x(1)+x(2)))) - x(2)*(1+x(1)^2);
F(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Solve the system of equations starting at the point [0,0].
fun = @root2d;
x0 = [0,0]; x = 0.3532 and 0.6061
x = fsolve(fun,x0)
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Example : Using 𝑥 (1) = 1.25, obtain the root of F(x) which is given below. Repeat
the calculations using 𝑥 (1) = 1.
1 𝑥
𝐹 𝑥 ≡𝑥− 𝑒 =0
3
function F = root2d(x) x = 1.5121 (1.25 initial point)
F = x – (1/3)*exp(x) x = 0.6191 (1.00 initial point)
“Convergence of solution thus depends on the initial estimate”.
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Example: Using 𝒙(1) = [0.5 0.5], solve F(x) in MATLAB
𝐹1 𝒙 ≡ 4 − 8𝑥1 + 4𝑥2 − 2𝑥13 = 0
𝐹2 𝒙 ≡ 1 − 4𝑥1 + 3𝑥2 + 𝑥22 = 0
The solution is:
x1 = 0.6652 and x2 = 0.4776
Few algorithms for solving a system of nonlinear equations – Successive
Substitution, Newton-Raphson, Quasi-Newton, etc.
Example: Chemical Engineering
Calculate the specific volume of a pure gas, at a given temperature and pressure, by
using the Soave-Redlich-Kwong equation of state -
Department of Chemical Engineering
Example: Chemical Engineering
Department of Chemical Engineering
Example: Chemical Engineering
pressure
(1 atm = 101325 Pa)
Calculate in MATLAB using ‘fsolve’.
Department of Chemical Engineering
Example: Chemical Engineering
First, write the function in the form: F x = 0
Now, create a user-defined function for F and the inputs include all the terms in
the function, including the constants.
function F = nonlinear_rk(V, T, R, P, alpha_constant, a, b)
F = (P * V^3) - (R * T * V^2) + (V * ((-R * T * b) + (a * alpha_constant) - (P
* b^2))) - (a * alpha_constant * b);
Solve F using ‘fsolve’ by calling this user-defined function in a new script file
(Note: Make sure the constants are defined before calling the function).
Department of Chemical Engineering
Example: Chemical Engineering
clc
clear
R = 8314;
TC = 425.2;
PC = 3797 * 10^3;
a = (0.4278 * R^2 * TC^2)/PC;
b = (0.0867 * R * TC)/PC;
acc_fac = 0.1931;
S = (0.48508) + (1.55171 * acc_fac) - (0.15613 * acc_fac^2);
T = 500;
alpha_constant = (1 + (S * (1 - sqrt(T/TC))))^2;
for i = 1:40
P = i*101325;
V0 = (R * T)/P;
fun = @(V)nonlinear_rk(V,T,R,P,alpha_constant,a,b);
x(i) = fsolve(fun,V0)
end
Department of Chemical Engineering
Example: Chemical Engineering
P (atm) V (m^3) P (atm) V (m^3)
1 40.8110923232392 21 1.73031532064666
2 20.2975039148817 22 1.64110961166522
3 13.4593803581419 23 1.55962719446042
4 10.0401216624851 24 1.48490319188554
5 7.98840788335256 25 1.41612719525585
6 6.62046575645629 26 1.35261356719420
7 5.64324964602943 27 1.29377834455646
8 4.91023678236088 28 1.23912109146854
9 4.34002567476643 29 1.18821050767103
10 3.88377547722951 30 1.14067291598296
11 3.51040570050605 31 1.09618297880836
12 3.19919606862790 32 1.05445615612753
13 2.93580182670618 33 1.01524253560739
14 2.70997677442401 34 0.978321752370961
15 2.51420708660363 35 0.943498780521284
16 2.34285745154380 36 0.910600426930179
17 2.19161857279320 37 0.879472394436782
18 2.05713883751100 38 0.849976809555899
19 1.93677230061047 39 0.821990131295789
20 1.82840227420672 40 0.795401374342716
Department of Chemical Engineering
Example: Chemical Engineering
Consider a CSTR operated isothermally with
negligible volume change due to reaction, in
Outlet overflow mode with constant fluid volume V,
and with the following two chemical reactions
(assumed elementary). Calculate the outlet
Reactor concentrations of A, B, C and D at steady state.
inlet
Department of Chemical Engineering
Example: Chemical Engineering
(All in SI units)
Department of Chemical Engineering
Example: Chemical Engineering
Let,
At steady state,
4 Nonlinear equations,
4 unknowns.
Department of Chemical Engineering
Example: Chemical Engineering
function F = cstr(x)
F(1) = (1 - x(1)) - 100 * x(1) * (x(2));
F(2) = (2 - x(2)) + 100 * (-(x(1) * x(2)) - (x(2) * x(3)));
F(3) = (0 - x(3)) + 100 * ((x(1) * x(2)) - (x(2) * x(3)));
F(4) = (0 - x(4)) + 100 * (x(2) * x(3));
clc
clear
fun = @cstr;
x0 = [0,0,0,0];
x = fsolve(fun, x0)
Department of Chemical Engineering
Example: Chemical Engineering
(All in SI units)
Department of Chemical Engineering
Solve a system of nonlinear algebraic equations in MATLAB
Single Variable Newton-Raphson Technique
Example 3: Using 𝒙(1) = 1.25, solve F(x) using the Newton-Raphson technique -
1 𝑥
𝐹 𝑥 =𝑥− 𝑒 =0
3
The algorithm is
(𝒌+𝟏) (𝒌) 𝑭(𝒙)
𝒙 =𝒙 −
𝑭ˈ(𝒙)
|𝐅(𝐱 (𝐤+𝟏) ) − 𝐅 𝐱 𝐤 | < 𝐭𝐨𝐥𝐞𝐫𝐚𝐧𝐜𝐞
Department of Chemical Engineering
Single Variable Newton-Raphson Technique
function functionvalue = example(x)
Step 1: Write a function for
functionvalue = x - ((1.0/3.0) * exp(x));
calculating the function value.
end
Step 2: Write a function for
calculating the derivative of the
given function.
function firstDerivative = examplederivative(x)
firstDerivative = 1 - ((1.0/3.0) * exp(x));
end
Department of Chemical Engineering
Single Variable Newton-Raphson Technique
clc
clear
tolerance = 1.0e-6;
iterNum = 0;
error = 0;
x = input("Enter the initial value");
while((error > tolerance) | (iterNum <= 1)) Step 3: Write a subprogram for
iterNum = iterNum + 1; the newton Raphson method by
slope = examplederivative(x); taking initial guess of 1.25
xPrevious = x; (tolerance = 1e-6)
x = x - (example(x) / slope);
error = abs(x - xPrevious);
disp([iterNum, x, error]);
xPrevious = x;
end
disp(["Solution = ", x])
Department of Chemical Engineering
Output - Enter an initial guess for x: 1.25
Iteration Root Error
1 1.7795 1.7795
2 1.5785 0.2011
3 1.5178 0.0607
4 1.5122 0.0056
5 1.5121 0.0000
6 1.5121 0.0000
Root = 1.51213455
Rapid Convergence, 6 iterations
Department of Chemical Engineering
Example: Chemical Engineering
Calculate the specific volume of a pure gas, at a given temperature and pressure, by
using the Soave-Redlich-Kwong equation of state -
Department of Chemical Engineering
Example: Chemical Engineering
Department of Chemical Engineering
Example: Chemical Engineering
pressure
(1 atm = 101325 Pa)
Calculate in MATLAB using Newton Raphson method.
Calculate the friction factor (f) from the Colebrook equation for flow of a fluid in a pipe where
𝜀, D and NRe are roughness, inside diameter of the pipe and Reynolds number respectively.
𝜀
Given: = 10−4 , NRe = 105 and Initial guess for f = 0.1.
𝐷
Department of Chemical Engineering