Numerical Inverse Z Transform?

Is there a way to compute the inverse z-transform in Mathematica numerically? I'm trying to compute the following...

f = InverseZTransform[(z*μ*(λ^3*μ - z*λ^3*μ + 3*λ^2*μ^2 - 4*z*λ^2*μ^2 + z^2*λ^2*μ^2 + 3*λ*μ^3 - 3*z*λ*μ^3 + μ^4))/(λ^4*μ - z*λ^4*μ + 4*λ^3*μ^2 - 6*z*λ^3*μ^2 + 2*z^2*λ^3*μ^2 + 6*λ^2*μ^3 - 9*z*λ^2*μ^3 +
3*z^2*λ^2*μ^3 + 4*λ*μ^4 - 4*z*λ*μ^4 + μ^5) /. {λ -> 0.75, μ -> 1.}, z, t];


UPDATE

So, I expanded the above polynomial in a power-series and then used the following code to compute the coefficients:

 I did try a power series expansion using the following code. The probability function does sum to 1 using the following: f = Total[

Table[ Last[
CoefficientList[
Normal[Series[ (
z \[Mu] (\[Lambda]^3 \[Mu] - z \[Lambda]^3 \[Mu] +
3 \[Lambda]^2 \[Mu]^2 - 4 z \[Lambda]^2 \[Mu]^2 +
z^2 \[Lambda]^2 \[Mu]^2 + 3 \[Lambda] \[Mu]^3 -
3 z \[Lambda] \[Mu]^3 + \[Mu]^4))/(\[Lambda]^4 \[Mu] -
z \[Lambda]^4 \[Mu] + 4 \[Lambda]^3 \[Mu]^2 -
6 z \[Lambda]^3 \[Mu]^2 + 2 z^2 \[Lambda]^3 \[Mu]^2 +
6 \[Lambda]^2 \[Mu]^3 - 9 z \[Lambda]^2 \[Mu]^3 +
3 z^2 \[Lambda]^2 \[Mu]^3 + 4 \[Lambda] \[Mu]^4 -
4 z \[Lambda] \[Mu]^4 + \[Mu]^5), {z, 0, k}]], z]], {k, 1,
100}]] /. {\[Lambda] -> 0.75, \[Mu] -> 1.0, k -> 4}

(* 1. *)


So, using the above, it appears to be a proper probability function as the terms generated by the function above sums to 1 if a sufficient number of terms are used.

The mean can be computed and turns out to be 3.05078.

This agrees with my simulation results.

• $$\frac{z \left(\frac{9 z^2}{16}-\frac{315 z}{64}+\frac{343}{64}\right)}{\frac{81 z^2}{32}-\frac{2793 z}{256}+\frac{2401}{256}}$$ is non causal. Feb 20, 2021 at 10:58

Assuming the problem is well-posed...

For a rational function, SeriesCoefficient will likely work:

w = (z*μ*(λ^3*μ - z*λ^3*μ + 3*λ^2*μ^2 - 4*z*λ^2*μ^2 +
z^2*λ^2*μ^2 + 3*λ*μ^3 - 3*z*λ*μ^3 + μ^4))/(λ^4*μ -
z*λ^4*μ + 4*λ^3*μ^2 - 6*z*λ^3*μ^2 + 2*z^2*λ^3*μ^2 +
6*λ^2*μ^3 - 9*z*λ^2*μ^3 + 3*z^2*λ^2*μ^3 + 4*λ*μ^4 -
4*z*λ*μ^4 + μ^5) /. {λ -> 3/4, μ -> 1};

SeriesCoefficient[w, {z, Infinity, t},
Assumptions -> t >= 0 && t ∈ Integers]
Table[%, {t, 0, 10}] // N
(*
{-0.986626, -2.95864, -9.09661, -28.2456, -88.0384, -274.805,
-858.257, -2681.03, -8375.66, -26166.8, -81749.8}
*)


To achieve a numerically-computed result, the exponent t must be given a value. The results are equivalent to the SeriesCoefficient results.

f = InverseZTransform[w, z, t];

Table[f /.
InverseZTransform[u_, z_, n_] :>
NIntegrate[u*z^n*I /. z -> 4 Exp[I*z], {z, 0, 2 Pi},
Method -> "Trapezoidal", MaxRecursion -> 20]/(2 Pi I),
{t, 0, 10}] /. w_Complex :> Re[w] /; Abs[Im[w]/w] < 10^-14

(*
{-0.986626, -2.95864, -9.09661, -28.2456, -88.0384, -274.805,
-858.257, -2681.03, -8375.66, -26166.8, -81749.8}
*)


The path z -> 4 Exp[I z] was chosen to encirle the poles:

Cases[f, InverseZTransform[u_, z_, n_] :>
NSolve@FunctionSingularities[u, z], Infinity]

(*  {{{z -> 1.18598}, {z -> 3.1242}}}  *)


That said, I don't know about this area and I assumed the path is in the ROC. Cesareo made a comment that seems relevant. I don't understand its significance, the comment has been ignored up til now. Perhaps it's not significant. The OP made a comment that the InverseZTransform should be a PDF, which it is clearly not above. Together, the comments make me wonder whether the problem, in particular w, has been defined correctly.

• For context...this polynomial represents the number of customers served during a busy period in a M/M/1/K queue
– PiE
Feb 21, 2021 at 16:09

Perhaps you might find useful to know how to numerically integrate ANY complex-valued function in Mathematica. I'm not familiar with Z-transforms so won't comment on that specifically, but in general, if you have any function $$f(z)$$ and you wish to integrate it over some contour $$z(t)$$, from $$t_0$$ to $$t_e$$ then simply evaluate: $$\int_{t_0}^{t_e} f(z(t))\frac{dz}{dt}$$ So in your case for $$f(z)$$, once you determine the contour say $$z(t)=re^{it}$$, numerically integrate it as:

NIntegrate[f[z] (r I Exp[I t])/.z->z[t],{t,0,2 Pi};


There are of course issues such as convergence, poles, branch-cuts, etc.

• My function is discrete...so I am not sure if this applies...it could...I just don't know
– PiE
Feb 20, 2021 at 12:25

This is not a answer your question,but how compute InverseZTransform symbolically,not numeric.

func = (z*μ*(λ^3*μ - z*λ^3*μ +
3*λ^2*μ^2 - 4*z*λ^2*μ^2 +
z^2*λ^2*μ^2 + 3*λ*μ^3 -
3*z*λ*μ^3 + μ^4))/(λ^4*μ -
z*λ^4*μ + 4*λ^3*μ^2 -
6*z*λ^3*μ^2 + 2*z^2*λ^3*μ^2 +
6*λ^2*μ^3 - 9*z*λ^2*μ^3 +
3*z^2*λ^2*μ^3 + 4*λ*μ^4 -
4*z*λ*μ^4 + μ^5);

F = SeriesCoefficient[func /. z -> 1/z, {z, 0, n}] // Simplify(*Analytical Inverse Z-Transform*)

(*Piecewise[{{μ/(2*λ + 3*μ), 1 + n == 0}, {-(((λ + μ)*(λ^2 + 4*λ*μ +
5*μ^2))/(λ*(2*λ + 3*μ)^2)), n == 0},
{(-((λ*μ*(λ + μ)^3*(2*λ + 3*μ)*((((λ + μ)^2*(λ + 4*μ - Sqrt[λ^2 +
4*μ^2]))/(λ*μ*(2*λ + 3*μ)))^n -
(((λ + μ)^2*(λ + 4*μ + Sqrt[λ^2 + 4*μ^2]))/(λ*μ*(2*λ + 3*μ)))^n))/2^n) +
2^(-1 - n)*λ*(λ^2 + 4*λ*μ + 3*μ^2)*(λ*μ*(2*λ + 3*μ)*(((λ + μ)^2*(λ + 4*μ -
Sqrt[λ^2 + 4*μ^2]))/
(λ*μ*(2*λ + 3*μ)))^(1 + n) - λ*μ*(2*λ + 3*μ)*(((λ + μ)^2*(λ + 4*μ +
Sqrt[λ^2 + 4*μ^2]))/
(λ*μ*(2*λ + 3*μ)))^(1 + n)) - 2^(-2 - n)*λ^2*μ*
(λ*μ*(2*λ + 3*μ)*(((λ + μ)^2*(λ + 4*μ - Sqrt[λ^2 + 4*μ^2]))/(λ*μ*(2*λ +
3*μ)))^(2 + n) -
λ*μ*(2*λ + 3*μ)*(((λ + μ)^2*(λ + 4*μ + Sqrt[λ^2 + 4*μ^2]))/(λ*μ*(2*λ +
3*μ)))^(2 + n)))/
(λ^2*(λ + μ)^2*(2*λ + 3*μ)*Sqrt[λ^2 + 4*μ^2]), n >= 1}}, 0] *)

Plot[F /. {λ -> 3/4, μ -> 1}, {n, 0, 2}]


• Not sure what you are doing here. You have substituted z with 1/z. In terms of z transform one is doing something similar to Laplace and inverse Laplace transforms but with discrete systems. You seem to have moved from the inside of a circle to the outside or something similar. Can you elaborate on what you are trying to do? Why are you taking the inverse of z?
– Hugh
Feb 20, 2021 at 13:31
• @Hugh See here: reference.wolfram.com/language/ref/InverseZTransform.html. Properties & Relations last line of the code. Feb 20, 2021 at 14:33

If I look at your expression in the z-plane we have

a = (1. z (5.359375 - 4.921875 z + 0.5625 z^2))/(
9.37890625 - 10.91015625 z + 2.53125 z^2)


Now I look to see its structure by taking it apart

Apart[a]


which gives

    (* -0.986626 +  0.222222 z
- 2.88291/(-3.1242 + 1. z) - 0.0757378/(-1.18598 + 1. z) *)


This has two real poles but you are going to run into trouble with the first two terms. If we take the inverse z transform term by term the first will give you a single value at the start and the second the derivative of a DiscreteDelta. I am not sure this is defined for a z transform so I am going by the equivalent in a Fourier transform.

The InverseZTransform of the poles works fine

InverseZTransform[-(
2.882905166707944/(-3.1242041285198123 + 1. z)) -
0.07573776119938679/(-1.1859810566653726 + 1. z), z,
n] // Simplify
`

Giving

We could try to do this numerically but all we would get is a huge value at the start. Can you define your problem better?

Hope that helps.

• If it exists, the inverse z transform should be a discrete probability distribution and should sum to 1
– PiE
Feb 20, 2021 at 12:28
• Probability is not my area. I use z transforms for time series analysis. However, can we get to your problem by taking the limit of something that turns your probability distribution from a continuous function into a discrete distribution? As I suggest it looks like you have a delta function and the derivative of a delta function at t = 0.
– Hugh
Feb 20, 2021 at 12:41