Is there a better method to plot the inverse Laplace of a function?
조회 수: 51 (최근 30일)
이전 댓글 표시
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
채택된 답변
Paul
2021년 6월 23일
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))
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
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.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Calculus에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!