MATLAB Answers

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

13 views (last 30 days)
Philip Mathew
Philip Mathew on 22 Jun 2021
Commented: Paul on 23 Jun 2021
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.

Accepted Answer

Paul
Paul on 23 Jun 2021
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)))
  2 Comments
Paul
Paul on 23 Jun 2021
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.

Sign in to comment.

More Answers (0)

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!