Unpredictable behavior of a function

Question

The Hermitian matrix mat is used to construct a function fun[x, y, z] as follows:

  mat = {{1/2 + (Sqrt[1 - 4 x^2] Sinh[1/2 Sqrt[1 - 4 x^2] z])/(
     8 x^2 - 2 Cosh[1/2 Sqrt[1 - 4 x^2] z]), (
    2 I x Sinh[1/4 Sqrt[1 - 4 x^2] z]^2)/(
    4 x^2 - Cosh[1/2 Sqrt[1 - 4 x^2] z])}, {(
    2 I x Sinh[1/4 Sqrt[1 - 4 x^2] z]^2)/(-4 x^2 + 
     Cosh[1/2 Sqrt[1 - 4 x^2] z]), 
    1/2 (1 + (
       Sqrt[1 - 4 x^2] Sinh[1/2 Sqrt[1 - 4 x^2] z])/(-4 x^2 + 
        Cosh[1/2 Sqrt[1 - 4 x^2] z]))}};


t1 = D[mat, x]; t2 = mat.t1;

fun[x_, z_] = Tr[t1.t1] + (1/Det[mat]) Tr[t2.t2];

Here, $0<x<0.5$, $z\ge0$. The problem is that fun[x,z] starts getting weird after some value of z, however it is expected to show smooth behavior.

Answer 1

Your expression is numerically unstable. Luckily, simplifying it gives better results:

fun[x_, z_] = Tr[t1 . t1] + (1/Det[mat]) Tr[t2 . t2] // FullSimplify;

Plot[Re[fun[0.3, z]], {z, 0, 100}, PlotRange -> All]

Compare with the un-simplified original expression:

Update

For even larger values of $z$, an asymptotic expansion can help. It's a bit tricky in this case because of the exponential behavior, so a regular series expansion won't do.

The dominant term for large $z$ is

f1[x_, z_] = Assuming[-1/2 < x < 1/2,
  Asymptotic[fun[x, z], z -> ∞]]
(*    4/(1-4*x^2)    *)

Adding the next correction term,

f2[x_, z_] = f1[x, z] + Assuming[-1/2 < x < 1/2,
  Asymptotic[fun[x, z] - f1[x, z], z -> ∞]]
(*    4/(1-4*x^2) - 32*Exp[-Sqrt[1-4*x^2]*z/2]*x^2*z/Sqrt[1-4*x^2]    *)

For large $z$ this function $f_2(x,z)$ gives an excellent approximation if $x$ is not too close to $\pm\frac12$.

More systematically, let's set $a=\exp(z\sqrt{\frac14-x^2})$ and series-expand around $a=\infty$:

Assuming[a > 0, 
  Series[fun[x, Log[a]/Sqrt[1/4-x^2]], {a, ∞, 2}] // FullSimplify]

$$ -\frac{4}{4 x^2-1}+\frac{16 \left(4 x^2 \log (a)-8 x^2+1\right)}{a \left(4 x^2-1\right)}+\frac{16 \left(-4 x^2 \log ^2(a)+64 x^4 \log (a)-112 x^4+24 x^2-1\right)}{a^2 \left(4 x^2-1\right)}+O\left(a^{-3}\right) $$ This approximation should give very good results for large $z$, specifically for large $\log(a)=z\sqrt{\frac14-x^2}$.

Answer 2

Here's another way to refactor fun[] (the strategy is similar to how we teach computing limits at infinity: divide by the highest-order growth term):

mat = {{1/2 + (Sqrt[1 - 4 x^2] Sinh[1/2 Sqrt[1 - 4 x^2] z]) /
         (8 x^2 - 2 Cosh[1/2 Sqrt[1 - 4 x^2] z]),
        (2 I x Sinh[1/4 Sqrt[1 - 4 x^2] z]^2) / 
         (4 x^2 - Cosh[1/2 Sqrt[1 - 4 x^2] z])},
       {(2 I x Sinh[1/4 Sqrt[1 - 4 x^2] z]^2)/
         (-4 x^2 + Cosh[1/2 Sqrt[1 - 4 x^2] z]), 
        1/2 (1 + (Sqrt[1 - 4 x^2] Sinh[1/2 Sqrt[1 - 4 x^2] z])/
         (-4 x^2 + Cosh[1/2 Sqrt[1 - 4 x^2] z]))}
       } /. {
     Cosh[1/4 Sqrt[1 - 4 x^2] z] -> c, 
     Sinh[1/4 Sqrt[1 - 4 x^2] z] -> s,
     Cosh[1/2 Sqrt[1 - 4 x^2] z] -> c^2 + s^2, 
     Sinh[1/2 Sqrt[1 - 4 x^2] z] -> 2 c s} /.
   s^2 -> c^2 - 1;

t1 = D[mat, x]; t2 = mat . t1;

ClearAll[fun];
funexpr = Tr[t1 . t1] + (1/Det[mat]) Tr[t2 . t2] // 
    Factor // 
   # /. s^2 -> c^2 - 1 /. c -> 1/c & // 
  Simplify
fun[x0_?NumericQ, z0_?NumericQ] := 
  Block[{x = x0, z = z0, c = 1/Cosh[1/4 Sqrt[1 - 4 x0^2] z0], s(*=
    Sinh[1/4 Sqrt[1-4 x0^2] z0]*)}, funexpr];

(*  funexpr: c == 1/Cosh[1/4 Sqrt[1 - 4 x^2] z]
  (16 (-1 +   c^2) ( 1 + c^2 (-1 + 4 x^2))  ) /
    ( (-1 + 4 x^2) (-2 + c^2 ( 1 + 4 x^2))^2)
*)

Plot[Re[fun[0.3, z]], {z, 0, 2000}, PlotRange -> All]

There is a hard limit at machine precision, but it turns out not to be significant for fun[]. The input to Cosh needs to be under around 700. There would be a similar limit if Cosh were rewritten in terms of Exp. This limit applies to any expression used for fun[x, z] that is in terms of Cosh[]. For z at x == 0.3, the value of z can go up to around 3500. For values above this limit, c (or c^2) underflows to 0., but that does not introduce a significant error in the value of fun[0.3, z].

Cosh[u] == $MaxMachineNumber // NSolve

(*  {{u -> -710.476}, {u -> 710.476}}  *)

Cosh[1/4 Sqrt[1 - 4 (0.3)^2] z] == $MaxMachineNumber // NSolve

(*  {{z -> -3552.38}, {z -> 3552.38}}  *)

The following shows the precision loss for x = 0.3 as a function of z. It's very good:

Plot[1000 - Precision@fun[0.3`1000, SetPrecision[z, 1000]],
 {z, 0, 2000}, PlotRange -> All]

The OP's original expression for fun loses about 35 digits of precision for each increase of 100 in the value of z, which is an extremely high rate: