Is there a better method to plot the inverse Laplace of a function?

Hi everyone! I am facing an issue when plotting the inverse Laplace of a function. Here's what I am doing right now,

1. First I compute the inverse laplace of a function, say with a simple code.

syms s
U = 1/(s+1)
ut = ilaplace(U)

2. I then copy-paste the output into a new function and plot it.

syms t
t = 0:0.1:2
ut_plot = exp(-t)
plot(t, ut_plot)

Voila! I get the graph as expected. The reason I do this is because when I try plotting 'ut' directly it doesn't go through. But just for more context the function I am plotting is as follows.

ut_plot = 1000 - 134217728000*symsum((exp(t*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k))*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2)/(2*(201326592*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2 + 7829367466666667*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k) + 33554431999999995805696)), k, 1, 3) - 67108863999999991611392000*symsum(exp(t*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k))/(2*(201326592*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2 + 7829367466666667*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k) + 33554431999999995805696)), k, 1, 3) - 7829367466666667000*symsum((exp(root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)*t)*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k))/(2*(7829367466666667*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k) + 201326592*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2 + 33554431999999995805696)), k, 1, 3);

Hence it's more complicated and contains more roots.

Summarising my question: Is there a more efficient way for me to plot 'ut_plot' without having to copy-paste everytime? If not, is there a method to automate this copy-pasting? (Since I need to do it for at least a thousand values)

Thanks in advance for your answers! Please let me know if the question isn't clear, and apologies for any confusions.

Answer 1

[Paul]

The copy/paste approach shouldn't be needed. For many functions, fplot() usually does the trick:

syms s t
U(s)=1/(s+1);
u(t)=ilaplace(U(s));
figure;
fplot(u(t),[0 3])

If you want, you can turn u(t) into a function that can be evaluated (though there are restrictions, as in your actual case):

ufunc = matlabFunction(u(t))

ufunc = function_handle with value:
    @(t)exp(-t)

tvec = 0:.01:3;
figure
plot(tvec,ufunc(tvec))

For your particular case, vpa() seems to do the trick. Does the plot look like what you expect?

syms s6 k
ut_plot(t) = 1000 - 134217728000*symsum((exp(t*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k))*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2)/(2*(201326592*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2 + 7829367466666667*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k) + 33554431999999995805696)), k, 1, 3) - 67108863999999991611392000*symsum(exp(t*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k))/(2*(201326592*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2 + 7829367466666667*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k) + 33554431999999995805696)), k, 1, 3) - 7829367466666667000*symsum((exp(root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)*t)*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k))/(2*(7829367466666667*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k) + 201326592*root(s6^3 + (7829367466666667*s6^2)/134217728 + (7999999999999999*s6)/16 + 4166666666666666491904, s6, k)^2 + 33554431999999995805696)), k, 1, 3);
tvec = 0:1e-8:1e-6;
plot(tvec,vpa(ut_plot(tvec)))

[Philip Mathew]

This works perfectly! Just what I was looking for. Thank you so much, you saved my thesis :)

[Paul]

Depending on your needs you may be interested in converting the sym object U(s) into a tf object in the Control System Toolbox, and then using impulse() and other functions as needed. There should be several Answers that show how to do that conversion. Good luck.