More laplace inverse problems

1 view (last 30 days)
Robert
Robert on 9 Oct 2016
Commented: Robert on 9 Oct 2016
I can not get a distinct answer from the ilaplace function almost ever This is incredibly annoying.
This transfer function is not that complicated
There is no built in ramp function in matlab which means in order to plot it in matlab i need to either take the transfer function in s and multiple it by 1/s^2 ---> take the inverse of it and plot it
or
take the inverse of the transfer function and multiple that by t and plot that
Neither of which matlab can accomplish Not being able to plot a ramp response without simulink is getting extremely old Can someone help me with a solutions
num = [10 10];
dem = [1, 6, 10 ,18, 10];
sys = tf(num,dem)
tf_inverse = ilaplace(tf2sym(sys)) %just ignore the tf2sym functions its %just converting tf type to symbolic
ramp_response = t*(tf_inverse)
ezplot(ramp_response)
I just can not wrap my head around the fact that something as powerful as matlab can not take the inverse of a 4th order laplace

Accepted Answer

Walter Roberson
Walter Roberson on 9 Oct 2016
ezplot(vpa(ramp_response))
We went through this before. The inverse laplace of transfer function with maximum degree N in the denominator generally involves the sum of N roots of the polynomial that is the denominator. By default the symbolic engine leaves the root in symbolic form for degree 4 or higher.
The exact value of tf_inverse is
(1 / 189864) * (( - (30 * 1i) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + (10386 * 1i) * 3^(1 / 2) + ( - (363 * 1i) * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + (37752 * 1i) * 3^(1 / 2)) / (1774 + 54 * 1465^(1 / 2))^(1 / 3)) * (( - 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) + (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + (210 * ((1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) + 5082 * 3^(1 / 2) * ( - 2615 / 1694 + ( - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 1350 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 210 + (15 * (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 1560) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + ((121 * 1i) * (((((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 90 * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 104 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)))^(3 / 2) + 121 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(3 / 2)) * 3^(1 / 2)) * exp( - (1 / 6) * (1i * 3^(1 / 2) * (( - 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) + (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) - 3^(1 / 2) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + 9) * t) + (1 / 189864) * (((30 * 1i) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - (10386 * 1i) * 3^(1 / 2) + ((363 * 1i) * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) - (37752 * 1i) * 3^(1 / 2)) / (1774 + 54 * 1465^(1 / 2))^(1 / 3)) * (( - 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) + (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + (210 * ((1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) + 5082 * 3^(1 / 2) * ( - 2615 / 1694 + ( - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 1350 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 210 + (15 * (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 1560) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + ( - (121 * 1i) * (((((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 90 * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 104 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)))^(3 / 2) + 121 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(3 / 2)) * 3^(1 / 2)) * exp((1 / 6) * (1i * 3^(1 / 2) * (( - 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) + (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + 3^(1 / 2) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 9) * t) + (1 / 189864) * (( - 30 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 10386 * 3^(1 / 2) + (363 * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 37752 * 3^(1 / 2)) / (1774 + 54 * 1465^(1 / 2))^(1 / 3)) * ((14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + (210 * ((1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) - 5082 * 3^(1 / 2) * ( - 2615 / 1694 + ( - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + 1350 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 210 + (15 * (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 1560) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + (121 * (( - (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 90 * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + 104 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)))^(3 / 2) - 121 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(3 / 2)) * 3^(1 / 2)) * exp( - (1 / 6) * (3^(1 / 2) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - 3^(1 / 2) * ((14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + 9) * t) - (121 / 189864) * exp( - (1 / 6) * (3^(1 / 2) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + 3^(1 / 2) * ((14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + 9) * t) * (( - (30 / 121) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - (10386 / 121) * 3^(1 / 2) + (3 * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) - 312 * 3^(1 / 2)) / (1774 + 54 * 1465^(1 / 2))^(1 / 3)) * ((14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 90 * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2))^(1 / 2) + ( - (210 / 121) * ((1774 + 54 * 1465^(1 / 2))^(1 / 3) - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) + 42 * 3^(1 / 2) * ( - 2615 / 1694 + ( - (1 / 14) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 52 / 7) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) - (1350 / 121) * 3^(1 / 2) / (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + 210 / 121 + ( - (15 / 121) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 1560 / 121) / (1774 + 54 * 1465^(1 / 2))^(1 / 3) + ((( - (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(2 / 3) + 90 * 3^(1 / 2) * (1774 + 54 * 1465^(1 / 2))^(1 / 3) + 14 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2) + 104 * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)) / ((1774 + 54 * 1465^(1 / 2))^(1 / 3) * (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(1 / 2)))^(3 / 2) + (((1774 + 54 * 1465^(1 / 2))^(2 / 3) + 7 * (1774 + 54 * 1465^(1 / 2))^(1 / 3) - 104) / (1774 + 54 * 1465^(1 / 2))^(1 / 3))^(3 / 2)) * 3^(1 / 2))
Well that's pretty clear, right?
  3 Comments
Walter Roberson
Walter Roberson on 9 Oct 2016
vpa() is only going to give an approximation. That approximation is probably going to be close enough for graphics. If you were to laplace the result of the vpa of the ilaplace, you would get something that was not so nice. For example, if you vpa to 50 digits, you would get [2.10758034974676e+31 3.125e+64 3.125e+64] [3.125e+63 1.875e+64 3.125e+64 5.625e+64 3.125e+64] and with the default digits (32) you get a complex-valued system. The breakpoint appears to be vpa to 36 digits.
Maple's invlaplace expresses in terms of RootOf, and fortunately Maple's laplace is able to detect that RootOf and transform it back nicely.
Doing some checking, but not proving mathematically, it appears that laplace(ramp_response) you would get for the system (N1*s+N0)/Q(s) will be 1/Q(s)^2 * (N0*Q(s) + N1*L(Q(s)))
Here L is an arbitrary name for the transformation on the polynomial sum(a(i+1)*s^i) for i=0 to D (the maximum degree) to the polynomial sum((i-1)*a(i+1)*s^i) for i = 0 to D . For example,
a(5)*s^4 + a(4)*s^3 + a(3)*s^2 + a(2)*s^1 + a(1)*s^0
to become
(4-1)*a(5)*s^4 + (3-1)*a(4)*s^3 + (2-1)*a(3)*s^2 + (1-1)*a(2)*s^1 + (0-1)*a(1)*s^0
which is
3*a(5)*s^4 + 2*a(4)*s^3 + a(3)*s^2 - a(1)
This transformation can easily be done numerically to create the coefficients of the new numerator in transfer function form. There is probably a nice numeric way to square the denominator as well, but it is not immediately obvious to me.
Robert
Robert on 9 Oct 2016
After a little testing it appears that taking the inverse and multiplying it by t and plotting it is not giving me the correct answer
Taking the tf in S and multiplying that by 1/s^2 then taking the inverse and using the vpa function is giving the correct answer at least it appears so

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!