MATLAB Modeling in Engineering Problems

Dimitrios Meimaroglou, UFAZ PSE course in the framework of P.C.C.E. Msc

General Introduction

During this course, a general introduction to MATLABTM will be provided, followed by a series
of targeted applications that will allow comprehending the basic usage principles and functions
for chemical/process engineering applications. During the introductory session, the following
principles will be reviewed:

• Main components of MATLAB, use of the toolbox, access to help.

• Working in the command window:

Variable manipulations.

Basic operators, logical operators and relational operators.

Elementary mathematical functions ("elfun").

Linear Algebra: working with vectors/matrices, creation, basic operations, extraction

of part of a vector/matrix, concatenation, useful functions ("matfun","elmat").

• Working with .m files:

Scripts/user-functions: creation, execution, rules.

PSE - MATLAB Modeling Page 1

Problem 1: Linear Regression

Generalities

The following expression describes a linear relationship between variables y and x:

y = â + b̂ · x (1)

â and b̂ being the estimators of the coefficients, which can be readily calculated as:

n · ni=1 (xi · yi ) − ( ni=1 xi ) · ( ni=1 yi )

P P P
b̂ = (2)
n · ni=1 x2i − ( ni=1 xi )2
P
P

and

â = y − b̂ · x (3)

In the above expressions, n denotes the number of measurements and z is the average (i.e.,
mean) value of z.

Work to be done

1. Read the external text file "Regression.dat" and load the data to appropriate variable(s).

2. Calculate the estimators of the coefficients â and b̂, using loop structures.

3. Plot the data (using distinct points) and the resulting regression line (using a continuous
curve) on the same figure. Add the appropriate axis titles and legend.

4. Repeat without using any loop.

5. Explore the possibilities of specifying the size, type and color of the markers & curves.
Use filled red squares for the data and a blue dash-dotted line for the regression. Specify
appropriate size and width for the markers and line. Use a font size equal to 12 and bold
format for the axis titles and legend. Save the figure as a .png file.

Useful commands: clc; clear vars; load; mean; sum; length; size; plot; xlabel; ylabel; legend.
Expected results: â = −9.2381; b̂ = 7.7766

PSE - MATLAB Modeling Page 2

Problem 2: Numerical Integration

Generalities

For the numerical integration of:

Z b
I= f (x)dx (4)
a

the following numerical methods can be considered:

Trapezoidal Rule Let us consider that the domain is subdivided into a number of intervals,
as follows:
b−a
xi = a + i · h; i = 0, 1, ..., N ; h= . (5)
N
The trapezoidal rule is then applied to each interval as:
h
Ii ≃ · (f (xi ) + f (xi+1 ) (6)
2
Accordingly, the implementation over the whole domain will result to:
N −1  
X f (a) f (b)
I≃ Ii = h · + f (a + h) + ... + f (b − h) +
2 2
i=0
(7)
N −1   !
b−a X b−a
= · f (a) + f (b) + 2 · f a+i
2N N
i=1

Simpson Rule N is chosen to be pair, a interval is defined as by: [x2i , x2i+1 , x2i+2 ], i =
0, 1, ..., N/2 − 1. Application of Simpson’s rule on every interval will give:
h b−a
Ii ≃ · (f (x2i ) + 4 · f (x2i+1 + f (x2i+2 ) ; h= (8)
3 N
Accordingly, the implementation over the whole domain will result to:
N/2−1
X h
I≃ Ii = (f (a) + 4 · f (a + h) + 2 · f (a + 2h) + 4 · f (a + 3h)... + 2 · f (b − 2h) + 4 · f (b − h) + f (b))
3
i=0
 
N/2−1 N/2−1
h X X
= f (a) + f (b) + 2 · f (a + 2 · i · h) + 4 · f (a + (2 · i + 1) · h)
3
i=1 i=1

(9)

Quadrature Gauss-Legendre Using orthogonal polynomials of order n, it is possible to

approximate the value of the integral by:
n
b−aX
I≃ wi f (zi ) (10)
2
i=0

where wi are the weights of the polynomials Legendre (i.e., given as function of the order, n).
Finally, the values of xi correspond to specific roots, zi of the polynomials (also defined according
to the order n). The relation allowing to switch from xi to zi and vice versa is:
2zi − (a + b)
xi = (11)
b−a

PSE - MATLAB Modeling Page 3

Work to be done
R3
1. Apply the above techniques for the calculation of the integral I = −2 x
4 dx.

For trapezoidal and Simpson rules, apply both interval-wise and complete domain
expressions.

For Gauss-Legendre*, try orders between n = 1 and n = 7. Comment on the result.

Compare the results with the analytical solution and comment.

2. Repeat the calculation of the integral using the internal numerical integration function of
MATLAB, "quad", and compare with the previously obtained results.
R3
3. Repeat for the calculation of the integral I = 0 x · exp(x2 )dx.

Useful commands: load, linspace, quad. The use of user-functions will be necessary, at least
for the definition of the expression of f (x).
*The roots and weights of Gauss-Legendre can be found here.
R3
Expected results for I = −2 x4 dx: Itrap,N =50 = 55.1167; ISimp,N =50 = 54.9957; IGL,n=2 =
37.6389; IGL,n=3 = 55.0000.
R3
Expected results for I = 0 x · exp(x2 )dx: Itrap,N =50 = 4097.1029; ISimp,N =50 = 4051.5383;
IGL,n=2 = 959.3304; IGL,n=7 = 4048.9208.

PSE - MATLAB Modeling Page 4

Problem 3: Iterative Solution of a System of Linear Algebraic Equa-
tions

Generalities

A system of linear algebraic equations can be formulated as:

Ax = b (12)
Among the possible numerical methods that can be implemented to solve this system, the meth-
ods of "Jacobi" and "Gauss-Seidel" propose an iterative, easy-to-implement, approach.

Jacobi The sequence of values for xi , i denoting the iteration number, is calculated accord-
ing to the following expression (subject to an initial value x0 ):

(i+1) (i)
X
ajj xj + ajk xk = bj , j = 1, . . . , n (13)
k̸=j
Gauss-Seidel Similarly for xi :
(i+1) (i+1) (i)
X X
ajk xk + ajj xj + ajk xk = bj , j = 1, . . . , n (14)
k<j k>j

Matricial Solution Should the calculation of matricial calculations be easily accessible, (12)
is subject to the following iterative solution:

xi+1 = I − B −1 A xi + B −1 b (15)

where B is the diagonal matrix of A, all other elements being equal to 0, for the method of
Jecobi. For Gauss-Seidel, B is the lower triangular matrix of A, all other elements set to 0.

Work to be done

1. Apply the above techniques to the following system (x0 = [0; 0; 0]):

 7 x1 + 2 x2 + x3 = 21


3 x1 − 8 x2 − x3 = 15 (16)


 x − 2 x + 6 x = 17
1 2 3

2. Use 15 iterations and compare the convergence speed of the two methods.

3. Test the same example after changing the coefficient a22 from -8 to -1. Comment.

4. Repeat the calculation using the internal function of MATLAB "linsolve". Investigate the
different proposed options and compare the results.

Useful commands: inv, linsolve.

Expected solution: x = [2.4764; −1.1990; 2.0209].

PSE - MATLAB Modeling Page 5

Problem 4: Iterative Solution of a System of Non-Linear Algebraic
Equations

Generalities

A system of non-linear algebraic equations can be formulated as:

 
f1 (x1 , . . . , xn )
..
(17)
 
f (x) = 
 . =0

fn (x1 , . . . , xn )
Such systems are commonly solved either by the iterative method of "Newton-Raphson", or by
employing optimization techniques. The first approach will be tested here.

Newton-Raphson Starting again from an initial point, x0 , the sequence of values for xi can
be calculated by:

x(i+1) = x(i) − (Df (x(i) ))−1 f (x(i) ) (18)

where Df denotes the Jacobean matrix (i.e., matrix of first derivatives) of f .

Work to be done

1. Apply the method of Newton-Raphson to the following system:

ln(x21 + x22 ) + 0.5 sin(x1 x3 ) = 0.45

q
x1 (x21 + x23 ) = 0.5 (19)
1
x3 exp(x1 x2 ) = 0.3
starting from: (0.1, −0.5, 1). The forms of f and Df will be calculated in the same user-
function (Df will be calculated analytically).

2. At a first stage, carry out 5 iterations and stop.

3. Propose a pertinent convergence criterion to stop automatically the iterations.

4. Repeat the calculation using the internal function of MATLAB "fsolve". Investigate the
different proposed options and compare the results.

Useful commands: fsolve.

Expected solution: x = [0.1829; −1.0960; 2.7279].

PSE - MATLAB Modeling Page 6

Problem 5: Differential Equations

Generalities

A system of Ordinary Differential Equations (ODEs) is commonly formulated as:

ẋ = f (x, t) (20)

where x represents the vector of the state variables. Several methods exist for solving such
systems numerically. MATLAB contains a vast toolbox of such methods that can be readily
employed via the use of internal functions. These functions are named "odexx" or "odexxs",
""xx" standing for different numbers (eg., 45, 15, etc.).

Work to be done

Case of a single Continuous Stirred Tank Reactor (CSTR) in dynamic state The temporal
evolution of the molar concentration of species A, which is supposed to participate as reactant
in a chemical reaction of order 1, with respect to A, in a CSTR, is described by the following
expression:
dCA 1
= (Ca,in − CA ) − kCA
dt τ
where Ca,in and k denote the molar concentration of A in the inlet stream of the reactor and the
kinetic rate constant of the reaction, respectively. k is considered to be constant throughout the
reaction (i.e., isothermal conditions).

1. Integrate the above ODE from 0 to 10 min, using: Ca,in = 0, 4 mol/l; τ = 2 min and k = 0, 5
min−1 . Test the cases of C(t = 0) = 0 and C(t = 0) = Ca,in .

2. Repeat the above integration, for the case of a variable entry flowrate: Ca,in = 0, 4 mol/l
for t ∈ [0, 20) min; Ca,in = 0, 5 mol/l for t ≥ 20 min. The final integration time will be set
equal to 40 min in this case.

3. Plot the evolution of the molar concentration of A with respect to time. Add the proper
figure elements (titles, legends, etc.).

Case of a cascade of three CSTRs in series in dynamic state

dCA1 = 1 (C
dt τ a,in − CA1 ) − kCA1
dCA2 = 1 (C − C ) − kC
dt τ A1 A2 A2

dCA3 = 1 (C − C ) − kC
dt τ A2 A3 A3

In the above system, Ai is the molar concentration of A in the reactor "i".

1. Repeat the previous cases for constant and variable entry flowrate.

2. Plot the dynamic evolution of the molar concentration of A in the different reactors, in a
single figure.

PSE - MATLAB Modeling Page 7

 Case of a Plug Flow Reactor (PFR) in dynamic state The temporal-spatial evolution of
the concentration of species A in a PFR can be described by a Partial Differential Equation of
the form:
∂CA ∂CA
+u + k CA = 0 (21)
∂t ∂z
subject to the following boundary condition:

C(z = 0, t = 0) = Ca,in

An alternative solution to solving the above PDE is to consider a cascade of infinite CSTRs
in series.

1. Extend the system of ODEs to N CSTRs using: u = 2 m/min and an overall space time
for the PFR equal to 10 min. At a first approximation, N can be considered equal to 50.
These elements allow calculating the overall reactor length, L, as well as the correspond-
ing length, dz, and space time, τi , of each CSTR.

2. Plot the dynamic evolution of the molar concentration of A in any given reactor ’j’.

3. Plot the spatial evolution of the molar concentration of A along the overall length of the
PFR.

PSE - MATLAB Modeling Page 8

Problem 6: Optimization

Generalities

MATLAB proposes in its toolbox a variety of functions for solving both constrained and un-
constrained optimization problems. Among these functions, "fminsearch" and "fmincon" are
specifically powerful and can be employed on the basis of different numerical methods.

Work to be done

The Rosenbrock function Let us first consider the function of Rosenbrock (a.k.a. the "banana"
function, due to its shape), defined by Eq. 22:

f (x) = (1 − x1 )2 + 100(x2 − x21 )2 (22)

1. Use the internal function "fminsearch" to locate the minimum of f and the corresponding
value of x.

2. We wish to constrain the search domain to the following intervals:

• x ∈ [−2, 2]x[−3, 3].

• x ∈ [−2, 0]x[−3, 0].

Use "fmincon" to achieve this. What do you observe?

3. Extend the unconstrained optimization problem to the general form of the same function:

f (x) = (α − x1 )2 + 100(x2 − x21 )2 (23)

Test different values of α, which should be defined in the main script.

The Himmelblau function The Himmelblau function is defined by Eq. 24:

f (x) = (x21 + x2 − 11)2 + (x1 + x22 − 7)2 (24)

1. Plot the function in a 3D plot for x ∈ [−5, 5]x[−5, 5]. What do you observe with respect to
the searched minima?

2. Use fminsearch to minimize f . Several initial values will be tested in order to increase the
possibility to locate all local minima. In this respect, a number of random initial points will
be created automatically (e.g., 50).

3. Plot on a 2D figure the coordinates of the located minima.

Parametric identification The vapor pressure of a fluid can be expressed as a function of

the temperature, according to a correlation of the form 25:

a2
ln(Psat ) = a1 + + a4 T + a5 T 2 + a6 (25)
a3 + T

PSE - MATLAB Modeling Page 9

 The text file "Pdata.res" contains experimental values of the vapor pressure with respect
to temperature (T in Kelvin, Psat in Pa).

1. Determine the values of the coefficients ai using different starting points.

2. Set the values of coefficients a4 and a5 equal to 0 and repeat the above procedure. Pro-
pose different means to reduce the correlation. Comment.

3. Repeat, specifying the following starting point: a = [11, 1000, 1, 0].

4. Plot the data and the correlation on the same plot.

Useful commands: surf, rand. Expected solutions: Global minimum of the Rosenbrock function
@[α, α2 ]. Local minima of the Himmeblau function @[−2.80, 3.13]; [3.0, 2.0]; [−3.78, −3.28]; [3.58, −1.85].
T −15.1 .
1324.3
The optimal reduced correlation should be: ln(Psat ) = 16.6 −

PSE - MATLAB Modeling Page 10