## Executive Summary:

Getting an accurate answer often depends on setting the `WorkingPrecision`

high enough. Once it is high enough, though, I would expect that its exact value shouldn’t much matter. But a dot at (x,y) in this graph indicates that a call of `NDSolveValue`

in version 11.3 of Mathematica failed with a `PrecisionGoal`

of `Infinity`

, an `AccuracyGoal`

of x, and a `WorkingPrecision`

of y (with y >= x):

So the `WorkingPrecision`

must exceed the `AccuracyGoal`

by at least a few digits; and accuracies greater than 33 are not achievable. But trying for an `AccuracyGoal`

of just 11 fails when the `WorkingPrecision`

is either 46 or 65 — which is puzzling. For an `AccuracyGoal`

of 30, the `WorkingPrecision`

has to be both large and even, while accuracies of 31 or 33 require the `WorkingPrecision`

to be large and odd.

What is going on? How can I avoid selecting a value of `WorkingPrecision`

that will fail, even though smaller values would succeed?

## The details of the example:

I am using a technique of Gauss’s to compute isothermal coordinates on a Riemannian 2-manifold whose first fundamental form (E, F, G), in (u, v) coordinates, is

```
mE = (75 + 225 u^2 - 15 u^4 + 27 u^6 + 110 v^2 + 490 u^2 v^2 -
22 u^4 v^2 - 18 u^6 v^2 - 45 v^4 + 153 u^2 v^4 +
9 u^4 v^4 + 3 u^6 v^4)/2;
mF = u v (-55 - 50 u^2 - 27 u^4 - 50 v^2 + 20 u^2 v^2 +
6 u^4 v^2 - 27 v^4 + 6 u^2 v^4 + u^4 v^4);
mG = (75 + 110 u^2 - 45 u^4 + 225 v^2 + 490 u^2 v^2 +
153 u^4 v^2 - 15 v^4 - 22 u^2 v^4 + 9 u^4 v^4 + 27 v^6 -
18 u^2 v^6 + 3 u^4 v^6)/2;
```

Gauss suggests that I define a complex-valued function `q(r)`

of a real variable `r`

that is determined by a first-order ODE. (For a fuller description of the technique of Gauss’s, see the Wikipedia page on the Beltrami equation.) The following function `test`

solves the ODE to determine the isothermal coordinates of the point (s, t), setting the`PrecisionGoal`

, `AccuracyGoal`

, and `WorkingPrecision`

to `Infinity`

, x, and y:

```
test(s_, t_, x_, y_) :=
NDSolveValue(
{q'(r) == (((-mF - I Sqrt(mE mG - mF^2))/mE) /.
{u -> q(r), v -> r}),
q(t) == s}, q, {r, 0, t},
PrecisionGoal -> Infinity, AccuracyGoal -> x,
WorkingPrecision -> y)(0)
```

I then call `test`

for lots of pairs (x,y), recording which pairs generate warning messages:

```
bads(s_, t_) := Quiet(ListPlot(
Reap(Do(Check(test(s, t, x, y), Sow({x, y})),
{x, 5, 40}, {y, x, 80}))((2, 1)),
PlotRange -> All))
```

The graph above resulted from the call `bads(1/4, 1/16)`

.

I am not an expert in numerical analysis, and Mathematica’s machinery for solving differential equations is quite complicated. But the apparent existence of “unlucky” values for the `WorkingPrecision`

is both mysterious and annoying. Can I do something better than, for example, trying both y and y+1 as values for the `WorkingPrecision`

, to see if either succeeds?