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

Question

Philip Mathew 2021-6-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/862590-is-there-a-better-method-to-plot-the-inverse-laplace-of-a-function

评论： Paul 2021-6-23

采纳的回答： Paul

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2021-6-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/862590-is-there-a-better-method-to-plot-the-inverse-laplace-of-a-function#answer_731060

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Philip Mathew 2021-6-23

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

Paul 2021-6-23

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.

请先登录，再进行评论。