# Is it possible to find the transfer function of these three differential equations using Mathematica

Imagine you are attempting to obtain a model of a device with two inputs [u1, u2] and three outputs [O1, O2, x] from the Lagrangian.

After some of Lagrangian work, you find the following three unsimplified equations:

E1 = -u1+ Kb r^2 (O1[t] - O2[t]) + Kb r (1/2 r (O1[t] - O2[t]) -
Cos[alpha]x[t]) + b1 O1'[t] + I1 O1''[t] == 0

E2 = -u2 - Kb r^2 (O1[t] - O2[t]) - Kb r (1/2 r (O1[t] - O2[t]) -
Cos[alpha]x[t]) + b2 O2'[t] + I2 O''2[t] == 0

E3 = Ks x[t] - 2 Kb Cos[alpha] (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) +
bs x'[t] + m x''[t] == 0

Kb, r, B1, I1, Ks, B2, I2, and m are known constants.
alpha is a nonlinear element, which, if needed, may be linearized around 30 degrees.


Regrerttably, probably due to misunderstanding the syntax, I cannot successfully apply LaplaceTransform on any of the three - but it is probably totally my fault, I'm very new to Mathematica.

What would you do next to isolate the systems transfer functions? Is there a way?

• Is that O''[t] in the def of E2 a typo for O1''[t] (or O2''[t])? Commented Dec 14, 2017 at 10:28
• Definitely a typo for O2''[t], thanks for pointing it out Commented Dec 14, 2017 at 11:16
• If the constants are known include their definitions in your code. The example you post should be self-contained. Commented Dec 14, 2017 at 11:58

I have had a look at your equations. Since u1 and u2 are inputs I assume they depend on time. I have thus changed your equations to

E1 = -u1[t] + Kb r^2 (O1[t] - O2[t]) +
Kb r (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) + b1 O1'[t] +
I1 O1''[t] == 0;

E2 = -u2[t] - Kb r^2 (O1[t] - O2[t]) -
Kb r (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) + b2 O2'[t] +
I2 O2''[t] == 0;

E3 = Ks x[t] -
2 Kb Cos[alpha] (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) +
bs x'[t] + m x''[t] == 0;


I assume you want transfer functions in symbolic form so we proceed with Laplace transforms and then solve for the Laplace transforms of O1, O2 and x;

e1 = LaplaceTransform[E1, t, s];
e2 = LaplaceTransform[E2, t, s];
e3 = LaplaceTransform[E3, t, s];
eqns = {LaplaceTransform[x[t], t, s], LaplaceTransform[O1[t], t, s],
LaplaceTransform[O2[t], t, s]} /.
First@Solve[{e1, e2, e3}, {LaplaceTransform[x[t], t, s],
LaplaceTransform[O1[t], t, s],
LaplaceTransform[O2[t], t, s]}] // Simplify;


You then have three rather complicated expressions for your three output variables. I am assuming that you don't want to look at the transient response but at the steady state response when the inputs are harmonic i.e. you wish to set s = I ω. You must therefore set the initial conditions to zero which may be done by using

ic = {O1[0] -> 0, O1'[0] -> 0, O2[0] -> 0, Derivative[1][O2][0] -> 0,
x[0] -> 0, x'[0] -> 0};


You have a matrix of transfer functions. If we set the input u2 to zero as well as the initial conditions we get

 h11 = 1/LaplaceTransform[u1[t], t, s] eqns[[1]] /.
Join[{LaplaceTransform[u2[t], t, s] -> 0}, ic] // Simplify


which gives

(2 Kb r (b2 + I2 s) Cos[ alpha])/((b1 + I1 s) ((3 Kb r^2 + 2 s (b2 + I2 s)) (Ks + s (bs + m s)) + 4 Kb (Kb r^2 + s (b2 + I2 s)) Cos[alpha]^2) + Kb r^2 (b2 + I2 s) (2 Kb + 3 (Ks + s (bs + m s)) + 2 Kb Cos[2 alpha]) )

This is one of your transfer functions. You could put in values and start plotting if you wished. You can find the others by setting u1 to zero and looking at the other equations. I also note that the denominator is fifth order in s so there are 5 poles. Not easy to work with unless you have values for your parameters.

Hope that helps.

• Hey Hugh, Thanks for the great solution! However, when I try to recreate it, I get the following tningy: (2 Kb r (b2 u1 + I2 s u1 - (b1 + I1 s) u2) Cos[alpha])/(s ((b1 + I1 s) ((3 Kb r^2 + 2 s (b2 + I2 s)) (Ks + s (bs + m s)) + 4 Kb (Kb r^2 + s (b2 + I2 s)) Cos[alpha]^2) + Kb r^2 (b2 + I2 s) (2 Kb + 3 (Ks + s (bs + m s)) + 2 Kb Cos[2 alpha])) LaplaceTransform[u1[t], t, s]) Any idea what the LaplaceTransform[u1[t], t, s] bit is about? Commented Dec 16, 2017 at 10:41
• Your differential equations have outputs and inputs. When you take the Laplace transform of the differential equations you get algebraic equations that contain the Laplace transforms of the outputs and inputs. These algebraic equations are linear equations and may be expressed in matrix form so that the vector of outputs equals a matrix times a vector of inputs. The matrix is the matrix of transfer functions. Thus the algebraic equations will have inputs like LaplaceTransform[u1[t],t,s] . The coefficients of these terms are the transfer functions. Does that make sense?
– Hugh
Commented Dec 18, 2017 at 22:26
• Hey Hugh, revisiting the project after the holidays, and it makes perfect sense. Thank you for your kindness and helpfulness! One last thing: how would you find out how many transfer functions are there?I am probably wrong, but if it's a 2 input 4 output MIMO system there should be either 2 or 4 transfer functions, and it seems like there are 3 Commented Jan 19, 2018 at 9:44
• Typically you have a matrix of transfer functions. So if there are 2 inputs and 4 outputs you have 8 transfer functions making up the matrix. Does that help?
– Hugh
Commented Jan 19, 2018 at 14:52
• Ouch. Yup, now to see what relates to what - thanks again Hugh! Commented Jan 19, 2018 at 16:46

It can be obtained by first converting it to a state-space representation.

E1 = -u1[t] + Kb r^2 (O1[t] - O2[t]) +
Kb r (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) + b1 O1'[t] +
I1 O1''[t] == 0;

E2 = -u2[t] - Kb r^2 (O1[t] - O2[t]) -
Kb r (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) + b2 O2'[t] +
I2 O2''[t] == 0;

E3 = Ks x[t] -
2 Kb Cos[alpha] (1/2 r (O1[t] - O2[t]) - Cos[alpha] x[t]) +
bs x'[t] + m x''[t] == 0;

tfm = TransferFunctionModel[
StateSpaceModel[{E1, E2, E3}, {x[t], O1[t], O2[t]}, {u1[t],
u2[t]}, {x[t], O1[t], O2[t]}, t], s]


The expressions are very involved. This is the term at {1, 1}.

tfm[s][[1, 1]]


• This is giving the same as my answer. The TransferFunctionModel saves you doing the LaplaceTransform which is useful.
– Hugh
Commented Dec 14, 2017 at 18:52
• This is less work. And there is an internal function that gives the result without having to compute a state-space representation which will be even faster. ControlDEqnstransferfunctionForm[{E1, E2, E3}, {{x[t], 0}, {O1[t], 0}, {O2[t], 0}}, {{u1[t], 0}, {u2[t], 0}}, {x[t], O1[t], O2[t]}, t, s] Commented Dec 14, 2017 at 20:02
• @SubaThomas I wish the documentation for state space in terms of the API could be improved. For example last API shown here StateSpaceModel.html is very incomprehensible. One must rely on looking at example to understand how to call many of these functions and using trial and error to figure how the calls work. I wish someone at WRI could take the time to write easier to understand API description. I find Matlab's API much easier to understand, even though I prefer the symbolic nature of Mathematica more. Commented Dec 14, 2017 at 20:12
• @Nasser, thanks for your feedback. I agree that some of the documentation for Control Systems can be improved. However, I believe that the usage messages of StateSpaceModel went through a couple of rounds of revisions and polishing. It's sometimes a balancing act between getting something that is both succinct (for the function templates, etc) and useful. May I ask what explicitly would you suggest as an improvement to that last usage message. Commented Dec 14, 2017 at 20:33
• Well, for example, the last API says StateSpaceModel[eqns, {{x1, x10}, \[Ellipsis]}, {{u1, u10}, \[Ellipsis]}, {g1, \[Ellipsis]}, \[Tau]] It is very hard to figure what to plugin what and where. Given the example for The state-space model of a set of ODEs below. only can guess that {{u[t],0}} comes from the input to the first ODE in the set. What is the 0 there? And what if the second ODE also has some input say u1[t]? Will it becomes {{u[t],0},{u1[t],0}} or {{u[t],u2[t]},{0,0}}` and so on. I mean the whole API is very cryptic. I start from examples, and reverse engineer them.. Commented Dec 14, 2017 at 20:46