MATLAB: Inverse Laplace problems – Matlab R2010a

Question

So I'm trying to get Matlab to find the time-response of a MDOF dynamic system (vibrations problem). I start by defining the transfer function in s-domain and then ask Matlab to find the inverse Laplace – it doesn't quite do it:

M=3*[1 0 0;0 1 0;0 0 1];
K=39*[2 -1 0;-1 2 -1;0 -1 2];
C=0.78*[2 -1 0;-1 2 -1;0 -1 2];
fs=[0;1;0];
syms s;
xs=(M*s^2+C*s+K)\fs;
x2=ilaplace(xs)

Matlab then returns (3×1 matrix):

x2 =
    sum((16250*exp(r5*t) + 325*r5*exp(r5*t))/(5000*r5^3 + 3900*r5^2 + 130338*r5 + 16900), r5 in RootOf(s5^4 + (26*s5^3)/25 + (65169*s5^2)/1250 + (338*s5)/25 + 338, s5))/3
    sum((32500*exp(r4*t) + 650*r4*exp(r4*t) + 1250*r4^2*exp(r4*t))/(5000*r4^3 + 3900*r4^2 + 130338*r4 + 16900), r4 in RootOf(s4^4 + (26*s4^3)/25 + (65169*s4^2)/1250 + (338*s4)/25 + 338, s4))/3
    sum((16250*exp(r3*t) + 325*r3*exp(r3*t))/(5000*r3^3 + 3900*r3^2 + 130338*r3 + 16900), r3 in RootOf(s3^4 + (26*s3^3)/25 + (65169*s3^2)/1250 + (338*s3)/25 + 338, s3))/3

I had a similar problem in PTC Mathcad 14 M020, which was solved by introducing a partial-fraction intermediate step. As far as I know, Matlab and Mathcad both use MuPAD, but I wasn't able to find symbolic partial fraction decompositions in Matlab.

Any ideas?

Answer 1

Hi Tamir,

you can apply a partial fraction via

   xs= feval(symengine, 'map', xs, 'partfrac', s)

This, however, does not help in this example, because the denominator of xs cannot be factorized.

If you do insist on avoiding symbolic sums over RootOfs in the ilaplace result, you can apply MuPAD's numeric::partfrac, which does a numerical factorization and a corresponding partial fraction expansion:

M=3*[1 0 0;0 1 0;0 0 1];
K=39*[2 -1 0;-1 2 -1;0 -1 2];
C=0.78*[2 -1 0;-1 2 -1;0 -1 2];
fs=[0;1;0];
syms s;
xs=(M*s^2+C*s+K)\fs
xs= feval(symengine, 'map', xs, 'partfrac', s)
x2=ilaplace(xs)
xs= feval(symengine, 'map', xs, 'numeric::partfrac', s)
x2=ilaplace(xs)
% beautification
xs= feval(symengine, 'map', xs, 'numeric::complexRound')
xs= feval(symengine, 'map', xs, 'numeric::rationalize')
x2=ilaplace(xs)
feval(symengine, 'float', x2)

Hope this helps,

Walter