Manipulació simbòlica amb Matlab: Transformada de Laplace

Una utilització del Matlab diferent de la usual (numèrica) és la manipulació simbòlica, que consisteix en fer que el Matlab tracti les equacions i les variables de forma analítica, com ho faríem a mà. El Matlab ens permet calcular la transformada de Laplace d'una funció

de forma simbòlica (com una funció)

Manipulació simbòlica amb Matlab

Per utilitzar les variables sense que hagin de prendre un valor numèric, hem de dir-li al Matlab explícitament, es pot fer:

de una en una: x=sym('x')
varies juntes: syms x t s

clear

t=sym('t');

a=sym('a');

f = exp(-a*t) %les funcions que tenen variables simbòliques ho son també

f =

Com a funció, la podem derivar

derf = diff(f,t) %indiquem respecte a quina variable derivem

derf =

segderf = diff(f,t,2) %indiquem que derivem dues vegades

segderf =

També podem calcular la seva primitiva

intf=int(f)

intf =

Podem avaluar numèricament donant valors a les variables

subs(f,[a t],[1 0]) %es substitueix respectivament a=1, t=0

ans =

% Si volem tenir només una variable, podem fer a=ctant i

fa = subs(f,a,2) %ara és una funció d'una única variable

fa =

fplot(fa,[-3,3]) %dibuixa la funció en l'interval especificat

Transformada de Laplace (TL)

El Matlab ja te implementat el càlcul simbòlic de la transformada de Laplace d'una funció. La variable per defecte és la s.

f = exp(-a*t)

f =

F=laplace(f) %calculem la TL de la funció f anterior

F =

% qualsevol altra funció es defineix de forma anàloga (ha de ser simbòlica)

g=sin(a*t);

G=laplace(g)

G =

pretty(G) %mostrem les fraccions de la TL de forma més 'maca'.....

   a
-------
 2    2
a  + s

% fent la transformada inversa recuperem la funció original g(t):

ilaplace(G) %és la funció g

ans =

Funcions especials: La funció Delta de Dirac i la funció Heaviside

Dintre de les funcions especials per calcular les transformades de Laplace hi ha dues de molt importants:

funció delta de Dirac:

funció esglaó o funció de Heaviside

delta=dirac(t);

U=heaviside(t);

fplot(U)

axis([-5,5,-0.1,1.1])

diff(U) %la derivada és la delta de dirac

ans =

Calculem les seves transformades de Laplace

lapHeaviside=laplace(U)

lapHeaviside =

lapdDirac=laplace(dirac(t-3))

lapdDirac =

Solució d'EDO amb valors inicials

Com aplicació de la transformada de Laplace, veiem com s'aplica a resoldre equacions diferencials lineals amb coeficients constants.

Com un exemple, considerem la EDO de segon ordre :

syms t s y Y; %definim les variables simbóliques

% y serà la funció solució que busquem

% t és la variable independent de manera que y = y(t)

% Y és la seva Tranformada de Laplace (TL)

% s és també una variable independent de manera que Y = Y(s)

%la idea és aplicar la TL a tota la EDO per convertir-la en una equació

%algebràica de Y(s)

f=sin(2*t);

F=laplace(f,t,s);

y0 = -2 % valor inicial de y(t)

y0 = -2

dy0 = 3; % valor inicial de y'(t)

Y1=s*Y-y0; %transformada de Laplace de la derivada y', amb y(0)=-2

Y2=s^2*Y-s*y0-dy0; %transformada de Laplace de la derivada segona y'', amb y'(0)=3

Y = solve(Y2 + 2*Y1 + Y - F, Y) %calulem la transformada de tota l'equació i aïllem Y(s)

Y =

y = ilaplace(Y,s,t) %la solució és la antitransformada d=Y(s)

y =

comprovem que satisfà l'equació diferencial

diff(y,2)+2*diff(y,1)+y - sin(2*t)

ans =

fplot(y,[-1,0]) %dibuixem la funció solució

% Comprovem també les condicions inicials y(0)=-2; y'(0)=2

subs(y,t=0)

ans =

subs(diff(y,1), t=0)

ans =

Exercici: Amortidor de la roda d'un cotxe

Considerem el sistema format per l'amortidor d'un cotxe, que essencialment es pot modelar per una part elàtica (molla) i un amortidor (veieu figura). Es tracta d'una EDO de segon ordre on a partir de les forces exercides, es calculen els desplaçaments verticals de la carrosseria.

L'equació diferencial que simula aquest sistema és:

on m és la massa del vehicle, b el coeficient d'amorteiment i k el coeficient d'elasticitat de la molla i f és la força externa.

A partir dels valors: m = 1163, b = 1000, k = 10000 i amb x(0)=0, x'(0)=6

P1. Utilitzant la transformació i antitranformació de Laplace, calculeu el desplaçament vertical de l'amortidor

si hi simulem dos casos:

(1)

(2) f(x) com una força externa associada al perfil del sot que simulem amb la funció

h=10000;

f=h*(U(t-1) - U(t-2));

% per veure el tipus de funció

% fem el cas h=1

syms t;

f=heaviside(t-1)-heaviside(t-2);

fplot(f)

axis([0,3,-0.1,3])

hold off;

Visualitzeu les gràfiques de la solució amb la mateixa finestra per

% dades del problema

syms t s x X XX;

m = 1163; b = 1000; k = 10000; h=10000;

% Solucionem el primer cas f=0

f = 0;

F = laplace(f,t,s);

x0 = 0; % valor inicial de x(t)

dx0 = 6; % valor inicial de x'(t)

X1 = s*X-x0; %transformada de Laplace de la derivada x', amb x(0)=0

X2 = s^2*X-s*x0-dx0; %transformada de Laplace de la derivada segona x'', amb x'(0)=6

X = solve(m*X2 + b*X1 + k*X - F, X) %calulem la transformada de tota l'equació i aïllem X(s)

X =

x1 = ilaplace(X,s,t) %la solució és la antitransformada d=X(s)

x1 =

%ara fem el mateix pel segon cas f=h*(U(t-1)-U(t-2))

ff = h*(heaviside(t-1)-heaviside(t-2));

FF = laplace(ff,t,s);

x0 = 0; % valor inicial de x(t)

dx0 = 6; % valor inicial de x'(t)

XX1 = s*XX-x0; %transformada de Laplace de la derivada x', amb x(0)=0

XX2 = s^2*XX-s*x0-dx0; %transformada de Laplace de la derivada segona x'', amb x'(0)=6

XX = solve(m*XX2 + b*XX1 + k*XX - FF, XX) %calulem la transformada de tota l'equació i aïllem X(s)

XX =

x2 = ilaplace(XX,s,t) %la solució és la antitransformada d=X(s)

x2 =