034IN "Fundamentals of Automatic Control" - Introduction to MATLAB

Part 4L-1 - The Laplace Transform

Introduction

In this live script, we illustrate how to use the Matlab Symbolic Toolbox to determine the Laplace transform of a causal function of time and the inverse Laplace transform corresponding to a given rational, strictly proper or simply proper function of the complex variable s.

• How to compute the Laplace transform and the inverse transform using the Symbolic Math Toolbox.
• What are the Laplace transform main properties and theorems.

The Laplace Transform

Introduction and Motivations

• A mapping of a mathematical function from one domain to another.
• A change in perspective, not a change of the function.

• Some mathematical problems are difficult to solve in their natural domain.
• Transform to and solve in a new domain, where the problem is simplified.
• Solved the problem in the new domain, transform the solution back into the original domain.
• Trade off the extra effort of transforming/inverse-transforming for simplification of the solution procedure.
• An example: the phasor method.

• An integral transform mapping functions from the time domain to the Laplace domain or s-domain.

• Time-domain functions are functions of time t.
• Laplace-domain functions are functions of the complex frequency s

Motivations

[refer to the class material Part 2, pp. 28-30]

• Differential and integro-differential equations in the time domain: they can be difficult to solve.
• Apply the Laplace transform: transform to the s-domain.
• Differential and integro-differential equations become algebraic equations: easy to solve!
• Transform the s-domain solution back to the time domain.

Causal Functions

• We are interested in functions in the time domain defined for t ≥ 0, assuming they are null for all negative times, i.e
• Time domain functions that have this property are called causal functions.

Example: The Heaviside Function

clear

 

sympref('default');

disp(heaviside(sym([-10, -1/10000, 0, +1/10000, 10])))

tVALs = (-2:1e-2:2); % refer to Part 3 of "Introduction to MATLAB"

plot(tVALs, heaviside(sym(tVALs)), 'LineStyle', 'none', 'Marker','o', 'MarkerSize', 6, ...

'MarkerEdgeColor','r', 'MarkerFaceColor','r')

grid on; xlabel('time t [s]'); ylabel('f(t)');

title('Heaviside Function')

ylim([-0.25, 1.25])

Note the function value at (recall the default configuration of the Symbolic Math Toolbox, analyzed at the beginning of the live script Introduction to MATLAB Part 4 - Symbolic Math Computations ).

sympref('HeavisideAtOrigin', 1); % forcing the setting for the Heaviside function

 

disp(heaviside(sym([-10, -1/10000, 0, +1/10000, 10])))

Now the behavior of the Heaviside function is in accordance with the definition .

% tVALs = (-2:1e-2:2);

plot(tVALs, heaviside(sym(tVALs)), 'LineStyle', 'none', 'Marker','o', 'MarkerSize', 4, ...

'MarkerEdgeColor','b', 'MarkerFaceColor','b')

grid on; xlabel('time t [s]'); ylabel('f(t)');

title('Heaviside Function')

ylim([-0.25, 1.25]);

Another Example

Usually we will use the Heaviside function to emphasize that we are considering a causal function:

syms t

f_causal(t) = 3*sin(4*t+pi/3)*heaviside(t)

f_causal(t) = 

tVALs = (-2:1e-3:2);

plot(tVALs, f_causal(sym(tVALs)), 'LineStyle', 'none', 'Marker','o', 'MarkerSize', 4, ...

'MarkerEdgeColor','g', 'MarkerFaceColor','g')

grid on; xlabel('time t [s]'); ylabel('f(t)');

title('A Causal Function')

ylim([-3.25, 3.25]);

The Unilateral or One-Sided Laplace Transform

Definition

Given a causal function , the Laplace transform of is the function defined by

for those for which the integral makes sense.

• F is a complex-valued function of complex numbers;
• s is called the (complex) frequency variable, with units ; t is called the time variable (meas. units );
• the lower limit of integration is : the transform is called the unilateral or one-sided Laplace transform.
• Common notation convention: lower case letter denotes the time-domain function ; capital letter denotes its Laplace transform

Existence

• Fact 1: is continuous for ; if not, must contain no impulses at .
• Fact 2: any continuous, bounded, causal function admits Laplace transform.
• Fact 3: any piecewise continuous, causal function admits Laplace transform.

is piecewise continuous on .

clear t tVALs

syms t

 

assume(t>0)

f_pc(t) = symfun(piecewise( (0<=t & t<=1), t^2, ...

(1<t & t<=2), 3-t, ...

2<t , t+1 ), t);

 

tVALs = (0:1e-3:3); % 3001 values in [0, 3]

plot(tVALs, f_pc(sym(tVALs)), 'LineStyle','none', 'Marker','square',...

'MarkerSize',4, 'MarkerFaceColor','b', 'MarkerEdgeColor','b');

grid on;

xlabel('time t [s]'); ylabel('f_{pc}(t)');

title('A Piecewise Continuous Function')

is NOT piecewise continuous on .

f_Npc(t) = symfun(piecewise( (0<=t & t<=2), 1/(2-t), ...

2<t , t+1 ), t);

 

tVALs = (0:1e-2:3);

plot(tVALs, f_Npc(sym(tVALs)), 'LineStyle','-',...

'Color', 'm', 'LineWidth', 3);

grid on;

xlabel('time t [s]'); ylabel('f_{Npc}(t)');

title('A NOT Piecewise Continuous Function')

Remark on Laplace Transform and Impulses

• Indeed, it is possible to cope with impulses, but you have to modify the definition of the Laplace transform slightly.
• This topic is out of the scope of this course!

Examples of Laplace Transforms

Exponential Function

clear variables

 

syms t s

 

ft = exp(t)*heaviside(t);

 

Fs = laplace(ft, t, s)

Fs = 

Unit Step Function - the Heaviside Function

clear variables

 

syms t s

 

ht = heaviside(t);

 

Hs = laplace(ht, t, s)

Hs = 

Important Remark: On Using the Heaviside Function Notation For Causal Functions

• You can avoid inserting the Heaviside function explicitly to mark a function as causal unless you are dealing with delayed causal functions.

clear variables

 

syms t s

 

ft1(t) = exp(t)*heaviside(t)

ft1(t) = 

ft2(t) = exp(t)

ft2(t) = 

 

Fs1 = laplace(ft1, t, s)

Fs1 = 

Fs2 = laplace(ft2, t, s)

Fs2 = 

Laplace Transforms Properties

Linearity

• The Laplace transform is linear.
• If and are any causal functions, for which the Laplace transform exists, and a is any scalar, we have

clear variables

 

syms t s

syms f(t) g(t)

syms F(s) G(s)

syms a

 

 

omogeneityL = laplace(a*f(t))

omogeneityL = 

omogeneityL = subs(omogeneityL, laplace(f(t)), F(s))

omogeneityL = 

superpositionl = laplace(f(t)+g(t), t, s);

superpositionl = subs(superpositionl, [laplace(f(t)), laplace(g(t))], [F(s), G(s)])

superpositionl = 

clear variables

 

syms t s

f(t) = 1*heaviside(t);

g(t) = exp(t);

 

F(s) = laplace(f(t), t, s);

G(s) = laplace(g(t), t, s);

disp([F, G])

L(s) = laplace(3*f(t)+2*g(t), t, s)

L(s) = 

L(s) = collect(L, s)

L(s) = 

L(s) = simplify(L)

L(s) = 

Transform of a Derivative

If the causal function is continuous at , then

clear variables

 

syms t s

syms f(t) dot_f(t)

syms F(s)

 

dot_f(t) = diff(f(t), 1); % the 1st derivative

Dot_F(s) = laplace(dot_f, t, s)

Dot_F(s) = 

Dot_F(s) = subs(Dot_F(s), laplace(f(t), t, s), F(s))

Dot_F(s) = 

Let's verify the continuity of at :

F(s) = laplace(f, t, s) % Laplace transform of f(t)

F(s) = 

dot_f(t) = diff(f,t,1) % the 1st derivative of f(t)

dot_f(t) = 

DotF(s) = laplace(dot_f(t), t, s) % the Laplace transform of dot_f

DotF(s) = 

% the Laplace transform of the derivative of f(t)

DotF_ALT(s) = laplace(diff(f, t, 1), t, s)

DotF_ALT(s) = 

% the Laplace transform, applying the derivative property "by hand"

DotF_ALT2(s) = collect( (s*F(s)-f(0)), s)

DotF_ALT2(s) = 

Higher-Order Derivatives

• Applying derivative formula twice yields

• A similar formula holds for the k-th order derivative

clear variables

 

syms t s

syms f(t) dot_f(t) ddot_f(t)

 

dot_f(t) = diff(f, 1); % the 1st derivative

ddot_f(t) = diff(f, 2); % the 2nd derivative

Dot_F(s) = laplace(dot_f, t, s)

Dot_F(s) = 

DDot_F(s) = laplace(ddot_f, t, s)

DDot_F(s) = 

clear variables

 

syms t s

syms f(t) dot_f(t) ddot_f(t)

 

f(t) = 2*(exp(3*t)-1)

f(t) = 

dot_f(t) = diff(f, 1) % the 1st derivative

dot_f(t) = 

ddot_f(t) = diff(f, 2) % the 2nd derivative

ddot_f(t) = 

DDot_F(s) = laplace(ddot_f, t, s)

DDot_F(s) = 

 

DDotF_ALT(s) = s^2*laplace(f, t, s)-s*f(0)-dot_f(0)

DDotF_ALT(s) = 

collect(DDot_F, s)

ans(s) = 

Transform of an Integral

Let be the running integral of a causal function :

where .

clear variables

 

syms t s tau

 

assume(tau >=0)

assumptions

ans = 

 

syms f(t) g(t)

 

f(t) = heaviside(t);

F(s) = laplace(f, t, s)

F(s) = 

assumeAlso(t>0)

g(t) = int(f(tau), tau, 0, t)

g(t) = 

t

G(s) = simplify(laplace(g(t), t, s))

G(s) = 

Multiplication by t

Let be a causal function and define . Then , where

clear variables

 

syms t s

f(t) = exp(t);

g(t) = t*exp(t);

 

F(s) = laplace(f, t, s);

G(s) = laplace(g, t, s)

G(s) = 

 

G_ALT = -diff(F, s, 1)

G_ALT(s) = 

clear variables

 

syms t s

f(t) = exp(t);

F(s) = laplace(f, t, s);

f1(t) = t*exp(t);

g1(t) = t^2*exp(t);

F1(s) = laplace(f1, t, s);

G1(s) = laplace(g1, t, s)

G1(s) = 

 

G1_ALT = -diff(F1, s, 1)

G1_ALT(s) = 

 

G1_ALT_ALT = -diff(-diff(F, s, 1), s, 1)

G1_ALT_ALT(s) = 

Exponential Scaling (Frequency Shifting)

Let be a causal function, with Laplace transform , and let g(t) be defined as ( a is a scalar value)

• given then

clear variables

 

syms t s omega

 

f(t) = cos(omega*t)

f(t) = 

 

g(t) = exp(-t)*f(t)

g(t) = 

 

F(s) = laplace(f, t, s)

F(s) = 

 

G(s) = laplace(g, t, s)

G(s) = 

subs(F(s), omega, 1)

ans = 

collect(subs(G(s), omega, 1))

ans = 

Time Shifting (Time Delay)

• Let be a causal function, with Laplace transform .
• Define the causal function as

• Then:

clear variables

syms t s

 

a = sym(2); % the first delay term

b = sym(5); % the second delay term

 

f(t) = heaviside(t-a)-heaviside(t-b)

f(t) = 

% let's plot the pulse

tVALs = (-2:1e-3:10);

plot(tVALs, f(sym(tVALs)), 'LineStyle', 'none', 'Marker', 'square',...

'MarkerSize', 4, 'MarkerFaceColor', '#0072BD', 'MarkerEdgeColor', '#0072BD');

 

% fplot(f, [-2, 10], 'LineWidth', 3);

grid on

title('a rectangular pulse from t=2 to t=5')

xlabel('time t [s]');

ylabel('pulse f(t)');

ylim([-0.2, 1.2]);

 

% the corresponding Laplace transform

F(s) = collect(laplace(f, t, s), s)

F(s) = 

clear variables

 

syms t s

f(t) = ( 3 * exp(-0.75*t) * sin(16*t+pi/3) )*heaviside(t)

f(t) = 

 

T = sym(1.5); % the delay

 

% the delayed function

g(t) = f(t-T) % <<<<<<<-- Pay attention to the notation! <<<<<<<--

g(t) = 

g_ALT(t) = ( 3*exp(-0.75*(t-T)) * sin(16*(t-T)+pi/3) )* heaviside(t-T)

g_ALT(t) = 

% let's draw both the functions f and g

tVALs = (-1:1e-3:5);

plot(tVALs, f(sym(tVALs)), 'LineStyle', 'none', 'Marker', 'square',...

'MarkerSize', 4, 'MarkerFaceColor', '#0072BD', 'MarkerEdgeColor', '#0072BD');

hold on; grid on;

plot(tVALs, g(sym(tVALs)), 'LineStyle', 'none', 'Marker', 'square',...

'MarkerSize', 2, 'MarkerFaceColor', '#D95319', 'MarkerEdgeColor', '#D95319');

hold off

xlabel('time t [s]');

ylabel('f(t) and g(t)');

legend('f(t)', 'g(t)','Location', 'best', 'FontSize', 14)

title('A causal function vs. a delayed causal function')

ylim([-2.75, 3.25])

Convolution

The convolution of causal functions and is the causal function , defined as follows:

• Note: same as , i.e

• The (very great) importance of this property will soon become clear (refer to Part 2, pp. 27-32 of the course material).
• Laplace transform turns convolution in time domain into multiplication in the frequency domain!

The Inverse Laplace Transform

Given a Laplace transform , the unique corresponding causal function can be recovered via:

where σ is large enough that is defined for .

Important Remark

The functions in the time domain , with which we will deal, will always be

• causal functions,
• piecewise continuous functions,
• they contain no impulses at the time

Thus, the corresponding Laplace transforms will always be rational polynomial functions

Furthermore, the degree of the numerator polynomial will be less than that of the denominator polynomial: .

Laplace Transform: Zeros and Poles

Given a Laplace transform as a rational polynomial function, taking in evidence the roots of both the polynomials, we define

• zeros: the roots of the polynomial in the numerator

• poles: the roots of the polynomial in the denominator

Computation of the Inverse Laplace Transform "by Hand": Sum of Partial Fractions

Starting from a Laplace transform that is a strictly proper rational function, it is always possible to express the function as a sum of partial fractions:

where and .

How to calculate the coefficients?

An Example using partfrac

clear variables

 

syms t s

 

F(s) = (5*s+3)/(s^3+5*s^2+8*s+4)

F(s) = 

 

Fpart_frac = partfrac(F, s)

Fpart_frac(s) = 

The Inverse Laplace Transform using MATLAB: Examples

A First Example

determine the corresponding causal function (i.e., the inverse Laplace transform)

clear variables

 

syms t s

 

F(s) = (5*s+3)/(s^3+5*s^2+8*s+4)

F(s) = 

 

F(s) = simplify(F)

F(s) = 

f(t) = ilaplace(F, s, t)

f(t) = 

Important Remark on Notation

More Examples

determine the corresponding causal function (i.e., the inverse Laplace transform).

clear variables

 

syms t s

 

F(s) = (2/(s^2+2*s+2))*exp(-2*s)

F(s) = 

 

f(t) = ilaplace(F, s, t)

f(t) = 

determine the corresponding causal function

clear variables

 

syms t s

 

F(s) = (sym(2*sym(sqrt(3))+4)*s^3+sym(6*sym(sqrt(3))+20)*s^2+sym(14*sym(sqrt(3))+92)*s+sym(10*sym(sqrt(3))+108))/(s^4+4*s^3+26*s^2+44*s+85)

F(s) = 

 

f(t) = ilaplace(F, s, t)

f(t) = 

 

ft1 = simplify(ilaplace(F, s, t), 10)

ft1 = 

 

ft2 = simplify(ilaplace(F, s, t), 100)

ft2 = 

Zeros and Poles of a Rational Laplace Transform

Consider a Laplace transform that is a rational function of the complex variable s as, for example

clear variables

 

syms t s

 

F(s) = (5*s+3)/(s^3+5*s^2+8*s+4)

F(s) = 

To find the poles of simply sue the command poles:

pSET = poles(F, s)

pSET = 

[pSET, Np_SET] = poles(F, s) %

pSET = 

Np_SET = 

zSET = solve(F, s)

zSET = 

Important Remark on Zeros of a Rational Laplace Transform

G(s) = 3*((s+3)^3*(s-1))/((s+10)^2*(s-5)^3)

G(s) = 

 

[pSET, nSET] = poles(G, s)

pSET = 

nSET = 

 

zSET = solve(G, s, "MaxDegree",4)

zSET = 

Theorems on the Laplace Transform

• Can we determine the initial value of a time-domain function from its Laplace transform?
• Can we determine the steady-state value of a time-domain function from its Laplace transform?

Initial Value Theorem

Given a Laplace transform that is a strictly proper rational function, the following property holds

• You can analyse the time behavior of at directly from its Laplace transform , without computing the inverse Laplace transform!

An Application Example

Given the Laplace transform , determine without computing the inverse Laplace transform:

clear variables

syms s t

 

F(s) = (5*s+3)/(s+1)/(s+2)^2

F(s) = 

 

f0 = limit(s*F, s, +Inf)

f0 = 

0

f(t) = ilaplace(F,s, t)

f(t) = 

f(0)

ans = 

0

Final Value Theorem

Given a Laplace transform that is a rational function,

• if all the poles of are in the left open half-plane, with at most one simple pole in the origin, then the following identity holds and the two limits assume the same finite value.

Important Remark

Application Example

Given the Laplace transform , determine, if possible,

• determine the poles of ;
• verify that all the poles belong to the left open half-plane, with at most one simple pole at the origin.

[pF, n_pF] = poles(F, s) % poles and corresponding multiplicity

pF = 

n_pF = 

% --- 1st check: is there a pole at s=0? ---

pole0ID = isAlways(pF==sym(0));

if any(pole0ID)

pole0 = true;

p0M = isAlways(n_pF(pole0ID)==1);

else

pole0 = false;

p0M = false;

end

if (pole0 && p0M) || ~(pole0)

p0OK = true;

else

p0OK = false;

end

% --- 1st check: is there a pole at s=0? ---

 

 

% ---- 2nd check: Do all the remaining poles have negative real part? ----

if pole0

poles2analyse = not(pole0ID);

else

poles2analyse = true(size(pF));

end

pN_OK = all(isAlways(real(pF(poles2analyse))<sym(0)));

% final check

 

if (pN_OK && p0OK)

disp('The Final Value Theorem''s assumptions are met!');

disp('You can apply the Final Value Theorem!')

 

limf_infty = limit(s*F, s , 0)

 

% let's verify, by comparison with the inverse Laplace transform

f(t) = ilaplace(F, s, t);

f2infty = limit(f, t, +Inf)

 

else

disp('*** The Final Value Theorem''s assumptions are NOT met! ***');

disp('*** You CANNOT apply the Final Value Theorem! ***')

end

The Final Value Theorem's assumptions are met!

You can apply the Final Value Theorem!

limf_infty = 

f2infty =