Evaluate a precision-sensitive function

I want to define a nice function Φ in 3D that is computed as follows

$$\phi (x, y, z) = - z + \frac{\sqrt{32}}{\pi} \sqrt{\cosh \nu - \cos u} \sum_{n = 1}^\infty A_n \, p_n (\nu) \sin n u$$ for \begin{aligned} \nu &= \frac{1}{2} \log \left( \frac{\left( \sqrt{x^2 + y^2} + a \right)^2 + z^2}{\left( \sqrt{x^2 + y^2} - a \right)^2 + z^2} \right) \\ u &= \arctan \left( x^2 + y^2 + z^2 - a^2, 2 a z \right) \end{aligned} functions $p$ and $q$ are real combinations of half-integer Legendre polynomials $$p_n (\nu) \equiv P_{n-1/2} (\cosh \nu), \quad q_n (\nu) = Q_{n-1/2} (\cosh \nu) + \frac{i \pi}{2} P_{n-1/2} (\cosh \nu)$$

And coefficients $A$ are computed from the following infinite system of linear equations $$\Big( \sinh \nu_0 \; p_1 (\nu_0) + 2 \cosh \nu_0 \; p_1^\prime (\nu_0) \Big) A_1 + p_2^\prime (\nu_0) A_2 = \sinh \nu_0 \; q_1 (\nu_0) + 2 \cosh \nu_0 \; q_1^\prime (\nu_0) + q_2^\prime (\nu_0)$$ (for $n = 1$) \begin{multline} \sinh \nu_0 \left( - n \: q_n (\nu_0) + A_n \, p_n (\nu_0) \right) + 2 \cosh \nu_0 \left[ - n \: q_n^\prime (\nu_0) + A_n \, p_n^\prime (\nu_0) \right] + \\ + \left( - (n+1) q_{n+1}^\prime (\nu_0) + A_{n+1} \, p_{n+1}^\prime (\nu_0) \right) - \left( - (n-1) q_{n-1}^\prime (\nu_0) + A_{n-1} \, p_{n-1}^\prime (\nu_0) \right) = 0 \end{multline} (for $n > 1$)

My main goal is also to compute $V$ $$V = - \vec{\nabla} \phi$$ in Cartesian basis.

See below for my not very effective and not very precise code.

There are three main aspects:

1. Precise calculations of terms in the sum. Mathematica returns wrong results for Legendre functions for larger arguments, so how to tell MMA "hey, watch (the precision of) what you return!"? Typical example is the function $q_n (\nu)$ (see below for the definition). Without Re and highly precise argument, there will be some imaginary part for $n = 8$ and $\nu = 3$, try to compare $q'$ and $dq$. For smaller arguments and smaller values of $n$, the results will be the same, but for higher values, there will be some significant imaginary part. How to deal with this and always be sure that the result is precise to several significant digits?

2. A possible bottleneck might be the conversion between $\nu, u$ and $x, y, z$ - only way I thought of to fix this is to somehow manually simplify (for example $\cosh (\nu)$ as a rational function of $x, y, z$) and plug it in, but this would explode the code. How to effectively write this kind of functions (where the known function is a function of some variable which has to be expressed in terms of $x$, $y$, $z$, so there might be still a lot of work for MMA to do internally like plugging in logarithms, atan2 etc)? I'm afraid, this might also be the a source of numerical errors, as $(1+x^2)/2x$ is not the same as $\cosh \log x$, numerically speaking.

3. How to effectively chop off the rest of the sum, when the magnitude is small enough (like 10^-10)?? Note that each calculations of coefficients $A$ requires another linear system to be solved (even if you change their number). So it would be beneficial to compute $A$ only once for given $\kappa$, but chop the sum when the terms are too small. Alternatively, the best situation would be to find a closed-form solution for $A_n$, but I couldn't figure it out.

I am sure that the sum is behaving nicely - coefficients go to zero when computed exactly, however, MMA's problems to compute certain things precisely enough produce errors.

So I defined $p$ and $q$ (along with derivatives)

Subscript[p, n_?IntegerQ][ν_?NumericQ] :=
With[{x = SetPrecision[ν, 100]},
Re[LegendreP[n - 1/2, 0, Cosh[x]]]];
Subscript[dp, n_?IntegerQ][ν_?NumericQ] :=
With[{x = SetPrecision[ν, 100]}, (n + 1/2)/
Sinh[x] Subscript[p, n + 1][x] - (n + 1/2) Coth[x] Subscript[p,
n][x]];
Subscript[q, n_?IntegerQ][ν_?NumericQ] :=
With[{x = SetPrecision[ν, 100]},
Re[LegendreQ[n - 1/2, 0, Cosh[x]] + (I π)/
2 LegendreP[n - 1/2, 0, Cosh[x]]]];
Subscript[dq, n_?IntegerQ][ν_?NumericQ] :=
With[{x = SetPrecision[ν, 100]}, (n + 1/2)/
Sinh[x] Subscript[q, n + 1][x] - (n + 1/2) Coth[x] Subscript[q,
n][x]];


Proceeded with $A$

a[κ_?NumericQ, n_?IntegerQ] :=
a[κ, n] = (t = SetPrecision[κ, 5 n];
With[{ν0 = Log[1/t + Sqrt[1/t^2 - 1]]},
vars = Subscript[A, #] & /@ Range[1, n];
eqn0 =
N[(Sinh[ν0] (-Subscript[q, 1][ν0] +
Subscript[A, 1] Subscript[p, 1][ν0]) +
2 Cosh[ν0] (-Subscript[dq, 1][ν0] +
Subscript[A, 1]
Subscript[dp, 1][ν0]) + (-2 Subscript[dq,
2][ν0] +
Subscript[A, 2] Subscript[dp, 2][ν0]) == 0), 5 n];
eqns =
N[(Sinh[ν0] (-k Subscript[q, k][ν0] +
Subscript[A, k] Subscript[p, k][ν0]) +
2 Cosh[ν0] (-k Subscript[dq, k][ν0] +
Subscript[A, k]
Subscript[dp, k][ν0]) + (-(k + 1) Subscript[dq,
k + 1][ν0] +
Subscript[A, k + 1]
Subscript[dp, k + 1][ν0]) - (-(k - 1) Subscript[
dq, k - 1][ν0] +
Subscript[A, k - 1] Subscript[dp, k - 1][ν0]) ==
0), 5 n] /. k -> # & /@ Range[2, n];
PrependTo[eqns, eqn0];
eqns = Chop[eqns /. Subscript[A, n + 1] -> 0];
sol =
vars /. First@Solve[eqns, vars, Reals, WorkingPrecision -> 10];
N[sol, 10]]
);


Executing a[κ, n] returns $n$ coefficients (so the infinite system is cut at $n$ setting $A_{n+1} = 0$). Those converge to zero quickly.

Now I tried to compute $\phi$ in the following way

ϕ[x_?NumericQ, y_?NumericQ,
z_?NumericQ, κ_?NumericQ] := ϕ[x, y, z, κ] = (
ξ = SetPrecision[x, 20];
η = SetPrecision[y, 20];
γ = SetPrecision[z, 20];
N[-z +
Sqrt[32]/π Sqrt[Cosh[ν] - Cos[u]]
Sum[a[κ, 10][[n]] Subscript[p, n][ν] Sin[n u], {n,
1, 10}] /. {ν ->
1/2 Log[((Sqrt[ξ^2 + η^2] +
1)^2 + γ^2)/((Sqrt[ξ^2 + η^2] -
1)^2 + γ^2)],
u -> ArcTan[ξ^2 + η^2 + γ^2 - 1, 2 γ]},
10]
);


There are four arguments: $x$, $y$, $z$ and $\kappa$. There is the memoization so running $\phi (numbers)$ twice for the same numbers will return the result that was once computed.

At first, I had a problem with the precision, so I converted the input to higher precision and calculated the sum (cut at some finite $n$) in prescribed precision. Running it without these would yield inprecise result (because $p_n (\nu)$ grows very quickly with $\nu$ and the growth is faster with bigger $n$, while coefficients $a$ go to zero quicker than $p_n$ goes to infinity).

There is a huge problem here:

First /@ (a[0.1, #] & /@ Range[5, 25, 5])


Produces

{-0.006306233872, -0.006306233872, Subscript[A, 1], Subscript[A, 1],
Subscript[A, 1]}


together with the error that solve might not return a value for all variables (which it did not). Notice that the first two values are practically identical, so $n = 5$ is obviously sufficient for $\kappa = 0.1$. However, raising $\kappa$ shows that the number of coefficients needed for precise result is not constant

First /@ (a[0.5, #] & /@ Range[5, 25, 5])
{-0.2597087477, -0.2597073309, -0.2597073309, -0.2597073309, -0.2597073309}

First /@ (a[0.9, #] & /@ Range[5, 25, 5])
{-5.010880757, -4.983022342, -4.983071722, -4.983071657, -4.983071657}

First /@ (a[0.95, #] & /@ Range[5, 25, 5])
{-11.80882736, -11.64878489, -11.64949546, -11.64949310, -11.64949311}


This produces automatically problem when evaluating the function $\phi$ - how many terms should we take? For lower $\kappa$ fewer terms will suffice, in fact, increasing their number will confuse Solve and it won't return a reasonable result. But for $\kappa$ close to 1 many terms are needed to compute it precisely.

Finally, my main goal is to compute gradient in the direction $x$, $y$ and $z$ which I tried to do in the following way

Vν[x_?NumericQ, y_?NumericQ, z_?NumericQ, κ_?NumericQ] :=
Vν[x, y, z, κ] = (
ξ = SetPrecision[x, 20];
η = SetPrecision[y, 20];
γ = SetPrecision[z, 20];
N[(Cosh[ν] -
Cos[u]) (-((
Sinh[ν] Sin[u])/(Cosh[ν] - Cos[u])^2) - (Sqrt[
32]/π Sinh[ν]/(2 Sqrt[Cosh[ν] - Cos[u]])
Sum[a[κ, 10][[n]] Subscript[p, n][ν] Sin[
n u], {n, 1, Length@a[κ, 10]}] +
Sqrt[32]/π Sqrt[Cosh[ν] - Cos[u]]
Sum[a[κ, 10][[n]] Subscript[p, n]'[ν] Sin[
n u], {n, 1, Length@a[κ, 10]}])) /. {ν ->
1/2 Log[((Sqrt[ξ^2 + η^2] +
1)^2 + γ^2)/((Sqrt[ξ^2 + η^2] -
1)^2 + γ^2)],
u -> ArcTan[ξ^2 + η^2 + γ^2 - 1, 2 γ]},
10]
);
Vu[x_?NumericQ, y_?NumericQ, z_?NumericQ, κ_?NumericQ] :=
Vu[x, y, z, κ] = (
ξ = SetPrecision[x, 20];
η = SetPrecision[y, 20];
γ = SetPrecision[z, 20];
N[(Cosh[ν] -
Cos[u]) (-((
1 - Cosh[ν] Cos[u])/(Cosh[ν] - Cos[u])^2) - (Sqrt[
32]/π Sin[u]/(2 Sqrt[Cosh[ν] - Cos[u]])
Sum[a[κ, 10][[n]] Subscript[p, n][ν] Sin[
n u], {n, 1, Length@a[κ, 10]}] +
Sqrt[32]/π Sqrt[Cosh[ν] - Cos[u]]
Sum[n a[κ, 10][[n]] Subscript[p, n][ν] Cos[
n u], {n, 1, Length@a[κ, 10]}])) /. {ν ->
1/2 Log[((Sqrt[ξ^2 + η^2] +
1)^2 + γ^2)/((Sqrt[ξ^2 + η^2] -
1)^2 + γ^2)],
u -> ArcTan[ξ^2 + η^2 + γ^2 - 1, 2 γ]},
10]
);
Vx[x_?NumericQ, y_?NumericQ,
z_?NumericQ, κ_?
NumericQ] := ((1 - Cosh[ν] Cos[u])/(Cosh[ν] - Cos[u])
Cos[φ] Vν[x, y, z, κ] - (
Sinh[ν] Sin[u])/(Cosh[ν] - Cos[u])
Cos[φ] Vu[x, y, z, κ]) /. {ν ->
1/2 Log[((Sqrt[x^2 + y^2] + 1)^2 +
z^2)/((Sqrt[x^2 + y^2] - 1)^2 + z^2)],
u -> ArcTan[x^2 + y^2 + z^2 - 1, 2 z], φ ->
ArcTan[x, y]};
Vy[x_?NumericQ, y_?NumericQ,
z_?NumericQ, κ_?
NumericQ] := ((1 - Cosh[ν] Cos[u])/(Cosh[ν] - Cos[u])
Sin[φ] Vν[x, y, z, κ] - (
Sinh[ν] Sin[u])/(Cosh[ν] - Cos[u])
Sin[φ] Vu[x, y, z, κ]) /. {ν ->
1/2 Log[((Sqrt[x^2 + y^2] + 1)^2 +
z^2)/((Sqrt[x^2 + y^2] - 1)^2 + z^2)],
u -> ArcTan[x^2 + y^2 + z^2 - 1, 2 z], φ ->
ArcTan[x, y]};
Vz[x_?NumericQ, y_?NumericQ,
z_?NumericQ, κ_?
NumericQ] := (-((Sinh[ν] Sin[u])/(Cosh[ν] - Cos[u]))
Vν[x, y, z, κ] - (1 - Cosh[ν] Cos[u])/(
Cosh[ν] - Cos[u]) Vu[x, y, z, κ]) /. {ν ->
1/2 Log[((Sqrt[x^2 + y^2] + 1)^2 +
z^2)/((Sqrt[x^2 + y^2] - 1)^2 + z^2)],
u -> ArcTan[x^2 + y^2 + z^2 - 1, 2 z]};


But there is the very same problem: precision. Is there any way to do this more effectively, write less code and have it precise enough, although it should return results fast (for 3D plotting)?

• TIP: don't use subscripts. – AccidentalFourierTransform May 3 '18 at 13:16
• Solve fails for the high precision cases because you specify a low working precision. (Just get rid on the working precision option to Solve and let it determine that based on input ). Note it might be fruitful to put your equations into suitable form to use LinearSolve – george2079 May 3 '18 at 13:36
• Okay, I got rid of 'WorkingPrecision' in Solve, that works for kappa = 0.1, n = 50, however for n = 60 it fails (First@a[0.1,60] returns A1). This wouldn't be that much of a problem - who needs n = 60 for kappa = 0.1 anyway (the coefficients are very very small so only about n = 5 is needed). But for a[0.99, 60] (0.99 is about the highest value I will need) the last coefficient is -0.00001114079598 which is still too large (because the last term in the sum for some x,y,z is p_60 (nu) where nu corresponds to kappa = 0.99 which gives -0.007749299883 - way too large). – user16320 May 3 '18 at 14:21
• What I wanted to say, this needs some universal way to pick value of n, depending on the input kappa and the point in space x,y,z, so that the algorithm somehow knows that A_n p_n (nu(x,y,z)) is small enough for n = ... so it executes a[kappa, n] with this particular n and not a way bigger n or too small n. Am I clear or I just confused anyone reading this :/ – user16320 May 3 '18 at 14:23
• Why use n_?IntegerQ when good old n_Integer will do the job and do it more quickly. – m_goldberg May 4 '18 at 4:20