3.

Mathematical Tools for Dynamic Modeling

3.1 Introduction
In Lect. 2, a number of mathematical dynamic models (linear/nonlinear ODEs) has been developed
to describe the dynamic operation of selected processes. The solution of such models (linear/nonlinear
ODEs) ,finding the output variables as functions of time for some change in the input variables, requires
either an analytical approach (to solve linear ODEs) or a numerical approach (to solve linear/nonlinear
ODEs) in order to integrate the differential equations (ODEs).

Linear systems (their dynamic models described by linear ODEs) represent the starting point for
many analysis techniques in process control. For linear ODEs, it is also possible to generate analytical
(such as Laplace Transforms) solutions. The important advantage of an analytical solution is that the
effects of changing the model structure or the model parameters are readily apparent. By contrast, for a
numerical solution, this information can only be obtained by recalculating the solution for each situation of
interest. The analytical solutions to ODEs can also be generated using software (such as MATLAB) for
symbolic problems.
1
 3.2 Laplace Transforms

3.2.1 Definition
The Laplace transform of a function f(t) is defined as
∞
ʆ𝒇 𝒕 = 𝑭 𝒔 = ‫𝒆 𝒕 𝒇 𝟎׬‬−𝒔𝒕 𝒅𝒕 (1)
where F(s) is the symbol for the Laplace transform, s is a complex independent variable, f(t) is some function of time to be
transformed, and ʆ is an operator, defined by the integral. The function f(t) must be piecewise continuous for 0< t <∞; this
requirement almost always holds for functions that are useful in process modeling and control. The inverse Laplace transform (ʆ-1)
operates on the function F(s) and converts it to f(t).
One of the important properties of the Laplace transform and the inverse Laplace transform is that they are linear operators; a
linear operator satisfies the Principle of Superposition:
ʆ(ax(t) + by(t)) = a ʆ (x(t)) + b ʆ (y(t)) = aX(s) + bY(s) (2)
Advantages:
• Standard notation in dynamics and control
• Converts mathematics to algebraic operations
• Advantageous for block diagram analysis
Uses:
• Solution of differential equations (linear)
• Analysis of linear control systems (frequency response)
• Prediction of transient response for different inputs

2
3.2.2 Laplace Transforms of Most Simple Functions
• Constant function: f(t) = a (where: a is constant)
∞ 𝒂 𝒂 𝒂
ʆ 𝒂 = ‫𝒆𝒂 𝟎׬‬−𝒔𝒕 𝒅𝒕 = (− ) 𝒆−𝒔𝒕 ∣∞
𝟎 = 0 - (− ) = (3)
𝒔 𝒔 𝒔

• Step function: S(t) = 0 for t<0 and S(t) = 1 for t≥0

(called unit step function, used frequently in process dynamics and control)
𝟏
ʆ S(t) = ʆ 1 = (4a)
𝒔
𝒂
If the step magnitude is a (i.e. aS(t)), ʆ aS(t) = (4b)
𝒔

• Exponential Functions: For an exponential, e−bt, with b > 0

∞ ∞ −𝟏 𝟏
ʆ e−bt = ‫ 𝟎׬‬e−bt𝒆−𝒔𝒕 𝒅𝒕 = ‫ 𝟎׬‬e−(b+s)t 𝒅𝒕 = ( )e−(b+s)t ∣∞
𝟎 = 0 - (− )
𝒃+𝒔 𝒃+𝒔
𝟏
= (5)
𝒔+𝒃

• Ramp Function: f(t)= t

∞
ʆ t = ‫׬‬0 t𝑒 −𝑠𝑡 𝑑𝑡 [integration by parts; u= 𝑡 and dv=𝑒 −𝑠𝑡 dt]
−𝒕𝒆−𝒔𝒕 ∞ 𝟏 ∞ 𝟏 𝟏
ʆ t =
𝒔
∣𝟎 +
𝒔
‫𝒆 𝟎׬‬−𝒔𝒕 𝒅𝒕 = 𝟎 + (𝒔𝟐)(𝟎 − (−𝟏))= 𝒔𝟐 (6)

3
• Rectangular Pulse Function:
An illustration of the rectangular pulse is shown below. The pulse has height h and width tw (area under the curve is htw).
Mathematically, the function f(t) is defined as:
f (t) = 0 t<0
= h 0 ≤ t < tw
= 0 t ≥ tw
∞ 𝒕 𝒉 𝒉 𝒕
F(s)=ʆ 𝒇(𝒕) = ‫𝒆)𝒕(𝒇 𝟎׬‬−𝒔𝒕 𝒅𝒕 = ‫𝒆𝒉 𝒘 𝟎׬‬−𝒔𝒕 𝒅𝒕 = (− ) 𝒆−𝒔𝒕 ∣𝒕𝟎𝒘 = (1-e- 𝒘 s) (7)
𝒔 𝒔

For a unit rectangular pulse, h = 1/tw and the area under the pulse is unity.

• Impulse Function:
A limiting case of the unit rectangular pulse is the impulse or Dirac delta function, which has the symbol δ(t). This function
is obtained when tw → 0 while keeping the area under the pulse equal to unity. A pulse of infinite height and infinitesimal
width results. Mathematically, this can be accomplished by substituting h = 1/tw into pulse Eqn.; the Laplace transform of δ(t)
𝟏
is: ʆ δ(t) = 𝒍𝒊𝒎 𝒕 𝒔 𝟏 − 𝒆−𝒕𝒘𝒔 (L’Hôpital Rule, which involves taking derivatives of both numerator and denominator)
𝒕𝒘→𝟎 𝒘

𝒔𝒆−𝒕𝒘𝒔
=𝒍𝒊𝒎 =1 (8)
𝒕𝒘→𝟎 𝒔

If the impulse magnitude (i.e., area twh) is a constant (a) rather than unity, then: (aδ(t)) = a)..
Rapid injection of dye or tracer into a fluid stream is an example of impulse function..

4
• Derivatives: df/dt, d2f/dt2, d3f/dt3, etc.
∞
ʆ df/dt = ‫ 𝟎׬‬df/dt 𝒆−𝒔𝒕 𝒅𝒕 [integration by parts; u= 𝒆−𝒔𝒕 and dv=(df/dt)dt]
∞
ʆ df/dt = ‫ 𝟎׬‬f(t)𝒆−𝒔𝒕 𝒔𝒅𝒕 + f(t) 𝒆−𝒔𝒕 ∣∞
𝟎 = s ʆ(f(t)) − f(0) = sF(s) − f(0) (9a)
ʆ d2f/dt𝟐 = ʆ d𝜱/dt = s𝜱(s) − 𝜱(0) (where: 𝜱= df/dt)
= s[sF(s) − f(0)] - f’(0) (where: f’=df/dt)
= s2F(s) - sf(0) - f’(0) (9b)
𝒏−𝟏 𝒏−𝟐 ( )( 𝒏−𝟐 𝒏−𝟏
For more general, ʆ dnf/dtn = 𝒔𝒏𝑭 𝒔 − 𝒔 𝒇 𝟎 −𝒔 𝒇 𝟏 𝟎) − · · · −𝒔𝒇 (𝟎) − 𝒇 (𝟎) (9c)
where f (i)(0) is the ith derivative evaluated at t = 0.

• Cosine and Sine functions: cos(wt), sin(wt)

𝒆𝒋𝒘𝒕 +𝒆−𝒋𝒘𝒕 𝟏 𝟏 𝟏 𝟏 𝒔−𝒋𝒘 𝒔+𝒋𝒘 𝒔
ʆ(cos(wt))= ʆ( )=( )( + )=( )( + )= j2= −1 (10a)
𝟐 𝟐 𝒔+𝒋𝒘 𝒔−𝒋𝒘 𝟐 𝒔𝟐 +𝒘𝟐 𝒔𝟐 +𝒘𝟐 𝒔𝟐 +𝒘𝟐

𝒆𝒋𝒘𝒕 −𝒆−𝒋𝒘𝒕 𝟏 𝟏 𝟏 𝟏 𝒔−𝒋𝒘 𝒔+𝒋𝒘 𝒘

ʆ(sin(wt))= ʆ( )=( )( - )=( )( - )= (10b)
𝟐𝒋 𝟐𝒋 𝒔+𝒋𝒘 𝒔−𝒋𝒘 𝟐𝒋 𝒔𝟐 +𝒘𝟐 𝒔𝟐 +𝒘𝟐 𝒔𝟐 +𝒘𝟐

Note: 𝒆−𝒋𝒘𝒕 = cos(wt) – jsin(wt) 𝒆𝒋𝒘𝒕 = cos(wt) + jsin(wt) j= −𝟏 (10c)

The Laplace Transforms for Various Time-Domain Functions can be seen in Table (3.1) (next 3 slides)
5
Table (3.1) Laplace Transforms for Various Time-Domain Functions

6
Table (3.1) (Contd.)
f(t) F(s)

cos

7
Table (3.1) (Contd.)
f(t) F(s)

8
3.2.3 Solution of Differential Equations by LAPLACE Transform Techniques
• Laplace Transform of Low Order Differential Equations
The procedure used to solve low order differential equation is quite simple. First, Laplace transform both sides of the
differential equation, substituting values for the initial conditions in the derivative transforms. Rearrange the resulting
algebraic equation, and solve for the transform of the dependent (output) variable. Finally, find the inverse of the
transformed output variable [matching one of the entries in table (3.1)]. The procedure can be explained by example 1.

Ex1: Solve the differential equation, 5dy/dt+ 4y = 2 y(0) =1 using Laplace transforms.
Soln.
L(5dy/dt+ 4y)=L(2) (take the Laplace transform of both sides of ODE):
5sY(s)-5y(0)+4Y(s)=2/s
(5s+4)Y(s)=5+(2/s) [collect all Y(s) terms to the left side and other terms to the right side]
5𝑠+2
Y(s)= [multiply both sides by s and solve for Y(s)]
𝑠(5𝑠+4)
5𝑠+2
L-1 [Y(s)]= L-1 [ ] [take the inverse Laplace transform of both sides and solve for y(t)]
𝑠(5𝑠+4)
𝑠+0.4
y(t)= L-1 [ ] [divide the numerator and denominator by 5 to put all factors in the s + b form]
𝑠(𝑠+0.8)

So, y(t)= 0.5 + 0.5e−0.8t [see table (entry11), (s+b3)/(s+b1)(s+b2), b3=0.4, b1=0.8, b2=0]

9
• Laplace Transform of High Order Differential Equations [Partial Fraction Expansion (PFE)]:
The high-order denominator polynomial in a Laplace transform solution arises from the differential equation terms (its
characteristic polynomial) plus terms contributed by the inputs. The factors of the characteristic polynomial correspond to
the roots of the characteristic polynomial set equal to zero. The input factors may be quite simple. Once the factors are
obtained, the Laplace transform is then expanded into partial fractions.
In general, there are 4 cases:

Case 1: real and distinct factors

𝑵(𝒔) 𝑵(𝒔) α𝒊
In this case, the following expansion formula can be used: 𝒀(𝒔) = = ς𝒏 = σ𝒏𝒊=𝟏 (11a)
𝑫(𝒔) 𝒊=𝟏(𝒔+𝒃𝒊) 𝒔+𝒃𝒊
𝑵(𝒔)
The ith coefficient can be calculated using the Heaviside expansion: αi =(𝒔 + 𝒃𝒊) | (11b)
𝑫(𝒔) s=−bi
𝑵′(𝒔) 𝑵′(𝒔) 𝜶′𝒊
Also, an expansion for real/distinct factors may be written as: 𝒀(𝒔) = = ς𝒏 = σ𝒏𝒊=𝟏 (12a)
𝑫′(𝒔) τ
𝒊=𝟏( i𝒔+𝟏) τi𝒔+𝟏
𝑵′(𝒔)
The ith coefficients can be calculated also using the Heaviside expansion: 𝜶′𝒊 = (τi𝒔 + 𝟏) 𝑫′(𝒔) |s=-1/τi (12b)

10
 𝑠+5
Ex2: If 𝑌 𝑠 = , find y(t)?
𝑠2+5𝑠+4
Soln.
PFE Procedure

The denominator can be factored into a product of first-order terms, (s + 1)(s + 4). This transform can be expanded into
the sum of two partial fractions:
𝑠+5 α α
Y(s)= = 1+ 2
(𝑠+1)(𝑠+4) 𝑠+1 𝑠+4

The Heaviside expansion method gives: α1 = 4/3, α2 = −1/3.

Invert each term individually using Table (3.1):
4/3 −1/3 1 1
L-1 [Y(s)] = L-1 [ + ] = (4/3) L-1 [ ] - (1/3) L-1 [ ]
s+1 s+4 s+1 s+4
y(t) = (4/3)e-t – (1/3)e-4t

11
Ex3: Solve the ordinary differential equation (ODE): d3y/dt3 + 6 d2y/dt2 + 11dy/dt+ 6y = 1 with initial conditions
y(0) = y′(0) = y′′(0) =0 where the primes denote derivatives.
Soln.
s3Y(s) + 6s2Y(s) + 11sY(s) + 6Y(s)=1/s (Take the Laplace transform of both sides of ODE)
Y(s)=1/[s(s3 + 6s2 + 11s + 6)] (Solve for Y(s))
Y(s)=1/[s(s+1)(s+2)(s+3)] (Factor the denominator into a product of first-order terms)
α1 α α3 α4
= + 2 + + (Expanded into the sum of partial fractions)
𝑠 𝑠+1 𝑠+2 𝑠+3
The Heaviside expansion method gives: α1 = 1/6, α2 = −1/2, α3 = 1/2, α4 = −1/6.
After the transform has been expanded into a sum of first-order terms, invert each term individually using Table (3.1):
1/6 −1/2 1/2 −1/6
L-1 [Y(s)] = L-1 [ + + + ]
𝑠 𝑠+1 𝑠+2 𝑠+3
= (1/6) L-1 [1/s] – (1/2) L-1 [1/(s+1)] + (1/2) L-1 [1/(s+2)] – (1/6) L-1 [1/(s+3)]
y(t) = (1/6) – (1/2)e-t + (1/2)e-2t – (1/6)e-3t

12
• The general procedure for solving both low and high order ODEs using Laplace Transforms consists of four steps, as
shown in the following figure:

13
Case 2: real and repeated factors
𝑵(𝒔)
A general method for repeated factors is described in this section. If a function Y(s) has the form: 𝒀 𝒔 = 𝑫(𝒔)
The required PFE has the general form (including both non-repeated factors and repeated factor (s+b)n):
𝑵(𝒔) 𝜶𝒏 𝜶𝒏 𝜶𝟏
𝒀 𝒔 = =··· + +···+ + [other partial fractions of non-repeated factors of D(s)] (13)
𝑫(𝒔) (𝒔+𝒃)𝒏 (𝒔+𝒃)𝒏−𝟏 𝒔+𝒃

The coefficients of non-repeated factors is calculated by the previous procedure (case 1). To calculate the coefficients of
repeated factor (𝛼1 , 𝛼2 , ……., 𝛼𝑛 ), the general PFE form will be multiplied by (s+b)n as below:
𝑵 𝒔
Q(s)= 𝒔 + 𝒃 𝒏 = (s + b)n−1 α1 + (s + b)n−2 α2 + · · ·+ αn + (s + b)n × [other partial fractions of non-repeated factors] (14)
𝑫 𝒔

The coefficients of repeated factors is calculated by the following general expression:

𝟏 𝐝(𝐢) 𝐐(𝐬)
𝛂𝐧−𝐢 = ∣s=-b (i=0, 1, 2,…….., n-1) (15)
𝐢! 𝐝𝐬 (𝐢)
𝒔+𝟏
Ex.4: Calculate y(t) if 𝒀 𝒔 =
𝒔 𝒔+𝟐 𝟐
Soln.
𝑠+1 𝛼 𝛼1 𝛼3 𝛼2 𝛼1 𝛼3 1 1 1
𝑌 𝑠 = = 2 + + So, y(t) =L-1 [ + + ] = 𝛼1 L-1 [ ]+ 𝛼2 L-1 [ ]+ 𝛼3 L-1 [ ]
𝑠 𝑠+2 2 (𝑠+2)2 (𝑠+2) 𝑠 (𝑠+2)2 (𝑠+2) 𝑠 (𝑠+2) (𝑠+2)2 𝑠
𝒔+𝟏
𝛼3 (coefficient of non-repeated factor) is calculated by the procedure of case1 (eqn. 11b): 𝜶𝟑 = │𝒔 = 𝟎 = ¼
(𝒔+𝟐)𝟐
𝟐 𝒔+𝟏
𝛼1 and 𝛼2 (coefficients of repeated factor) are calculated by the procedure of case2 (eqns. 14 and 15): 𝑸 𝒔 = 𝒀 𝒔 ∗ 𝒔 + 𝟐 =
𝒔
𝒔+𝟏 𝒔+𝟏 −𝟏
𝒊 = 𝟎, 𝜶𝟐 = │𝒔 = −𝟐 = 1/2 𝒊 = 𝟏, 𝜶𝟏 = d( )/ds │𝒔 = −𝟐 = │𝒔 = −𝟐 = -1/4
𝒔 𝒔 𝒔𝟐
𝟏 𝟏 𝟏 𝟏
y(t)= (-1/4) L-1 [ ]+ ( )L-1 [ ]+ (1/4)L-1 [ ] = (-1/4)e-2t +(1/2)te-2t + (1/4) (see table 3.1)
(𝒔+𝟐) 𝟐 (𝒔+𝟐)𝟐 𝒔

14
Case 3: complex factors
An important special case for PFE occurs when Y(s) [N(s)/D(s)] is factored and has a factor of the form:
𝒅𝟐𝟏
Y(s)= (c1s + c0)/(s2 + d1s + d0) with < d0 (16)
𝟒

Then the denominator of Eq. 16 has complex roots and Y(s) is said to have complex factors.
The general procedure for performing a PFE when there is a complex factor can be summarized as below. First,
convert Eq. 16 to the standard form (completing the square):
Y(s)= (c1s + c0)/(s2 + d1s + d0) = [c1(s + b) + (c0-c1b)]/[(s + b)2 + ω2] (17a)
𝒅𝟐𝟏
Where: b=d1/2 and w= d0 − 𝟒
(17b)

So, y(t) = L-1 [Y(s)] = c1 e-bt cos(wt) + [(c0-c1b)/w] e-bt sin(wt) (see table 3.1, entries 17 and 18) (18)
Ex.5: If y(s)=(s+2)/(s2+s+1), calculate y(t)?
Soln.
(note: c1=1, c0=2, d1=1, d0=1)
b=0.5, w= 𝟑/𝟐 (see eqn. 17b)
Y(s) = (s + 2)/[(s + 0.5)2 − 0.25 + 1] = [(s + 0.5) + 1.5]/[(s + 0.5)2 +( 𝟑/𝟐)𝟐] (see eqn. 17a)
So, y(t) = L-1 [Y(s)] = e-0.5t cos( 𝟑 𝒕/𝟐) + 𝟑 e-0.5t sin( 𝟑 𝒕/𝟐) (see eqn. 18)

15
Case 4: (both repeated and complex roots/factors show up):
This case can be explained by the following example.
Ex.6: If Y(s) = (s + 1)/[s2(s2 + 4s + 5), calculate y(t)?
Soln.
The denominator of Y(s) has two repeated roots at s = 0. Also, the roots of s 2 + 4s + 5 are complex because the inequality in Eq.
𝒅𝟐𝟏
16 ( < d0) is satisfied. The appropriate partial fraction (PFE) for Y(s) is:
𝟒
Y(s) = (s + 1)/[s2(s2 + 4s + 5) = α1/s+ α2/s2 + (α3s + α4)/(s2 + 4s + 5) (19a)
Multiply both sides of Y(s) by s2(s2 + 4s + 5) and collect terms:
s + 1 = (α1 + α3)s3 + (4α1 + α2 + α4)s2 + (5α1 + 4α2)s + 5α2 (19b)
Equate coefficients of like powers of s:
s3∶ α1 + α3 = 0 (20a) s2∶ 4α1 + α2 + α4 = 0 (20b) s1∶ 5α1 + 4α2 = 1 (20c) s0∶ 5α2 = 1 (20d)
Solving simultaneously gives: α1 = 0.04, α2 = 0.2, α3 = −0.04, and α4 = −0.36.
So, y(t)= L-1[0.04/s+ 0.2/s2 - (0.04s + 0.36)/(s2 + 4s + 5)]= L-1[0.04/s]+ L-1[0.2/s2]- L-1[(0.04s + 0.36)/(s2 + 4s + 5)] (21)
Note: completing the square for the last term must be done:
(0.04s + 0.36)/(s2 + 4s + 5)=(0.04s + 0.36)/(s2 + 4s + 4 + 5 - 4)= [0.04(s+2) + (0.36-0.08)]/[(s+2)2 + 1] )= [0.04(s+2) + 0.28]/[(s+2)2 + 1]
So, y(t)=0.04+0.2t-L-1[0.04(s+2) + 0.28]/[(s+2)2 + 1]
= 0.04+0.2t-0.04 e-2t cos(t)- 0.28 e-2t sin(t)

16
3.2.4 Other Laplace Transform Properties
• Final Value Theorem:
The asymptotic value of y(t) for large values of time y(∞) can be found from the following eqn.:
y(∞) = lim 𝒚 𝒕 = 𝒍𝒊𝒎 𝒔𝒀 𝒔 (22)
𝒕→∞ 𝒔→𝟎

• Initial Value Theorem:

Analogous to the Final Value Theorem, the Initial Value Theorem can be stated as:
y(0) = 𝒍𝒊𝒎 𝒚 𝒕 = lim 𝒔𝒀 𝒔 (23)
𝒕→𝟎 𝒔→∞

Ex.7: Apply the Initial and Final Value Theorems to this function: Y(s) = (5s + 2)/[s(5s + 4)]
Soln.
Final Value Theorem: y(∞) = 𝒍𝒊𝒎 𝒔𝒀 𝒔 = 𝒍𝒊𝒎 𝒔(5s + 2)/[s(5s + 4)] = 𝒍𝒊𝒎 (5s + 2)/(5s + 4) = 0.5 (using eqn. 22)
𝒔→𝟎 𝒔→𝟎 𝒔→𝟎
Initial Value Theorem: y(0) = 𝐥𝐢𝐦 𝒔𝒀 𝒔 = 𝐥𝐢𝐦 (5s + 2)/(5s + 4) = 1 (using eqn. 23)
𝒔→∞ 𝒔→∞

17
Ex.8: A process is described by a third-order ODE: d3y/dt3 + 6d2y/dt2 + 11dy/dt+ 6y = 4u , with all initial
conditions on y, dy/dt, and dy2/dt2 equal to zero. Show that for a step change in u of two units, the steady-state
result in the time domain is the same as applying the Final Value Theorem.
Soln.
(s3 + 6s2 + 11s + 6)Y(s) = 8∕s (Taking Laplace for both sides of ODE, u=2)
Y(s)=8/[s(s3 + 6s2 + 11s + 6)] (Solve for Y(s))
One of the benefits of the Final Value Theorem is that we do not have to solve for the analytical solution of y(t). Instead,
simply apply Eq. 22 to the transform Y(s) as follows: y(∞) = 𝒍𝒊𝒎 𝒔𝒀 𝒔 = 𝒍𝒊𝒎 𝒔𝒀 𝒔 = 8/(s3 + 6s2 + 11s + 6) = 4/3
𝒔→𝟎 𝒔→𝟎
The time domain solution obtained from a PFE is:
y(t) = L-1[ 8/[s(s3 + 6s2 + 11s + 6)] ] = L-1[ 8/[s(s+1)(s+2)(s+3] ] = L-1[a1/s+a2/(s+1)+a3/(s+2)+a4/(s+3)]
note: a1= 4/3, a2= -4, a3= 4, a4= -4/3
= 4∕3 − 4e−t + 4e−2t − 4∕3e−3t
As t → ∞, only the first term remains (4/3), which is the same result as the one obtained from the Final Value Theorem.
Check:
At steady state (u = 2 while all derivatives =0), so 6y =8 or y = 4∕3 (Steady-state value (final value) of y)

18
• Transform of an Integral:
Occasionally, it is necessary to find the Laplace transform of a function that is integrated with respect to time. The
𝒕
transform of an integral is given by: L{‫ 𝟎׬‬f (t∗)dt∗}= (1/s)F(s) (24)

• Time Delay (Translation in Time, θ):

Time delays play an important role in process dynamics and control. They commonly occur as a result of the transport
time required for a fluid to flow through a pipe (θ=L/v).
fd(t) = f (t − θ)S(t − θ) {fd(t) is a f(t) delayed by θ time units) (25)
L[f (t − θ)S(t − θ)]=e-θs F(s) (26a)
Or f(t − θ)S(t − θ) = L-1 [e-θs F(s)] (26b)

EX. 9: Find the inverse transform of Y(s) = (1 + e−2s)/[(4s + 1)(3s + 1)]

Soln.
Y(s)=1/[(4s + 1)(3s + 1)] + e−2s/[(4s + 1)(3s + 1)]
y(t)=L-1[ 1/[(4s + 1)(3s + 1)] ] + L-1 [ e−2s/[(4s + 1)(3s + 1)] ]
PFE: 1/[(4s + 1)(3s + 1) = a1/(4s+1)+a2/(3s+1), (a1=4, a2=-3)
y(t)= L-1[ 1/(s+1/4) – 1/(s+1/3) ] + L-1[e−2s /(s+1/4) – e−2s/(s+1/3) ]
= e−t/4 - e-t/3 + (e−(t-2)/4 - e-(t-2)/3) S(t-2)

19
EX. 10: A TRANSIENT RESPONSE EXAMPLE
CSTR is operated at a temp. below its reaction temp. (no reaction occurred) to check the exit concentration response when pulse
and impulse functions has been occurred. The Stirred-Tank Reactor System and Operating Conditions are the following: Volume =
4 m3, Flow rate q = 2 m3/min, Nominal feed concentration ci = 1 kg mol/m3
Find the exit concentration response if a pulse function is used for feed concentration as shown in the following figure?
Soln.
The dynamic model of this CSTR is:
Vdc/dt=q(ci-c) (ignored the reaction term)
4dc/dt+2c=2ci (c(0)=ci(0)=1, initially at steady state)
4sC(s)-4(1)+2C(s)=2Ci(s)
or (2s+1)C(s) = 2+Ci(s)
So, C(s) = 2/(2s+1) +Ci(s)/(2s+1)
Using Pulse input: Ci(s)=1/s+(5/s)(1-e-0.25s)
C(s) = 2/(2s+1) +1/[s(2s+1)]+{5/[s(2s+1)]}(1-e-0.25s)
=2/(2s+1) +6/[s(2s+1)]-{5/[s(2s+1)]}e-0.25s
c(t)= L-1[ 2/(2s+1) + 6/[s(2s+1)] - {5/[s(2s+1)]}e-0.25s ]
c(t)= e-t/2 + 6(1 − e−t∕2) – 5(1 − e−(t-0.25)∕2)S(t-0.25)

20
 3.3 Solving Symbolic Mathematical Problems Using MATLAB
Analytical solutions of dynamic model responses based on derivations using the Laplace transform method
have been considered in the previous section. Analytical solutions provide greater insight into the structure
of the solution and its sensitivity to changes in model parameters but they are restricted to linear models and
only a few special nonlinear models. By contrast, solutions of most nonlinear dynamic models can only be
obtained using numerical methods.

Analytical solutions to linear dynamic models can also be obtained using software for symbolic mathematical
problems such as MATLAB. MATLAB is capable of symbolic processing.

• Using MATLAB for Symbolic Processing:

To prepare MATLAB for symbolic operations, some variable names will be declared symbolic (rather than
numeric) using the syms command.

Ex: syms t s a b c d w

This means that MATLAB can give solution as a function of these variables: t, s, a, b, c, d, w

1
• Laplace Transforms of Simple Functions using MATLAB:
laplace command is used to calculate the Laplace (L) for different functions of time as shown below:

MATLAB Code for Heaviside Function [the unit step, heaviside(t)]:

>> syms s t
>> laplace (heaviside(t))
ans =
1/s

MATLAB Code for Dirac Function [the unit impulse function, dirac(t)]:
>> laplace (dirac(t))
ans =
1

MATLAB code for Exponential Function [exp(-b*t)]:

>> syms b
>> laplace (exp(-b*t))
ans =
1/(b + s)
2
MATLAB code for Ramp Function (t):
>> laplace(t)
ans =
1/s^2

MATLAB code for sin Function [sin(w*t)]:

>> syms w
>> laplace(sin(w*t))
ans =
w/(s^2 + w^2)

MATLAB code for cos Function [cos(w*t)]:

>> laplace(cos(w*t))
ans =
s/(s^2 + w^2)

3
• Inversion of Laplace Transforms Using MATLAB:
ilaplace command is used to calculate the Laplace inverse (L-1).

Ex: Use the ilaplace Command in the MATLAB Symbolic Math Toolbox to find the inverse Laplace transform
of Y(s) = (s + 1)/[(s + a)(s + b)] where a and b do not have numerical values.
Soln.
syms a
>> ilaplace((s+1)/((s+a)*(s+b)))
ans =
(exp(-a*t)*(a - 1))/(a - b) - (exp(-b*t)*(b - 1))/(a - b)

EX: Use the ilaplace Command in the MATLAB Symbolic Math Toolbox to find the inverse Laplace transform
of a function: Y(s) = (s + 2)/( s2 + s + 1)
Soln.
>> ilaplace((s + 2)/( s^2 + s + 1))
ans =
exp(-t/2)*(cos((3^(1/2)*t)/2) + 3^(1/2)*sin((3^(1/2)*t)/2))

4
EX: Use the ilaplace Command in the MATLAB Symbolic Math Toolbox to find the inverse Laplace transform
of a function: Y(s) = (s + 1)/(s2 + 4s + 4) = (s + 1)/( s + 2)2
Soln.
>> ilaplace((s + 1)/( s + 2)^2)
ans =
exp(-2*t) - t*exp(-2*t)
EX: Using MATLAB, Calculate c(t) if
a- C(s) = (1 + 3.25s)/[s(2s + 1)] b- C(s) =[1/(2s + 1)](2 + 6/s− 5e−0.25s/s)
Soln.
a- >> ilaplace((1 + 3.25*s)/(s*(2*s + 1)))
ans =
(5*exp(-t/2))/8 + 1

b- >> ilaplace((1/(2*s + 1))*(2 + 6/s-(5*exp(-0.25*s))/s))

ans =
5*heaviside(t - 1/4)*(exp(1/8 - t/2) - 1) - 5*exp(-t/2) + 6
5
• Solution of Differential Equations Using MATLAB Symbolic Processing:
MATLAB can solve differential equations symbolically using the DSOLVE command. The derivatives
are represented as Dx (first derivative), D2x (second derivative), etc.
EX: Solve the following ODE using MATLAB (dsolve command)
d3x/dt3 + 4d2x/dt2 + 5dx/dt + 2x=2 using the following initial conditions: x(0) = dx(0)/dt = d2x(0)/dt2 = 0
Soln.
>> x=dsolve('D3x+4*D2x+5*Dx+2*x=2','x(0)=0','Dx(0)=0','D2x(0)=0')
x=
1 - 2*t*exp(-t) - exp(-2*t)

• Simplifying Complex symbolic expressions in MATLAB:

Complex expressions can be simplified in MATLAB using commands such as collect, expand,
and simplify.
EX: Use expand command to collect the terms on the right side of y = 2(x + 1)2 + 4(x + 2)3 + 7
Soln.
syms x
y = 2*(x + 1)^2 + 4*(x + 2)^3 + 7
expand (y)
ans = 4*x^3 + 26*x^2 + 41
6
EX: Use simplify command to simplify the following: y=sin2x+cos2x
Soln.
>> syms x
>> simplify(sin(x)^2+cos(x)^2)
ans =
1

EX: Use collect command to collect the coefficients of the following expression (x2y+yx-x2-2x) in powers of
x and in powers of y.
Soln.
>> syms x y
>> coeffs_x=collect(x^2*y+y*x-x^2-2*x,x)
coeffs_x =
(y - 1)*x^2 + (y - 2)*x
>> coeffs_y=collect(x^2*y+y*x-x^2-2*x, y)
coeffs_y =
(x^2 + x)*y - x^2 - 2*x
7
• Partial Fractions Using MATLAB Symbolic Processing:
MATLAB functions (diff and int) can be used together to do the PFE as shown in the following example
EX: Using MATLAB, find the partial fraction expansion of
Soln.
>> syms x
>> x=(s^4-6*s^2+9*s-8)/(s*(s-2)*(s+1)*(s+2)*(s-1))
x=
(s^4 - 6*s^2 + 9*s - 8)/(s*(s - 1)*(s + 1)*(s - 2)*(s + 2))
>> diff(int(x))
ans =
2/(3*(s - 1)) + 11/(3*(s + 1)) + 1/(12*(s - 2)) - 17/(12*(s + 2)) - 2/s
>> pretty(ans)

2 11 1 17 2
+ + − −
3(𝑠 − 1) 3(𝑠 + 1) 12(𝑠 − 2) 12(𝑠 + 2) 𝑠