0
$\begingroup$

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.

$\endgroup$
0

2 Answers 2

5
$\begingroup$

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]

enter image description here

Compare with the un-simplified original expression:

enter image description here

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

$\endgroup$
6
  • $\begingroup$ +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. $\endgroup$
    – Bob Hanlon
    Commented Feb 10, 2021 at 20:16
  • 1
    $\begingroup$ @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. $\endgroup$
    – Roman
    Commented Feb 10, 2021 at 20:23
  • 1
    $\begingroup$ I agree, but note that you said "luckily"; for functions that do not simplify, high precision may be needed. $\endgroup$
    – Bob Hanlon
    Commented Feb 10, 2021 at 20:33
  • $\begingroup$ Thanks Roman. However, there is a problem! The function is showing an unexpected behavior in the range (roughly) z: 1000 - 2000. $\endgroup$
    – Mike
    Commented Feb 13, 2021 at 8:44
  • $\begingroup$ @Mike please see update. For $z=1000$ and above the exponential correction term is extremely small unless $x$ is close to $\pm\frac12$. $\endgroup$
    – Roman
    Commented Feb 13, 2021 at 12:40
0
$\begingroup$

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:

$\endgroup$

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.