discontinuous differential equations using inverse z transform?

I have a continuse system of sys, where I would like to get a discontinious system in differential equation format.

s=tf('s')
sys = s/(s-1)^2
sysz = c2d(sys,1)

so from this I get:

sysz =
 
     2.718 z - 2.718
  ---------------------
  z^2 - 5.437 z + 7.389

but I would like to get an answer like this:

  2.718 x[n-1] - 2.718 x[n-2] = y[n] - 5.437 y[n-1] + 7.389 y[n-2]

So, I tried to use iztrans by converting the transfer function to a fractional format but did not work!

 syms z
  sysznew = (2.718 *z - 2.718)/(z^2 - 5.437 *z + 7.389)
 iztrans(sysznew)

I get a very complicated response of:

ans =
(5436*(-1)^n*7389^n*cos(n*(pi - acos((5437*8210^(1/2))/492600))))/(5437*(30*8210^(1/2))^n) - (302*kroneckerDelta(n, 0))/821 - (25388838*(-1)^n*1000^(1 - n)*4969^(1/2)*(- 4969^(1/2)/2 - 5437/2)^(n - 1))/27016453 + (25388838*(-1)^n*1000^(1 - n)*4969^(1/2)*(4969^(1/2)/2 - 5437/2)^(n - 1))/27016453

Anyone who knows how this works?

Thanks,

Answer 1

[Paul]

It appears the goal is to derive a difference equation from the corresponding, discrete-time, transfer function.

Define the transfer function

syms z n y(n) x(n) Y(z) X(z)
H(z) = vpa((2.718 *z - 2.718)/(z^2 - 5.437 *z + 7.389)),pretty(H(z))

   2.718 z - 2.718
--------------------
 2
z  - 5.437 z + 7.389

If we just take the iztrans(H(z)), we get exactly that, i.e., the impulse response.

Deriving the difference equation representation of that transfer function takes a little bit of work.

Get the numerator and denominator of H(z) and define the algebraic relationship between X(z) and Y(z) defined by the transfer function

[num,den] = numden(H(z));
eqn = Y(z)*den == X(z)*num,disp(char(eqn))

Y(z)*(1000.0*z^2 - 5437.0*z + 7389.0) == X(z)*(2718.0*z - 2718.0)

Take the iztrans of both sides

eqn = iztrans(eqn,z,n),disp(char(eqn))

7389.0*iztrans(Y(z), z, n) + 1000.0*iztrans(z^2*Y(z), z, n) - 5437.0*iztrans(z*Y(z), z, n) == 2718.0*iztrans(z*X(z), z, n) - 2718.0*iztrans(X(z), z, n)

Define Y(z) and X(z) as being the z-transforms of y[n] and x[n] respectively and sub back into the equation

Y(z) = ztrans(y(n),n,z); X(z) = ztrans(x(n),n,z);
eqn = subs(eqn),disp(char(eqn))

1000.0*y(n + 2) - 5437.0*y(n + 1) + 7389.0*y(n) + 1000.0*y(0)*iztrans(z^2, z, n) - 5437.0*y(0)*iztrans(z, z, n) + 1000.0*y(1)*iztrans(z, z, n) == 2718.0*x(n + 1) - 2718.0*x(n) + 2718.0*x(0)*iztrans(z, z, n)

Set all initial conditions to 0

eqn = subs(eqn,[y(0),y(1),x(0)],[0, 0, 0]),disp(char(eqn))

1000.0*y(n + 2) - 5437.0*y(n + 1) + 7389.0*y(n) == 2718.0*x(n + 1) - 2718.0*x(n)

Reindex into the desired form

eqn = subs(eqn,n,n-2),disp(char(eqn))

7389.0*y(n - 2) - 5437.0*y(n - 1) + 1000.0*y(n) == 2718.0*x(n - 1) - 2718.0*x(n - 2)

Set the coeffiicient of y[n] to unity. There should be a more robust way to do this, but here we know we can divide by exactly 1000.

eqn = vpa(eqn/1000),disp(char(eqn))

7.389*y(n - 2) - 5.437*y(n - 1) + 1.0*y(n) == 2.718*x(n - 1) - 2.718*x(n - 2)