# Unpredictable behavior of a function

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, $$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.

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)    *)


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}$$.

• +1 The original expression works by using sufficiently high arbitrary-precision, e.g., Plot[Re[fun[0.350, z]], {z, 0, 100}, PlotRange -> All, WorkingPrecision -> 50, PlotPoints -> 75]. The simplification obviates the need for the high precision. Commented Feb 10, 2021 at 20:16
• @BobHanlon using higher working precision in the face of numerical instability was the first thing we learned NOT to do in numerical analysis lectures. There's almost always a better way. Commented Feb 10, 2021 at 20:23
• I agree, but note that you said "luckily"; for functions that do not simplify, high precision may be needed. Commented Feb 10, 2021 at 20:33
• Thanks Roman. However, there is a problem! The function is showing an unexpected behavior in the range (roughly) z: 1000 - 2000.
– Mike
Commented Feb 13, 2021 at 8:44
• @Mike please see update. For $z=1000$ and above the exponential correction term is extremely small unless $x$ is close to $\pm\frac12$. Commented Feb 13, 2021 at 12:40

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.31000, 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: