MATLAB Answers

Inverse Laplace problems - Matlab R2010a

11 views (last 30 days)
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?
  1 Comment
RahulTandon
RahulTandon on 7 Jul 2015
use the 'partfrac' function. Now new in version 2015a of MATLAB!

Sign in to comment.

Accepted Answer

Walter Oevel
Walter Oevel on 10 Mar 2011
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
  2 Comments
john
john on 19 Mar 2014
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)

Sign in to comment.

More Answers (2)

Shahin
Shahin on 30 Apr 2012
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)
  1 Comment
john
john on 19 Mar 2014
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)

Sign in to comment.


Walter Oevel
Walter Oevel on 1 Apr 2014
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.
  2 Comments
john
john on 8 Apr 2014
I have no idea how to symbolical solve ilaplace

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!