Inverse Laplace problems - Matlab R2010a

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 (3x1 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?

[RahulTandon]

use the 'partfrac' function. Now new in version 2015a of MATLAB!

Answer 1

[Walter Oevel]

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

[Tamir Duberstein]

Thanks Walter! Looks like the numerical factorization is what I needed.

[john]

I always get result in complex form. Why? for example:

xs=-(17000*(628000000*2^(1/2)*3^(1/2) + 2000000*2^(1/2)*s + 300*2^(1/2)*s^2 + 2^(1/2)*s^3 + 94200*2^(1/2)*3^(1/2)*s + 314*2^(1/2)*3^(1/2)*s^2))/((s^2 + 98596)*(7*s^3 + 19600*s^2 + 25200000*s + 20000000000))

xs= feval(symengine, 'map', xs, 'numeric::partfrac', s)
x2=ilaplace(xs)
feval(symengine, 'float', x2)

ans=2.200214881815407484640742403786/exp(1668.6981733835695764475415844609*t) + exp(314.0*t*i)*(2.2968284526268183887863094454911*i - 0.39524760709456622949419454901778) + exp(t*(- 1179.9312491214260645340516852799*i - 565.65091330821521177622920776954))*(0.37373997606394551823296182847589*i - 0.70485983381313751282617665287521) + exp(t*(1179.9312491214260645340516852799*i - 565.65091330821521177622920776954))*(- 0.37373997606394551823296182847589*i - 0.70485983381313751282617665287521) + 2.4041630560342615829628708311565*10^(-36)*dirac(t) + (1/exp(314.0*t*i))*(- 2.2968284526268183887863094454911*i - 0.39524760709456622949419454901778)

Answer 2

[Shahin]

you need to first simple the expression for xs and then convert it to double precision using vpa command as below.

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=simple(xs);

xs=vpa(xs,10); % ten digits precision

xt=ilaplace(xs)

xt =

 (4.9682060888811749440110359892427*10^(-37)*(cos(6.6473886206565217191026766304297*t) - 3.5684784608947646690759146132747*10^34*sin(6.6473886206565217191026766304297*t)))/exp(0.44384776310850235634421953414726*t) - (4.9682060888811749440110359892427*10^(-37)*(cos(2.758518538267630636964277414834*t) - 8.5992038062962428061290422100332*10^34*sin(2.758518538267630636964277414834*t)))/exp(0.076152236891497643655780465852739*t)
  (2.8583511709075181610190574710558*10^(-36)*(cos(2.758518538267630636964277414834*t) + 2.1137677058807934725481278701769*10^34*sin(2.758518538267630636964277414834*t)))/exp(0.076152236891497643655780465852739*t) - (2.8583511709075181610190574710558*10^(-36)*(cos(6.6473886206565217191026766304296*t) - 8.771666191057941295717388440763*10^33*sin(6.6473886206565217191026766304296*t)))/exp(0.44384776310850235634421953414726*t)
 (4.9682060888811749440110359892427*10^(-37)*(cos(6.6473886206565217191026766304297*t) - 3.5684784608947646690759146132747*10^34*sin(6.6473886206565217191026766304297*t)))/exp(0.44384776310850235634421953414726*t) - (4.9682060888811749440110359892427*10^(-37)*(cos(2.758518538267630636964277414834*t) - 8.5992038062962428061290422100332*10^34*sin(2.758518538267630636964277414834*t)))/exp(0.076152236891497643655780465852739*t)

[john]

simple, vpa doesn't help. What I'm doing wrong? I got result in complex form. Real command doesn't help :-(

>> ilaplace(vpa(simple(-(17000*2^(1/2)*s^3 + (5100000*2^(1/2) + 5338000*6^(1/2))*s^2 + (34000000000*2^(1/2) + 1601400000*6^(1/2))*s + 10676000000000*6^(1/2))/((s^2 + 98596)*(7*s^3 + 19600*s^2 + 25200000*s + 20000000000)))))

ans=2.2002148818154076690682970002036/exp(1668.6981733835695764475415844609*t) + exp(314.0*t*i)*(2.2968284526268184064188098725796*i - 0.39524760709456610146639250224196) + exp(t*(- 1179.9312491214260645340516852799*i - 565.65091330821521177622920776954))*(0.37373997606394535632565131956824*i - 0.70485983381313773306775599785985) + exp(t*(1179.9312491214260645340516852799*i - 565.65091330821521177622920776954))*(- 0.37373997606394535632565131956824*i - 0.70485983381313773306775599785985) + (1/exp(314.0*t*i))*(- 2.2968284526268184064188098725796*i - 0.39524760709456610146639250224196)

and also for above xt I got: xt =

   exp(t*(- 0.44384776310850235634421953414726 + 6.6473886206565217191026766304297*i))*(2.4841030444405874720055179946213e-37 + 0.0088644682087293467821068982264565*i) + exp(t*(- 0.076152236891497643655780465852739 + 2.758518538267630636964277414834*i))*(- 2.4841030444405874720055179946213e-37 - 0.021361308354985584586744545272149*i) + exp(t*(- 0.44384776310850235634421953414726 - 6.6473886206565217191026766304297*i))*(2.4841030444405874720055179946213e-37 - 0.0088644682087293467821068982264565*i) + exp(t*(- 0.076152236891497643655780465852739 - 2.758518538267630636964277414834*i))*(- 2.4841030444405874720055179946213e-37 + 0.021361308354985584586744545272149*i)
     exp(t*(- 0.44384776310850235634421953414726 + 6.6473886206565217191026766304297*i))*(- 1.4291755854537590805095287355279e-36 - 0.012536251164010178205593528271341*i) + exp(t*(- 0.076152236891497643655780465852739 + 2.758518538267630636964277414834*i))*(1.4291755854537590805095287355279e-36 - 0.030209451985654322420243569116655*i) + exp(t*(- 0.44384776310850235634421953414726 - 6.6473886206565217191026766304297*i))*(- 1.4291755854537590805095287355279e-36 + 0.012536251164010178205593528271341*i) + exp(t*(- 0.076152236891497643655780465852739 - 2.758518538267630636964277414834*i))*(1.4291755854537590805095287355279e-36 + 0.030209451985654322420243569116655*i)
   exp(t*(- 0.44384776310850235634421953414726 + 6.6473886206565217191026766304297*i))*(2.4841030444405874720055179946213e-37 + 0.0088644682087293467821068982264565*i) + exp(t*(- 0.076152236891497643655780465852739 + 2.758518538267630636964277414834*i))*(- 2.4841030444405874720055179946213e-37 - 0.021361308354985584586744545272149*i) + exp(t*(- 0.44384776310850235634421953414726 - 6.6473886206565217191026766304297*i))*(2.4841030444405874720055179946213e-37 - 0.0088644682087293467821068982264565*i) + exp(t*(- 0.076152236891497643655780465852739 - 2.758518538267630636964277414834*i))*(- 2.4841030444405874720055179946213e-37 + 0.021361308354985584586744545272149*i)

Answer 3

[Walter Oevel]

The complex entries are introduced by factorizing the denominator of xs numerically, which produces a product of terms with pairs of complex conjugate roots. The inverse Laplace transform consists of corresponding exp terms, involving these complex roots. In the overall combination of all these terms, however, the imaginary parts cancel (because to each entry there is a corresponding complex conjugate term).

You will see this by declaring

   syms t 'real'

and applying imag to xt (which produces zeros). So, you can apply xt = real(xt) and continue processing xt.

Hope this helps.

[john]

Hello Mr. Walter,

and what could I do with this? How can I solve it?

-((2^(1/2)*C*L2*U*sin((pi*phiU)/180)*s^3)/2 + ((2^(1/2)*C*R3*U*sin((pi*phiU)/180))/2 + (2^(1/2)*C*L2*U*omega*cos((pi*phiU)/180))/2)*s^2 + ((2^(1/2)*U*sin((pi*phiU)/180))/2 + (2^(1/2)*C*R3*U*omega*cos((pi*phiU)/180))/2)*s + (2^(1/2)*U*omega*cos((pi*phiU)/180))/2)/((omega^2 + s^2)*(R1 + R3 + L1*s - C*M12^2*s^3 + C*L1*L2*s^3 + C*L2*R1*s^2 + C*L1*R3*s^2 + C*L2*R3*s^2 + 2*C*M12*R3*s^2 + C*R1*R3*s))

This is analytical solving of electric circuit like before.

For numeric solution I use and it works:

s=ilaplace(result(i,1));
s=feval(symengine, 'float', s);
syms t 'real'
s=vpa(sqrt(2)*real(s),4);

Thank you

[john]

I have no idea how to symbolical solve ilaplace