Strange ilaplace answer?????
22 views (last 30 days)
Show older comments
Obviously I'm doing something wrong here
Why am i getting such a strange answer for an inverse laplace??
syms s
tf = 10/(s^4 + 6*s^3 + 8*s^2 + 10*s)
>> ilaplace(tf)
ans =
1 - 8*symsum(exp(root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)*t)/(3*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)^2 + 12*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k) + 8), k, 1, 3) - 6*symsum((root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)*exp(root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)*t))/(12*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k) + 3*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)^2 + 8), k, 1, 3) - symsum((exp(root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)*t)*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)^2)/(3*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k)^2 + 12*root(s3^3 + 6*s3^2 + 8*s3 + 10, s3, k) + 8), k, 1, 3)
The answer is still in terms of S???
It did not even do the inverse? Its only a 4th order function? I don't get it
This TF happens to be the the step response of a system.
My goal was to transfer it back to the time domain so i could plot it in time.
I wanted to compare the manual plot i made to the internal step function of matlab
step(tf_closed_loop)
Obviously they should of been the same
but i can't get a reasonable answer from ilapalce???
0 Comments
Accepted Answer
Walter Roberson
on 2 Oct 2016
The output you are seeing is derived from the fact that the solution involves the sum of a polynomial times an exponential, with the sum taken over the three roots of a third degree polynomial. If you let the k'th root of the polynomial be designated by alpha(k) then the expression looks like
1 - 8*(symsum(exp(alpha(k)*t)/(3*alpha(k)^2+12*alpha(k)+8), k, 1, 3)) - 6*(symsum(alpha(k)*exp(alpha(k)*t)/(3*alpha(k)^2+12*alpha(k)+8), k, 1, 3))) - symsum(exp(alpha(k)*t)*alpha(k)^2/(3*alpha(k)^2+12*alpha(k)+8), k, 1, 3)
Maple's invplace transforms it to a different form:
1 + (1/611) * (symsum((53*alpha(k)^2+258*alpha(k)-41)*exp(alpha*t), k, 1, 3))
where again alpha(k) represents the k'th root of s3^3 + 6*s3^2 + 8*s3 + 10
MATLAB is not expanding the polynomial roots, because the expression gets out of hand. A simplified version of Maple's invlaplace is
(1/10998)*(((-(1656*I)*3^(1/2)+1656)*(135+3*1833^(1/2))^(1/3)+((414*I)*3^(1/2)+414)*1833^(1/2)+((53*I)*(135+3*1833^(1/2))^(4/3)+10998*I)*3^(1/2)-53*(135+3*1833^(1/2))^(4/3)-3666*(135+3*1833^(1/2))^(2/3)+10998)*exp(-(1/6)*((12*I)*3^(1/2)/(135+3*1833^(1/2))^(1/3)-I*3^(1/2)*(135+3*1833^(1/2))^(1/3)-12/(135+3*1833^(1/2))^(1/3)-(135+3*1833^(1/2))^(1/3)+12)*t)+(((1656*I)*3^(1/2)+1656)*(135+3*1833^(1/2))^(1/3)+(-(414*I)*3^(1/2)+414)*1833^(1/2)+(-(53*I)*(135+3*1833^(1/2))^(4/3)-10998*I)*3^(1/2)-53*(135+3*1833^(1/2))^(4/3)-3666*(135+3*1833^(1/2))^(2/3)+10998)*exp((1/6)*((12*I)*3^(1/2)/(135+3*1833^(1/2))^(1/3)-I*3^(1/2)*(135+3*1833^(1/2))^(1/3)+12/(135+3*1833^(1/2))^(1/3)+(135+3*1833^(1/2))^(1/3)-12)*t)+(106*(135+3*1833^(1/2))^(4/3)-3666*(135+3*1833^(1/2))^(2/3)-828*1833^(1/2)-3312*(135+3*1833^(1/2))^(1/3)-21996)*exp(-(1/3)*((135+3*1833^(1/2))^(1/3)+12/(135+3*1833^(1/2))^(1/3)+6)*t)+10998*(135+3*1833^(1/2))^(2/3))/(135+3*1833^(1/2))^(2/3)
where I is sqrt(-1).
The MATLAB generated expression, expanded, can be shown to be exactly the same.
At the moment I do not know how to get MATLAB to expand the root() into explicit values.
2 Comments
Walter Roberson
on 2 Oct 2016
I do not know how it works. Looking at the code I see that it invokes stepplot() and that stepplot() invokes ltiplot() which there appears to be little documentation for. It mentions that it invokes @resppack which seems to be a class. resppack.timeplot is invoked and I lose the thread after that.
If an inverse laplace is being done, then it could be happening numerically. The plots are not required to have indefinite precision like the symbolic ilaplace is.
More Answers (1)
Karan Gill
on 5 Oct 2016
Edited: Walter Roberson
on 18 Oct 2017
Use "vpa" to solve for the roots and get a result in terms of "t" as documented here: <http://www.mathworks.com/help/symbolic/root.html#buui6g4-3>
>> syms s
tf = 10/(s^4 + 6*s^3 + 8*s^2 + 10*s);
inverse = ilaplace(tf);
numInverse = vpa(inverse)
ans =
1.0 - 0.88866527104257189494371986264317*exp(-0.61959108298623234041484930760945*t)*cos(1.3101856107106995194523812766458*t) - 0.82480942549695286520323525995076*exp(-0.61959108298623234041484930760945*t)*sin(1.3101856107106995194523812766458*t) - 0.11133472895742810505628013735683*exp(-4.7608178340275353191703013847811*t)
Now, plot the result by using "fplot" (16a and after) or ezplot (pre-16a) as
fplot(numInverse)
1 Comment
Star Strider
on 5 Oct 2016
Using the vpa fiunction certainly never occurred to me in the Answer I offered (and then deleted when it wasn’t Accepted).
+1
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!