# Find the transfer function of a differential equation symbolically

502 views (last 30 days)
Bill Tubbs on 6 Jun 2020
Edited: Bill Tubbs on 6 Jun 2020
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...

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)

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)
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?
##### 2 CommentsShowHide 1 older comment
madhan ravi on 6 Jun 2020