Find the transfer function of a differential equation symbolically

138 views (last 30 days)
As an exercise, I wanted to verify the transfer function for the general solution of a second-order dynamic system with an input and initial conditions—symbolically.
I found a way to get the Laplace domain representation of the differential equation including initial conditions but it's a bit convoluted and maybe there is an easier way:
syms u(t) y(t) dydt(t) t s Y(s) U(s) y0 dydt0 u0 omega_n z K
% Standard form of second-order system
eqn_t = ((1/omega_n^2)*diff(y(t), t, 2) + (2*z/omega_n)*diff(y(t), t) + y) / K == u(t);
% In Laplace domain
eqn_s = subs(laplace(eqn_t), [laplace(y(t), t, s), laplace(u(t), t, s), diff(y(t), t)], [Y(s), U(s), dydt(t)])
% Set initial conditions to zero to get transfer function
eqn_s0 = subs(eqn_s, [y(0), dydt(0)], [0, 0])
This produces:
eqn_s =
-(dydt(0) + s*y(0) - omega_n^2*Y(s) - s^2*Y(s) + 2*omega_n*z*(y(0) - s*Y(s)))/(K*omega_n^2) == U(s)
eqn_s0 =
(omega_n^2*Y(s) + s^2*Y(s) + 2*omega_n*s*z*Y(s))/(K*omega_n^2) == U(s)
Which I think is the correct expression.
Next I want to rearrange this to get an expression for the transfer function (G(s)=Y(s)/U(s)). How do I do that? I tried this:
G(s) = Y(s)/U(s);
solve(eqn_s0,G(s))
But this produces:
ans =
struct with fields:
s: [0×1 sym]
z: [0×1 sym]
What does this mean? Both s and z are 'Empty sym: 0-by-1'.
What I want is:
G(s) =
(K*omega_n^2)/(omega_n^2 + 2*z*omega_n*s + s^2)
I'm sure there is a much easier way of doing all this...

Accepted Answer

Bill Tubbs
Bill Tubbs on 6 Jun 2020
Edited: Bill Tubbs on 6 Jun 2020
@madhan ravi was right. Using the simplify command works in this case:
rhs(eqn_s0) / U(s) / (simplify(lhs(eqn_s0)) / Y(s))
ans =
(K*omega_n^2)/(omega_n^2 + 2*z*omega_n*s + s^2)

More Answers (1)

Bill Tubbs
Bill Tubbs on 6 Jun 2020
Edited: Bill Tubbs on 6 Jun 2020
Maybe the problem is something to do with simplifying the left hand side when Y(s) is a function.
If you substitute Y(s) for a new symbolic variable and dividing by U(s) after solving it seems to work:
syms Ytemp
solve(subs(eqn_s0,Y(s),Ytemp),Ytemp)/U(s)
This produces:
ans =
(K*omega_n^2)/(omega_n^2 + 2*z*omega_n*s + s^2)
ADDITIONAL COMMENT:
Maybe this boils down to a more fundamental question. If you take the following expression, MATLAB doesn't simplify it:
syms a b
(a*b + b)/b
ans =
(b + a*b)/b
(2*b + b)/b
ans =
3
Maybe there a reason for that. Anyone know?

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!