142 views (last 30 days)
Robert on 2 Oct 2016
Edited: Walter Roberson on 18 Oct 2017
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???

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.

Robert on 2 Oct 2016
Thanks for your attempt at solving. Im somewhat disappointed to not be able to see my plot comparison
To suffice my curiosity can you confirm the following for me?
When i use the internal step function matlab is essentially taking the transfer function with an input of unit step , finding the laplace inverse and plotting it against time t correct?
If this is true, it must be doing my above problem a different method because it has no problem finding the step response of this function. Its only have a tuff time when I'm trying to do it manually myself
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.

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