Clear Filters
Clear Filters

Plotting an Inverse Laplace Function

129 views (last 30 days)
Feliciana
Feliciana on 3 Apr 2024 at 17:34
Edited: David Goodmanson on 3 Apr 2024 at 21:19
I have an expression:
pretty(vsol)
(s + 100000) 147573952589676412928 14757395258967641292800000
---------------------------------- - --------------------------
#1 #1
2
(s + 100000 s + 10000000000) 2213609288845146 3 2
- ---------------------------------------------- + ((1106804644422573 s + 110680464442257300000 s
#1
+ 18446744073709550646400000 s + 737869762948382064640000000000) 6)/(s #1)
where
3 2
#1 == 5534023222112865 s + 700976274800962912928 s + 106991115627515394524800000 s
+ 5165088340638674452480000000000
Where when I do an inverse laplace transformation it turned into this:
pretty(ans)
/ 3 \
| --- |
| \ exp(t #2) |
| / --------- | 22136092888451460000000000
| --- #1 |
6 \ k = 1 /
- - ---------------------------------------------
7 7
-
/ 3 \
| --- |
| \ exp(#2 t) #2 |
| / ------------------------------------------------------------------------------- | 73786976294838187072
| --- 2 |
\ k = 11401952549601925825856 #2 + 16602069666338595 #2 + 106991115627515394524800000 /
-------------------------------------------------------------------------------------------------------------
7
/ 3 \
| --- 2 |
| \ exp(t #2) #2 |
| / ------------- | 2213609288845146
| --- #1 |
\ k = 1 /
- ---------------------------------------
7
where
2
#1 == 16602069666338595 #2 + 1401952549601925825856 #2 + 106991115627515394524800000
/ 2
| 3 700976274800962912928 z 21398223125503078904960000 z
#2 == root| z + ------------------------ + ----------------------------
\ 5534023222112865 1106804644422573
\
1033017668127734890496000000000 |
+ -------------------------------, z, k |
1106804644422573 /
How do I plot this? Where did the z come from?
  2 Comments
Torsten
Torsten on 3 Apr 2024 at 18:16
Your code to determine vsol ?
David Goodmanson
David Goodmanson on 3 Apr 2024 at 19:17
Edited: David Goodmanson on 3 Apr 2024 at 21:19
Hello Feliciana,
z is just a dummy variable to express a cubic polynomial, which has three roots. The index k denotes which root, those roots are used in the solution, at which point z is left behind.

Sign in to comment.

Answers (2)

Star Strider
Star Strider on 3 Apr 2024 at 19:47
If you have the Control System Toolbox —
% (s + 100000) 147573952589676412928 14757395258967641292800000
% ---------------------------------- - --------------------------
% #1 #1
% 2
% (s + 100000 s + 10000000000) 2213609288845146 3 2
% - ---------------------------------------------- + ((1106804644422573 s + 110680464442257300000 s
% #1
% + 18446744073709550646400000 s + 737869762948382064640000000000) 6)/(s #1)
% where
% 3 2
% #1 == 5534023222112865 s + 700976274800962912928 s + 106991115627515394524800000 s
% + 5165088340638674452480000000000
s = tf('s');
num(1) = ((s + 100000) * 147573952589676412928 - 14757395258967641292800000) - (14757395258967641292800000) - ((s^2 + 100000 * s + 10000000000) * 2213609288845146);
num(2) = ((1106804644422573 * s^3 + 110680464442257300000 * s^2 + 18446744073709550646400000 * s + 737869762948382064640000000000) * 6);
den(1) = 5534023222112865 * s^3 + 700976274800962912928 * s^2 + 106991115627515394524800000 * s + 5165088340638674452480000000000;
den(2) = s * den(1);
H = num(1)/den(1) + num(2)/den(2)
H = 2.45e31 s^6 + 6.37e36 s^5 + 1.296e42 s^4 + 1.622e47 s^3 + 1.405e52 s^2 + 8.548e56 s + 2.287e61 ---------------------------------------------------------------------------------------------------- 3.063e31 s^7 + 7.758e36 s^6 + 1.676e42 s^5 + 2.072e47 s^4 + 1.869e52 s^3 + 1.105e57 s^2 + 2.668e61 s Continuous-time transfer function.
figure
bp = bodeplot(H);
setoptions(bp, 'FreqUnits','Hz')
grid
figure
impulseplot(H)
figure
stepplot(H)
details = stepinfo(H)
Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.
details = struct with fields:
RiseTime: NaN TransientTime: NaN SettlingTime: NaN SettlingMin: NaN SettlingMax: NaN Overshoot: NaN Undershoot: NaN Peak: Inf PeakTime: Inf
This assumes that I correctly coded your results. It would be appropriate for you to check that to be certain.
.

Sam Chak
Sam Chak on 3 Apr 2024 at 20:16
Not exactly an answer, but rather an attempt to recreate the 'answer' that was displayed in your Command Window.
It appears that this is a 13th-order Transfer Function.
%% Visualizing the 13th-order Transfer Function
s = tf('s');
den = 5534023222112865*s^3 + 700976274800962912928*s^2 + 106991115627515394524800000*s + 5165088340638674452480000000000;
term1 = ((s + 100000)*147573952589676412928)/den;
term2 = (14757395258967641292800000)/den;
term3 = ((s^2 + 100000*s + 10000000000)*2213609288845146)/den;
term4 = ((1106804644422573*s^3 + 110680464442257300000*s^2 + 18446744073709550646400000*s + 737869762948382064640000000000)*6)/(s*den);
G = term1 - term2 - term3 + term4
G = 7.503e62 s^12 + 3.852e68 s^11 + 1.327e74 s^10 + 3.172e79 s^9 + 5.903e84 s^8 + 8.704e89 s^7 + 1.033e95 s^6 + 9.892e99 s^5 + 7.507e104 s^4 + 4.401e109 s^3 + 1.873e114 s^2 + 5.011e118 s + 6.1e122 ----------------------------------------------------------------------------------------------------------------------------------------------------------- 9.379e62 s^13 + 4.752e68 s^12 + 1.628e74 s^11 + 3.869e79 s^10 + 7.167e84 s^9 + 1.052e90 s^8 + 1.243e95 s^7 + 1.186e100 s^6 + 8.966e104 s^5 + 5.236e109 s^4 + 2.219e114 s^3 + 5.897e118 s^2 + 7.117e122 s Continuous-time transfer function.
%% Describe the transfer function symbolically
syms s
den = 5534023222112865*s^3 + 700976274800962912928*s^2 + 106991115627515394524800000*s + 5165088340638674452480000000000;
term1 = ((s + 100000)*147573952589676412928)/den;
term2 = (14757395258967641292800000)/den;
term3 = ((s^2 + 100000*s + 10000000000)*2213609288845146)/den;
term4 = ((1106804644422573*s^3 + 110680464442257300000*s^2 + 18446744073709550646400000*s + 737869762948382064640000000000)*6)/(s*den);
vsol = term1 - term2 - term3 + term4;
%% Convert the transfer function back to time-domain function
ft = ilaplace(vsol)
ft = 

Products

Community Treasure Hunt

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

Start Hunting!