Inverse Laplace problems - Matlab R2010a
3 views (last 30 days)
Show older comments
Tamir Duberstein
on 8 Mar 2011
Commented: RahulTandon
on 7 Jul 2015
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
Accepted Answer
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
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)
More Answers (2)
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
on 19 Mar 2014
Edited: 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)
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
on 2 Apr 2014
Edited: john
on 2 Apr 2014
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
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!