Plotting a ratio $G(s) = \frac{Y(s)}{X(s)}$ with an impulse input on MATLAB

Question

I am trying to simulate a system as such.

$G(s) = \frac{Y(s)}{X(s)}$

such that $G(s) = \frac{s^3 + 4s^2 +6s+8}{s^3 + 3s^2 +5s}$

With an impulse input of $sin(t)$

I am trying to do this in Matlab and was going about it in the following way

In the code I get up to LINE 91 before an error occurs as it can't compute a Laplace of the equation. I am unaware how to fix this and get past my errors. Is there an easier way of doing this?

I have also provided what I am trying to answer as further clarification towards the problem.

$ \dddot y + 3\ddot y +5\dot y = \dddot x+4\ddot x +6\dot x +8x $

Find $G(s) = \frac{Y(s)}{X(s)}$, when all initial conditions are zero.

when $x(t)=\delta(t)$, what is $y(t)$ as t $\implies \infty$

Plot the simulation (I am using Matlab).

Any help would be appreciated as I am now confusing myself in my own working.

Link to my other question involving the same problem: Laplace Transform: Find Ratio $G(s)$ and value of $y(t)$ as $t\to \infty$

Answer 1

You can find a staggered first-order system for this equation requiring no derivative of the input $x$ as \begin{align} y&=u_0+x,& \implies u_0'''+3u_0''+5u_0'&=x''+x'+8x, \\ u_0'&=u_1+x,& \implies u_1''+3u_1'+5u_1&=-2x'+3x, \\ u_1'&=u_2-2x,& \implies u_2'+3u_2+5u_1&=9x \\ u_2'&=-3u_2-5u_1+9x. \end{align}

This now can be symbolically or numerically solved.

For $x=\delta$ the system simplifies for $t>0$ to \begin{align} y&=u_0,& y_0(0)&\text{ does not exist} \\ u_0'&=u_1,& u_0(0^+)&=1 \\ u_1'&=u_2,& u_1(0^+)&=-2 \\ u_2'&=-3u_2-5u_1, &u_2(0^+)&=9. \end{align}

For $x(t)=\theta(t)\sin(t)$, $\theta(t)=0$ for $t<0$ and $=1$ for $t>0$, the right sides of the system are regular, so a differentiable solution exists for $u_0$, $y$ is not differentiable at $t=0$. This means that the first order system has initial conditions equal to zero.

def x(t): return sin(maximum(0,t))
def u_fun(u,t): xt=x(t); return [ u[1]+xt, u[2]-2*xt, -3*u[2]-5*u[1]+9*xt]
t=linspace(-1,20,301)
u = odeint(u_fun, [0,0,0], t)
plot(t,u[:,0],':r',label='$u_0(t)$')
plot(t,u[:,0]+x(t),'b',label='$y(t)$')

Answer 2

If you are allowed to use the Control System Toolbox, you can do the following:

g= tf([1 4 6 8],[1 3 5 0]);
t= 0:0.1:20;
lsim(g,sin(t),t);
figure
impulse(g)

This code produces the pictures: