Finding roots of differential equations

How can I find the roots/zeros of a differential equation such as

[Torsten]

What do you mean by "root of a differential equation" ?

Best wishes

Torsten.

[stevesy55]

@Torsten Solving the equation for s, roots being the x-intercepts of the differential function?

[Torsten]

What is s ? What is x ?

Best wishes

Torsten.

[John D'Errico]

Your question is fairly ambiguous, and can be interpreted in several ways.

Most logical to me would be to assume that since a differential equation implicitly defines a function y(t), then you are asking to solve for the roots of y(t), thus the set of values t such that y(t)=0.

So you need to explain what you are asking. As well, do the parameters B and C have known values? Are there known initial conditions on the differential equation? Is u(t) known?

[stevesy55]

Sorry for not explaining correctly, I'm new to differential equations, laplace and MATLAB.

The constant B = 7 and C = 0, there are no known initial conditions and u(t) is unknown. I need to find the roots, do a partial fraction expansion, an inverse Laplace transform and plot the function f(t). I then need to plot the poles and zeros of the differential equation.

I know some of the commands required are "roots", "residue", "ilaplace" and "plot".

Answer 1

[Star Strider]

If I understand your Question correctly, you have to transform your equation from the time domain into Laplace space (assume all initial conditions are 0), then solve for the roots in s-space. Any book on signal processing or differential equations will tell you how to do the transformation, and will have tables of Laplace transforms.

The transfer function, H(s), will then be given by:

H(s) = U(s)/Y(s)

The ‘roots’ of that function will be the zeros of the numerator. In this second-order system, this becomes a straightforward application of the quadratic formula. If you do not have numerical values for ‘B’ and ‘C’, you will necessarily get only a symbolic result.

[stevesy55]

Thank you very much for your reply Star Strider.

[Star Strider]

My pleasure.

You need to do everything you described in the Symbolic Math Toolbox. If you have R2015a and later, use the partfrac function to do the partial fraction expansion.

You can also isolate your transfer function coefficients and import them to the Signal Processing Toolbox (for discrete systems) or Control Systems Toolbox (continuous and discrete systems) to get the various plots, such as the pole-zero plots (the Control System Toolbox pzplot function). That is relatively straightforward to do. I will leave those details to you.

[stevesy55]

Thanks again Star Strider, much appreciated. May I just ask, is it better to use partfrac rather than residue? Also, once I've found the laplace of the differential equation in MATLAB should I then do the partial fraction expansion?

[Star Strider]

My pleasure.

The residue function will not work with symbolic objects.

Example:

n = sym([1 2]);
d = sym([3 4 5]);
[q,r] = residue(n,d);

throws the error:

Error using residue (line 97)
Inputs must be floats, namely single or double.

Error in ... (line ##)
[q,r] = residue(n,d);

So partfrac is the only way to solve this.

You first have to express your equation as a transfer function in Laplace s-space, then do the partial fraction expansion and the inversion.

[stevesy55]

I'm stuck on finding the laplace transform/transfer function of the differential equation. Can anyone help? This is what I've tried:

[Star Strider]

First, it is best to not use the single-quote syntax. That’s deprecated, and throws warnings. When you do that, you will also have to replace the single-equal ‘=’ with a double-equal ‘==’.

Then manually eliminate all the expressions with 0 arguments (such as ‘y(0)’ and the others), since all the initial conditions are 0 (by convention). (I have put in a request to have the initial conditions to the laplace function work similar to the way they’re entered in the dsolve function. That has not been implemented as of R2017a.)

Lastly, replace the ‘#1’ by ‘Y’ and ‘#2’ by ‘U’ (without the quotes, since I put them in just to make what I write more easily readable). You can use ‘Y(s)’ and ‘U(s)’ if you want. Remember to add these to your syms statement!

The Laplace and Inverse Laplace transform capabilities of the Symbolic Math Toolbox are mathematically correct, but inconvenient to work with. You have to do a bit of manual editing to make them work as you want them to.

That should get you past your current problem.

Your next step is to put your equation in the form of a transfer function. Then do the partial fraction expansion and the inversion.

[stevesy55]

@Star Strider When I leave out the single quotes I get an error:

[Star Strider]

Replace the single equal = with double equals ==.

That will solve the problem that threw that error. (I mentioned this in my previous Comment.)

[stevesy55]

@Star Strider I had tried the double equals aswell but was getting an error. I added symbolic variables u(t) and y(t) and it is now gone:

[Star Strider]

I feel your pain to the extent that I’m giving you my entire solution code:

syms s t U u(t) Y y(t) 
Dy = diff(y);
D2y = diff(y,2);
Du = diff(u);
D2u = diff(u,2);

Eqnt = D2y +6*Dy + 11*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)                                                    % Laplace Transform
Eqns = subs(Eqns, {'D(y)(0)', 'y(0)', 'D(u)(0)', 'u(0)'}, {0, 0, 0, 0})     % Zero Initial Conditions
H(s) = (6*s*laplace(y(t), t, s) + s^2*laplace(y(t), t, s) + 11*laplace(y(t), t, s))/(7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))
% H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {Y, U})
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1})
H(s) = partfrac(H, s)
h(t) = ilaplace(H)

figure(1)
fplot(h, [0 50])
grid

Go through it line by line to understand how it works.

I got frustrated with this, so I can imagine your situation. I wanted it to invert ‘U(s)’ to ‘u(t)’, and even simplification refused to transform the ‘ilaplace(U(s),s,t)’ to ‘u(t)’. The techniques I used here are those in the documentation on solving differential equations using Laplace transforms, so this seems to be the only option.

[stevesy55]

Thanks a million Star Strider! I'll go through your code and do my best to understand it.

[Star Strider]

My pleasure!

[stevesy55]

% declare symbolic variables needed by the program

syms s t U u(t) Y y(t)

% assign derivatives to variables

Dy = diff(y);
D2y = diff(y,2);
Du = diff(u);
D2u = diff(u,2)

Shouldn't there be a D3y = diff(y,3)?

[Star Strider]

There should.

The (Revised) Code —

syms s t U u(t) Y y(t) 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);

Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)                                                    % Laplace Transform
Eqns = subs(Eqns, {'D(D(y))(0)', 'D(y)(0)', 'y(0)', 'D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0})     % Zero Initial Conditions
H(s) = (11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s))/(7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))
% H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {Y, U})
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1})
H(s) = partfrac(H, s)
h(t) = ilaplace(H)

figure(1)
fplot(h, [0 50])
grid

[stevesy55]

Eqnt = D2y +6*Dy + 11*y == 7*D2u + 0*Du + u;                                % Time Domain Equation

Differential equation assigned to variable Eqnt

Eqns = laplace(Eqnt,t,s)                                                    % Laplace Transform

ℒ of differential equation Eqnt (t = independent variable, s = transformation variable) assigned to variable Eqns (Transfer Function)

Eqns = subs(Eqns, {'D(y)(0)', 'y(0)', 'D(u)(0)', 'u(0)'}, {0, 0, 0, 0})     % Zero Initial Conditions

Can you explain how the symbolic substitution works for setting zero initial conditions? Should there not be a D2u, D2y, D3y and 3 more zeros:

Eqns = subs(Eqns, {'D(D(D(y)))(0)','D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0, 0})

[Star Strider]

In my revised code, ‘D3y’ appears and is defined, so now:

Eqns =
11*s*laplace(y(t), t, s) - 11*y(0) - 6*D(y)(0) - D(D(y))(0) - 6*s*y(0) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) - s*D(y)(0) - s^2*y(0) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) - 7*s*u(0) - 7*D(u)(0) + laplace(u(t), t, s)

There are initial conditions ‘D(D(y))(0)’, ‘D(y)(0)’, and ‘y(0)’ on the LHS and ‘D(u)(0)’ and ‘u(0)’ on the RHS. The subs call sets all of them to 0. The result:

Eqns =
11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)

so they are replaced correctly.

The initial conditions terms derive from the integration in the Laplace transform. Differential equations textbooks all have a discussion of the Laplace transform. (If I remember correctly, they are constants-of-integration.) Please consult a reference to understand how the initial conditions derive from the transformation and are calculated.

[stevesy55]

Thanks Star Strider, so am I understanding this much correctly?

syms s t U u(t) Y y(t) 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);

Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)
Eqns = subs(Eqns, {'D(D(D(y)))(0)','D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0, 0})
H(s) = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) / 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)

[Star Strider]

My pleasure.

I believe so. It seems it’s possible to understand Laplace transforms and still not understand how to manipulate them efficiently in the Symbolic Math Toolbox!

[stevesy55]

@Star Strider This is the way I have the code now, can you please take a look if you don't mind and tell me where I'm going wrong? Not sure about the results I'm getting and the plot looks wrong. I'm entering the code exactly as it is below.

syms s t U u(t) Y y(t) 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);

Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)
Eqns = subs(Eqns, {'D(D(D(y)))(0)','D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0, 0})
% Eqns = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)
H(s) = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) / 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1})
H(s) = partfrac(H, s)
pretty(H)
h(t) = ilaplace(H)
pretty(h)
figure(1)
fplot(h, [0 50])
grid

[Star Strider]

The one problem I see is that this should be a ratio of two polynomials in ‘s’:

H(s) = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) / 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)

You need to enclose the numerator and denominator each in parentheses:

H(s) = (11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s)) / (7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))

Your code then produces the correct result.

[stevesy55]

@Star Strider Aha! Here are the results I'm now getting:

Trial>> syms s t U u(t) Y y(t) 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);

Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)
Eqns = subs(Eqns, {'D(D(D(y)))(0)','D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0, 0})
% Eqns = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)
H(s) = (11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s)) / (7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1})
H(s) = partfrac(H, s)
pretty(H)
h(t) = ilaplace(H)
pretty(h)
figure(1)
fplot(h, [0 50])
grid

Eqns =

11*s*laplace(y(t), t, s) - 11*y(0) - 6*D(y)(0) - D(D(y))(0) - 6*s*y(0) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) - s*D(y)(0) - s^2*y(0) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) - 7*s*u(0) - 7*D(u)(0) + laplace(u(t), t, s)

Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 
Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 
Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 
Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 
Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 
Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 
Warning: Support of character vectors that are not valid variable names or define a number will be removed in a
future release. To create symbolic expressions, first create symbolic variables and then use operations on them. 
> In sym>convertExpression (line 1559)
  In sym>convertChar (line 1464)
  In sym>tomupad (line 1216)
  In sym (line 179)
  In sym/subs>@(x)sym(x) (line 146)
  In sym/subs>normalize (line 146)
  In sym/subs>mupadsubs (line 137)
  In sym/subs (line 125) 

Eqns =

11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)

H(s) =

(11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s))/(7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))

H(s) =

(s^3 + 6*s^2 + 11*s + 3)/(7*s^2 + 1)

H(s) =

s/7 + ((76*s)/7 + 15/7)/(7*s^2 + 1) + 6/7

      76 s   15
      ---- + --
  s     7     7   6
  - + --------- + -
  7       2       7
       7 s  + 1

h(t) =

(76*cos((7^(1/2)*t)/7))/49 + (6*dirac(t))/7 + (15*7^(1/2)*sin((7^(1/2)*t)/7))/49 + dirac(1, t)/7

     / sqrt(7) t \                              / sqrt(7) t \
  cos| --------- | 76                sqrt(7) sin| --------- | 15
     \     7     /      6 dirac(t)              \     7     /      dirac'(t)
  ------------------- + ---------- + --------------------------- + ---------
           49                7                    49                   7

[Star Strider]

That looks correct to me. Those are the same results I got with my code.

Don’t worry about the ‘Warnings’. Using single quotes are the only ways it’s possible to substitute the initial conditions for the higher-order derivatives. Not using them throws an error about parentheses having to be last subscript reference. This is a problem with the way the Symbolic Math Toolbox interfaces with and returns results from MuPad.

The Symbolic Math Toolbox interface to MuPad is necessary for it to integrate with the rest of MATLAB. The problem is making it compatible with MATLAB syntax. I would not like to have a long symbolic expression of a Laplace transform and have to manually go through and replace the initial conditions vectors with zeros, since that will be the situation when the single-quote syntax is deprecated, unless the initial conditions syntax similar to that with dsolve is implemented first.

Fortunately, it’s possible to use 's'-domain expressions in the Control Systems Toolbox (and some others), so robust control system modeling and analysis of continuous systems are not dependent on the Symbolic Math Toolbox functions.

[stevesy55]

@Star Strider That's great, thank you so much you've really really saved me here.. If not for you I'd probably still be trying to find the laplace of the differential equation! So thank you, it's much appreciated. I've added some code to show the roots of the numerator of H(s):

Trial>> p = [1 6 11 3];
Trial>> r = roots(p)

r =

    -2.8358 + 1.0469i
    -2.8358 - 1.0469i
    -0.3283 + 0.0000i

Does this look ok to you?

syms s t U u(t) Y y(t) p r
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);
Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)
Eqns = subs(Eqns, {'D(D(D(y)))(0)','D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0, 0})
% Eqns = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)
H(s) = (11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s)) / (7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1})
p = [1 6 11 3];
r = roots(p) %roots of (s^3 + 6*s^2 + 11*s + 3)
H(s) = partfrac(H, s)
pretty(H)
h(t) = ilaplace(H)
pretty(h)
figure(1)
fplot(h, [0 50])
grid

[Star Strider]

That looks correct. The MATLAB convention is to have the coefficients of the highest-order term first, and the lowest last.

The Control System Toolbox gives the same results:

s = tf('s');
sys = (s^3 + 6*s^2 + 11*s + 3)/(7*s^2 + 1);
sysz = zero(sys)

If you were to plot them, note that the imaginary axis is the y-axis and the real axis is the x-axis. Plot the zeros with circles and the poles (roots of the denominator) with crosses.

[stevesy55]

@Star Strider Great stuff! Thanks for the advice, I need to plot the poles and zeros now so I'll check out the Control System Toolbox pzplot function that you suggested.

[stevesy55]

@Star Strider pzplot code added at the end:

Trial>> syms s t U u(t) Y y(t) 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);
Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u;                                % Time Domain Equation
Eqns = laplace(Eqnt,t,s)
Eqns = subs(Eqns, {'D(D(D(y)))(0)','D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0, 0})
% Eqns = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)
H(s) = (11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s)) / (7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s))
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1})
H(s) = partfrac(H, s)
pretty(H)
h(t) = ilaplace(H)
pretty(h)
figure(1)
fplot(h, [0 50])
grid
figure(2);
sys = tf([1 6 11 3],[7 0 1]);
P = pzplot(sys);
grid on;
z = getoptions(P);
z.Title.Color = [0 0 1];
setoptions(P,z);

[Star Strider]

You need to plot them in different figures:

s = tf('s');
sys = (s^3 + 6*s^2 + 11*s + 3)/(7*s^2 + 1);

figure(1)
fplot(h, [0 50])
grid

figure(2)
P = pzplot(sys);
set(gca, 'XLim',[-3 0.1])
axis equal
grid
z = getoptions(P);
z.Title.Color = [0 0 1];
setoptions(P,z)

I tweaked the pzplot call slightly.

[stevesy55]

Thanks a lot Star Strider I'll have a look.

[Star Strider]

My pleasure.

My tweaks just made the plot easier to read.

[stevesy55]

Yeah it looks great:

I'm not really sure what the tf command does to the plot, but the documentation says:

Use tf to create real- or complex-valued transfer function models (TF objects) or to convert state-space or zero-pole-gain models to transfer function form.

Just a quick question, does adding semi-colons at the end of certain commands mess things up i.e. suppress certain outputs?

[Star Strider]

I use the tf function when I have a system as a transfer function (as I do here, since all I needed to do was to copy the transfer function and paste it into my code). It’s an easy way to create a system object I can then explore.

Adding the semi-colons doesn’t mess things up, it speeds things up.

Using the semi-colons are usually best, because suppressing the Command Window output significantly speeds up code execution. I only leave them off with symbolic calculations, since I usually want to see those results. I put them back when I’m happy that my code is doing what I want it to.

I usually use them in non-symbolic code for the speed advantage, unless I want to see the details of a specific result, such as when I want to paste the result of a specific calculation to an Answer to demonstrate that my code works.

[stevesy55]

@Star Strider Thanks for explaining, I'll need to read up on system objects. So you made the transfer function into a system object by doing a tf on 's'? Are there any semi-colons I can include in the current code that won't suppress any calculation info in the command window output?

[Star Strider]

My pleasure.

The assignment:

s = tf('s');

creates a continuous-time transfer function object. The next line:

sys = (s^3 + 6*s^2 + 11*s + 3)/(7*s^2 + 1);

takes that a step further and defines the system object itself. It must be in s-domain transfer function form.

Your code does the same:

sys = tf([1 6 11 3],[7 0 1]);

but requires the extra step of creating vectors from the coefficients. (I’m lazy, especially if I have the transfer function already defined for me!)

You have to decide the results you want to display and those you want to suppress. I would suppress those that don’t need their results shown in the Command Window, and display those that do. (I’m not certain what the assignment wants you to show.)

[stevesy55]

I'm adding a title to the fplot. Should I use

f(t)= (76*cos((7^(1/2)*t)/7))/49 + (6*dirac(t))/7 + (15*7^(1/2)*sin((7^(1/2)*t)/7))/49 + dirac(1, t)/7

or

f(t) = y'''(t) + 6y''(t) + 11y'(t) + 3y(t) = 7u''(t) + u(t)

[Star Strider]

I would use the first one, since that’s the function you’re plotting.

If you want to be ambitious, you can use the TeX functions or even LaTeX: The Not So Short Introduction to LaTeX2e

[stevesy55]

Thanks again for explaining. So inverse laplace gives you the original function f(t) but not in derivative form? Sorry, I don't fully understand the math yet. I'll definitely have a look at LaTeX, thanks!

[Star Strider]

My pleasure.

The idea behind using the Laplace transform here is to use it to solve (or integrate) your differential equation. You are plotting the solution, so I would use the solution in the plot title.

You might be able to do everything you want in TeX. LaTeX offers a large number of additional formatting options. Use your pretty output as a formatting guide if you decide to use TeX or LaTeX. I mention those because they can make good plots into impressive plots. For example, instead of displaying dirac(t) you can actually display a lower-case ‘delta’ with \delta(t). See the documentation on the text function for details.

One other item is that you might want to use the vpa function to generate decimal fractions in ‘h(t)’ instead of all the fractions.

Also, since this is the time-domain expression of your transfer function, it is actually ‘h(t)’:

h(t)= 1.551*cos(0.37796*t) + 0.80992*sin(0.37796*t) + 0.85714*dirac(t) + 0.14286*dirac(1, t)

[stevesy55]

Bless you sir, thanks a mil.

[Star Strider]

Thank you. As always, my pleasure!

[stevesy55]

@Star Strider The final thing I need to do is to investigate the response of the transfer function. Plot the output via Simulink when the input to the transfer function is a unit step function. Any tips you could give would be great.

[Star Strider]

I haven’t done Simulink in a while, and when I tried it just now, I got a bunch of errors. It doesn’t like improper transfer functions, although the Control System Toolbox has no problem with them. (I may go back later to see if I can get the model running.)

The Control System Toolbox function for that is the step function.

[stevesy55]

Thanks Star Strider, I'm working on it now. I'll post whatever I come up with.

[stevesy55]

Final version of code and plots (Simulink to follow):

syms s t u(t) y(t) % Declare symbolic variables.
Dy = diff(y); % Assign derivatives to variables.
D2y = diff(y,2);
D3y = diff(y,3);
Du = diff(u);
D2u = diff(u,2);
Eqnt = D3y + 6*D2y + 11*Dy + 3*y == 7*D2u + 0*Du + u; % Define DE ("t-space" Time domain).
Eqns = laplace(Eqnt,t,s) % Laplace of DE gives transform of time domain variables("s-space" Frequency domain).
% Eqns = 11*s*laplace(y(t), t, s) - 11*y(0) - 6*D(y)(0) - D(D(y))(0) - 6*s*y(0) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) - s*D(y)(0) - s^2*y(0) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) - 7*s*u(0) - 7*D(u)(0) + laplace(u(t), t, s)
Eqns = subs(Eqns, {'D(D(y))(0)','D(y)(0)', 'y(0)', 'D(D(u))(0)','D(u)(0)', 'u(0)'}, {0, 0, 0, 0, 0, 0}); % Zero initial conditions
% Eqns = 11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s) == 7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)
H(s) = (11*s*laplace(y(t), t, s) + 6*s^2*laplace(y(t), t, s) + s^3*laplace(y(t), t, s) + 3*laplace(y(t), t, s)) / (7*s^2*laplace(u(t), t, s) + laplace(u(t), t, s)) % H(s) is transfer function. Divide '/'.
H(s) = subs(H, {laplace(y(t), t, s), laplace(u(t), t, s)}, {1, 1}) % Replace these instances with "1" for y and u.
H(s) = partfrac(H, s) % Partial fraction expansion of transfer function.
pretty(H); % Prettyprint to clean up answer.
h(t) = ilaplace(H) % Find time-domain expression of transfer function via ilaplace().
h = vpa(h,3) % Generate decimal fractions in ‘h(t)’ to 3 decimal places.
s = tf('s'); % Create continuous-time transfer function object.
sys = (s^3 + 6*s^2 + 11*s + 3)/(7*s^2 + 1); % Define system object in s-domain transfer function form.
figure(5); % Display figure number for plot.
fplot(h, [0 50], 'c'); % Create plot of h(t).
title(['\fontsize{10} \color[rgb]{0 0.5 0.5}h(t) = 1.55*cos(0.378*t) + 0.81*sin(0.378*t) + 0.857*\delta(t) + 0.143*\delta(1, t)']); % Create plot title.
xlabel('Time','fontweight','bold'); % Label axes.
ylabel('h(t)','fontweight','bold');
grid; % Grid toggle.
figure(6);
P = pzplot(sys); % Create pole-zero map of dynamic system model "sys". 
set(gca, 'XLim',[-3 0.1]); % Set limit for axes.
axis equal; % Set equal scale for x and y axes.
grid;
z = getoptions(P); % Set colour of plot title.
z.Title.Color = [0 0 1 1]; 
setoptions(P,z);

% Note: Some semi-colons have been ommitted to allow symbolic calculation
% results to be seen in the Command Window output i.e. the output is not suppressed.

[Star Strider]

My pleasure.

The problem I have when I ‘flip’ it so it becomes proper (in the Symbolic Math Toolbox) is that the inverse Laplace transform is complex in the time domain. Leaving it as an improper transfer function produces all real arguments and coefficients, and it appears stable, with its poles on the imaginary axis, which is legitimate.

Since complex coefficients are not appropriate in a time-domain function, the Symbolic and Control System realizations in the original solution (although ‘improper’) are correct. You may have to use something other than a ‘transfer function’ block in Simulink to get around Simulink not allowing improper transfer functions.

[stevesy55]

@Star Strider I'm just trying to absorb what you just said. I found this definition:

A system g(s) is improper if the order of the numerator polynomial exceeds the order of the denominator polynomial and proper otherwise.

[stevesy55]

and these: https://uk.mathworks.com/matlabcentral/answers/93570-how-can-i-implement-an-improper-transfer-function-without-delays-in-simulink

http://www.thebuchis.de/inversion/inversion_impossible.pdf

http://homepages.wmich.edu/~kamman/ME471ControlEffort.pdf

[Star Strider]

Thank you!

It didn’t occur to me to look there.

[stevesy55]

@Star Strider Please don't put yourself out to help me with this. I have just found out that there is a problem with some of the questions in this assignment (I'm guessing the improper TF is the problem) and so we've been given the option of alternative questions which I'm assuming give proper transfer functions:

[stevesy55]

@Star Strider for q1 above, should I include an initial condition for the 6 i.e. 6(0) or should it be left out?

[Star Strider]

This actually interests me!

That said, 2. is the only one that works, although 3. would be the easiest. The problem with 1. is that the inverse as I would have originally calculated it creates complex coefficients. Complex quantities simply do not exist in the time domain, at least in control theory, the reason I created it as an improper transfer function, since at least that is stable.

Do 2., since it’s the only one that’s stable when correctly calculated:

syms s t x(t) y(t) X Y 
Dy = diff(y);
D2y = diff(y,2);
Dx = diff(x);
Eqn = 9*D2y + 12*Dy + 13*y == 2*Dx + 3*x;
Eqns = laplace(Eqn, t, s);
Eqns = subs(Eqns, {'D(y)(0)', y(0), x(0), laplace(y(t), t, s), laplace(x(t), t, s)}, {0, 0, 0, Y, X});
% H = (Y*s^2 - 2*Y*s - 8*Y)/(4*U);
H = solve(Eqns, Y)/X;
H = partfrac(H, s);
h = ilaplace(H);
pretty(h)

figure(1)
fplot(h, [0 5])

[Hn, Hd] = numden(H);
Hnd = sym2poly(Hn);
Hdd = sym2poly(Hd);

sys = tf(Hnd, Hdd);
sysp = pole(sys);
sysz = zero(sys);

figure(2)
pzplot(sys)
grid
axis([-1.5  0.1    -1.1  1.1])
axis equal

figure(3)
step(sys)

Since 3. is unstable, it would likely not be a good option:

syms s t u(t) y(t) U Y 
Dy = diff(y);
D2y = diff(y,2);
Eqn = D2y - 2*Dy - 8*y == 4*u;
Eqns = laplace(Eqn, t, s);
Eqns = subs(Eqns, {'D(y)(0)', y(0), laplace(y(t), t, s), laplace(u(t), t, s)}, {0, 0, Y, U});
% H = (Y*s^2 - 2*Y*s - 8*Y)/(4*U);
H = solve(Eqns, Y)/U;
H = partfrac(H, s);
h = ilaplace(H);
pretty(h)

figure(1)
fplot(h, [0 5])

[Hn, Hd] = numden(H);
Hnd = sym2poly(Hn);
Hdd = sym2poly(Hd);

sys = tf(Hnd, Hdd);
sysp = pole(sys);
sysz = zero(sys);

figure(2)
pzplot(sys)

figure(3)
step(sys)

[stevesy55]

Thanks a lot for that Star Strider, I'm slightly lost and frustrated at this point as I thought I was nearly finished. I need to get 3 questions done by tomorrow evening.. Here are the final two options:

A = 4, B = 7, C = 0.

[Star Strider]

My pleasure.

With my code, this is just ‘plug & chug’. They’re both stable, but 2. is non-minimum-phase (a zero in the RHP).

Un-comment the one you want to solve and simulate:

syms s t x(t) y(t) X Y 
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
Dx = diff(x);
D2x = diff(x,2);
% Eqn = D3y + 4*D2y + 7*Dy + 18*y == 0*Dx + 3*x;              % 1.
% Eqn = 4*D2y + 7*Dy + 11*y == Dx - x;                        % 2.
Eqns = laplace(Eqn, t, s);
Eqns = subs(Eqns, {'D(D(y))(0)', 'D(y)(0)', y(0), 'D(x)(0)', x(0), laplace(y(t), t, s), laplace(x(t), t, s)}, {0, 0, 0, 0, 0, Y, X});
H = solve(Eqns, Y)/X;
H = partfrac(H, s);
h = ilaplace(H);
pretty(h)

figure(1)
fplot(h, [0 5])
grid

[Hn, Hd] = numden(H);
Hnd = sym2poly(Hn);
Hdd = sym2poly(Hd);

sys = tf(Hnd, Hdd);
sysp = pole(sys);
sysz = zero(sys);

figure(2)
pzplot(sys)
grid
axis([-4  0.1    -3  3])
axis equal

figure(3)
step(sys)
grid

figure(4)
impulse(sys)
grid

[stevesy55]

@Star Strider Why did you not set all the derivatives and functions to zero initial conditions in q2 and q3? Also, the use of Y, U and X with the zeroes confuses me.

Eqns = subs(Eqns, {'D(y)(0)', y(0), laplace(y(t), t, s), laplace(u(t), t, s)}, {0, 0, Y, U});
% H = (Y*s^2 - 2*Y*s - 8*Y)/(4*U);
H = solve(Eqns, Y)/U;

[Star Strider]

I did set all the initial conditions to zero with all of them. Only the initial conditions were set to zero. I also remember checking them to be sure the substitutions were correct. I did not set the Laplace transforms of the variables to zero, I assigned them as their upper-case single-letter counterparts instead. This is to get around some of the problems in the Symbolic Math Toolbox laplace and ilaplace functions.

The upper-case letters conventionally denote the Laplace transform of their lower-case time domain variables in anything involving Laplace transforms. So ‘Y = Y(s) = L{y(t}}’, and so for the others.

In the subs call, I substituted everything at once for efficiency.

[stevesy55]

I'm still trying to understand the initial conditons part of the code and solve. The transfer functions have a higher order numerator so they are improper? Will they work in Simulink? Thanks for your time.

[stevesy55]

For example:

In q2.

Eqns = subs(Eqns, {'D(y)(0)', y(0), x(0), laplace(y(t), t, s), laplace(x(t), t, s)}, {0, 0, 0, Y, X});

I don't understand why laplace(y(t), t, s) and laplace(x(t) are there, and why D(D(y))(0) and D(x)(0) are not there.

Similarly in q3.

Eqn = D2y - 2*Dy - 8*y == 4*u;
Eqns = laplace(Eqn, t, s);
Eqns = subs(Eqns, {'D(y)(0)', y(0), laplace(y(t), t, s), laplace(u(t), t, s)}, {0, 0, Y, U});

I'm wondering where D(D(y))(0) and u(0) are.

[Star Strider]

My pleasure.

The solve call creates the transfer function. Since the transfer function is defined as ‘Output(s)/Input(s)’, solving for ‘Y’ as the output (since it’s common to all of these), then dividing by either ‘U’ or ‘X’ as the input (depending on the system) creates the transfer function. (Note that the upper-case letters denote the Laplace transforms of the time-domain variables.)

The last two are stable. You should have no problems getting them to work in Simulnk. Use the ‘Hnd’ and ‘Hdd’ vectors as the numerator and denominator functions, respectively, as I did in the tf call.

I had problems with the first one (the problem you originally posted) since it gave such a weird solution. (I couldn’t imagine a homework problem that was not reasonably straightforward.) The others are well-behaved, even if not all are stable or minimum-phase.

[stevesy55]

I thought laplace created the transfer function, as in the code for the original problem I posted.

Ok I will focus on the last two and the other stable q2 and leave the original problem? Sorry for the confusion, it's just not clicking at the moment.

[stevesy55]

I think something clicked with me there. In the old code the laplace transform gave the transforms of the y and u time-domain variables on either side of the '==', so this wasn't the TF until we replaced '==' with '/' which would be Y(s)/U(s)?

[Star Strider]

q2: ‘I don't understand why laplace(y(t), t, s) and laplace(x(t) are there, and why D(D(y))(0) and D(x)(0) are not there.’

The laplace variables are there because I want to replace them with ‘Y’ and ‘X’ respectively, to make the inversion easier.

Since:

Eqn = 9*D2y + 12*Dy + 13*y == 2*Dx + 3*x;

has one second derivative on the left and one first derivative on the right, there are no higher derivatives with initial conditions.

q3: ‘I'm wondering where D(D(y))(0) and u(0) are.’

The highest derivative is a second derivative on the left, so only ‘D(y)(0)’ (the initial condition for the first derivative) and ‘y(0)’ appear, and on the right, ‘u(t)’ simply appears in the transformed equation as ‘U(s)’, or more simply in my code, ‘U’.

I’ll have to think about this one:

Eqn = D3y + 5*D2y + 17*Dy + 13*y == D2x + 5*Dx + 6;

You can’t ignore the 6, since it transforms to a ‘6/s’ step. That system may not have a solution.

[Star Strider]

‘so this wasn't the TF until we replaced '==' with '/' which would be Y(s)/U(s)?’

Correct!

However, we didn’t replace ‘==’ with ‘/’. We solved for ‘Y’ then divided that solution by ‘U’ to create ‘Y/U’ that we assigned to ‘H’ (or ‘H(s)’) as the transfer function.

[stevesy55]

So laplace solved for Y and U (transformed all the y and u time-domain variables) and we then divided? Can you explain exactly what the subs command is doing, and maybe how the initial conditions work? From what I understood, if there are no initial conditions, you put a zero in for each function/derivative in the differential equation via the subs command?

[stevesy55]

Looking at the last 2 questions I can see in the command window output that the Y and X in the subs command replace all the instances of laplace(y(t), t, s) and laplace(x(t), t, s) in the transform result giving the transform in this form:

Eqns = 4*Y*s^2 + 7*Y*s + 11*Y == X*s - X

I see it is then inverted and the Ys are gone:

H =

-(X - X*s)/(X*(4*s^2 + 7*s + 11))

And after partial fraction expansion:

H =

    (s - 1)/(4*s^2 + 7*s + 11)

For this one the there are no Xs or Ys left in the equation after inversion:

Eqns =

    Y*s^3 + 4*Y*s^2 + 7*Y*s + 18*Y == 3*X

      H =

      3/(s^3 + 4*s^2 + 7*s + 18)

      H =

      3/(s^3 + 4*s^2 + 7*s + 18)

[stevesy55]

I think I get the what's happening with the subs command, you're just getting rid of variables you don't want from the tranform result. So you put a 0 in for them, or replace them with Y and X for example?

[Star Strider]

‘So laplace solved for Y and U (transformed all the y and u time-domain variables) and we then divided?’

Yes.

‘Can you explain exactly what the subs command is doing, and maybe how the initial conditions work?’

The subs function does exactly what it says. It takes the first argument expression, looks up the variables in the second argument, and replaces them with the values in the third argument.

I refer you to a textbook on differential equations for a discussion of initial conditions, and to a discussion of Laplace transforms with respect to how they incorporate initial conditions. That is much too extensive a discussion to get into here.

‘I think I get the what's happening with the subs command, you're just getting rid of variables you don't want from the tranform result. So you put a 0 in for them, or replace them with Y and X for example?’

I’m setting the initial conditions to zero, because in control systems, by convention, all initial conditions are zero. It’s not that I don’t ‘want’ them, it’s that I want them to be zero.

I replace the laplace() expressions with simpler expressions, since that makes the subsequent inversion easier.

You’re correct about eliminating the variables in the transfer function calculation, since that’s how transfer functions are defined. Transfer functions must be expressions of polynomials only in s.

[stevesy55]

@Star Strider

The subs function does exactly what it says. It takes the first argument expression, looks up the variables in the second argument, and replaces them with the values in the third argument.

I’m setting the initial conditions to zero, because in control systems, by convention, all initial conditions are zero. It’s not that I don’t ‘want’ them, it’s that I want them to be zero.

I replace the laplace() expressions with simpler expressions, since that makes the subsequent inversion easier.

I'm clearer about this now, thank you. I'm not sure how you performed the inversion via solve. Regarding the intial conditions could you say that the laplace transform result provides you with all possible initial conditions and you include or omit them as is appropriate? Also, you converted improper transfer functions into proper transfer functions but they are now in the form of U(s)/Y(s) when going into Simulink? Many thanks again for your help Star Strider I owe you.

[Star Strider]

My pleasure.

The solve call doesn’t do the inversion. It solves for the output (here ‘Y’) and then after solving, divides by the input (here ‘U’ or ‘X’) to create the ratio of polynomials in s only that define the transfer function. Then I do the partial fraction decomposition with partfrac, and after that, the inversion.

‘Regarding the initial conditions could you say that the Laplace transform result provides you with all possible initial conditions and you include or omit them as is appropriate?’

Yes. Use the subs function to set them to be what you want. In control systems, they are by convention all zero.

‘Also, you converted improper transfer functions into proper transfer functions but they are now in the form of U(s)/Y(s) when going into Simulink?’

No. I did that for the first one only because I thought the definition was in error (the y and u terms were reversed), since the inversion of the function in the original form produced a complex result in the time domain. I cannot imagine that a complex result in the time domain would ever be appropriate.

The transfer function defined as Y(s)/U(s) is correct for all the others. They all appear proper (although not all are stable or minimum-phase), and you should have no problems using a transfer function block for them.

‘Star Strider I owe you.’

I appreciate your sentiment, and I thank you. You don’t owe me, actually. The best way you can reciprocate is to share what you learned in this thread with others. (Referring them to it and asking them to vote for it if they consider it helpful to them would be nice!)

[stevesy55]

Sorry Star Strider by inversion I meant the way you invert the numerator and denominator so as the transfer function is proper (numerator has a lower degree than denominator).

e.g.

Eqns =

Y*s^3 + 4*Y*s^2 + 7*Y*s + 18*Y == 3*X

H =

3/(s^3 + 4*s^2 + 7*s + 18)

No. I did that for the first one only because I thought the definition was in error (the y and u terms were reversed), since the inversion of the function in the original form produced a complex result in the time domain. I cannot imagine that a complex result in the time domain would ever be appropriate.

Which question was that again?

[Star Strider]

No worries.

That is inherent in the solution. Breaking it into two steps:

Eqns = Y*s^3 + 4*Y*s^2 + 7*Y*s + 18*Y == 3*X;
Ys = solve(Eqns, Y)
H = Ys/X

produces:

Ys =
(3*X)/(s^3 + 4*s^2 + 7*s + 18)

H =
3/(s^3 + 4*s^2 + 7*s + 18)

The algebra is straightforward. First, factor ‘Y’ from the LHS, then to solve for ‘Y’ on the LHS, divide both sides of the equation by ‘(s^3 + 4*s^2 + 7*s + 18)’. This produces:

Y*(s^3 + 4*s^2 + 7*s + 18)/(s^3 + 4*s^2 + 7*s + 18) == 3*X/(s^3 + 4*s^2 + 7*s + 18)

Y == 3*X / (s^3 + 4*s^2 + 7*s + 18)

To calculate the transfer function, divide through by ‘X’ to get:

H == 3 / (s^3 + 4*s^2 + 7*s + 18)

the result we want.

[stevesy55]

So even though the numerator is now X and denomiator is Y, the transfer function H(s) = Y(s)/X(s)?

[stevesy55]

So solve, aside from solving for Y, puts X over Y? Which is how the 3X here on the RHS

Eqns = Y*s^3 + 4*Y*s^2 + 7*Y*s + 18*Y == 3*X;

ends up on the LHS here

Ys =
(3*X)/(s^3 + 4*s^2 + 7*s + 18)

[Star Strider]

‘So even though the numerator is now X and denomiator is Y, the transfer function H(s) = Y(s)/X(s)?’

Yes.

‘So solve, aside from solving for Y, puts X over Y?’

No.

It puts the coefficient polynomial of ‘X’ over the coefficient ploynomial of ‘Y’. The transfer function is defined as ‘H=Y/X’.

The solve call first solves for the output, ‘Y’. It then divides the solution for ‘Y’ by the input, ‘X’, to create the transfer function, ‘Y/X’. The transfer function are the coefficient polynomial of ‘X’ divided by the coefficient polynomial of ‘Y’ to create the transfer function, ‘H = Y/X’.

I realise this must seem incredibly confusing if you’re seeing it for the first time. However the maths are straightforward and make sense when you analyse them. It helps to do this on paper with simple transfer functions so you can see how it all works. The transfer function calculations are actually straightforward beginning algebra, when you solve for a specific variable (in control systems this is the output ‘Y’) in terms of the input variable (here ‘X’ or ‘U’). It gets confusing when you then divide ‘Y’ by the input variable to create the transfer function as purely a ratio of polynomials.

Don’t worry. You’ll understand it, and it won’t take long before you’re comfortable with all the concepts. To the best of my knowledge, no one was ever born knowing it, so everyone who’s ever encountered it will have gone through the same process to understand it.

[stevesy55]

@Star Strider Thanks for the advice and the kind words, I'll stick at it. One more question..

' It puts the coefficient polynomial of ‘X’ over the coefficient ploynomial of ‘Y’. The transfer function is defined as ‘H=Y/X’.'

Do you mean whatever variables/values are attached to 'X' will be placed over the variables/values attached to'Y'?

[Star Strider]

My pleasure.

‘Do you mean whatever variables/values are attached to 'X' will be placed over the variables/values attached to 'Y'?’

Yes. That’s just the way the math works in the transfer function calculation.

[stevesy55]

That makes me happy. I can sleep now.

[Star Strider]

Me, too!

[stevesy55]

May I just say a big thank you again Star Strider, I'm a newbie but already have a very high opinion of this forum. Thanks!

[Star Strider]

As always, my pleasure!

I usually use the Control System Toolbox and the numerical differential equation solvers, not Laplace transforms, for these problems, so it was good to review them. It was something of a challenge for me to Answer your Questions and explain my answers since I’ve not revisited the details of this topic in some time. This has been quite definitely beneficial for both of us, and I hope for everyone else who reads this thread and learns from it.

I wish you well in your studies!

[MINATI]

Can it be possible to solve given eqn by using LAPLACE TRANSFORM with two independent (x, t) variables.

[Vignesh Ramakrishnan]

If this equation were nonlinear: say:

yddot=a*sin(y)+b*cos(y)*dy^2;

How do we solve for the "roots" of the solution obtained numerically, if we parameterize a,b in order to roughly get the non-linear equivalent of a root locus?

[Star Strider]

@MINATI — You may be describing a partial differential equation. One way to approach that is to take the Laplace transform with respect to time ‘t’, and then solve the transformed equation as an ordinary differntial equaation in ‘x’, and then invert it. I did that many years ago and got it to work, however I would be reluctant to replicate that. (It is not as straightforward as it might first appear by my description of it.)

@Vignesh Ramkrishnan — I believe that nonlinear differential equations cannot be transformed into Laplace space. That is limited to ordinary differetial equations with constant coefficients, although it is occasionally possible to use Laplace transforms to solve relatively simple partial differential equations (see my reply to MINATI above).

[Vignesh Ramakrishnan]

@Star Strider, I am not talking of using the Laplace domain(which is restricted to linear and time-invariant systems) here. Say I have the above non-linear differential system(which I am using to model a pendulum with some weird damping), for which, say a,b are chosen and fixed for obtaining the solution's equilibria by checking the zero crossing when I run this through ode45, and place all this in a loop where a,b are varied, that should be a rough equivalent of a 'root locus'. Something like

for a=0:0.1:10
    for b=0:0.1:10
        y_roots=roots(ddy(t,a,b));
    end
end

Is there an easier way to pull this off?

[Star Strider]

Probably the easiest way to detect zero crossings is to use an Events location and function. That should automatically record the times.

I leave the rest of that to you.