3
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ $$\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. $\endgroup$
    – Cesareo
    Feb 20, 2021 at 10:58

4 Answers 4

4
$\begingroup$

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.

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

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.

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

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}]

enter image description here

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Hugh
    Feb 20, 2021 at 13:31
  • $\begingroup$ @Hugh See here: reference.wolfram.com/language/ref/InverseZTransform.html. Properties & Relations last line of the code. $\endgroup$ Feb 20, 2021 at 14:33
1
$\begingroup$

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

Inverse z transform

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.

$\endgroup$
2
  • 1
    $\begingroup$ If it exists, the inverse z transform should be a discrete probability distribution and should sum to 1 $\endgroup$
    – PiE
    Feb 20, 2021 at 12:28
  • $\begingroup$ 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. $\endgroup$
    – Hugh
    Feb 20, 2021 at 12:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.