Plotting an Inverse Laplace Function

4 Ansichten (letzte 30 Tage)
Feliciana
Feliciana am 3 Apr. 2024
Bearbeitet: David Goodmanson am 3 Apr. 2024
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 Kommentare
Torsten
Torsten am 3 Apr. 2024
Your code to determine vsol ?
David Goodmanson
David Goodmanson am 3 Apr. 2024
Bearbeitet: David Goodmanson am 3 Apr. 2024
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.

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Star Strider
Star Strider am 3 Apr. 2024
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 am 3 Apr. 2024
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 = 

Kategorien

Mehr zu Plot Customization finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by