Possible bug in hypergeometric function AppellF1

In the context of sums over Legendre polynomials ((1), (2)) I stumbled upon the interesting hypergeometric function AppellF1[].

Unfortunately, the implementation in Mathematica appears to have difficulties when it comes to numerical evaluation, including plots.

This behaviour spoils to a large degree the "success" of having obtained closed form expressions.

I would consider this a bug.

$Version (* Out[2944]= "10.1.0 for Microsoft Windows (64-bit) (March 24, 2015)" *)  $$\\$$ Example 1: Plot3D, real arguments As a compact example, consider this plotting command Plot3D[{Re[AppellF1[-(1/2), -(1/2), -(1/2), 1/2, x, y]], Im[AppellF1[-(1/2), -(1/2), -(1/2), 1/2, x, y]]}, {x, -1, 1}, {y, -1, 1}]  We see that the "corners" are cut out despite the fact that AppellF1[] is well defined for Abs[x] < 1 and Abs[y] < 1. $$\\$$ Example 2: numerical value at a specific points 2a. What is the numerical value of this expression (which appears naturally in the summation problem mentioned in the beginning) AppellF1[-(1/2), -(1/2), -(1/2), 1/2, E^((2 I π)/3), E^(-((2 I π)/3))] // N (* Out[2973]= AppellF1[-0.5, -0.5, -0.5, 0.5, -0.5 + 0.866025 I, -0.5 - 0.866025 I]  $$\\$$ 2b. Particular results in the "buggy" region of Plot3D AppellF1[-(1/2), -(1/2), -(1/2), 1/2, 0.9, 0.9] // N (* Out[3091]= 1.9 *)  There seems to be no problem. Generally, these special cases do not have difficulties AppellF1[-(1/2), -(1/2), -(1/2), 1/2, x, x] // FullSimplify (* Out[3093]= 1 + x *) AppellF1[-(1/2), -(1/2), -(1/2), 1/2, x, -x] // FullSimplify (* Out[3094]= Hypergeometric2F1[-(1/2), -(1/4), 3/4, x^2] *)  But consider this case with x = 0.9 and y = 0.8: AppellF1[-(1/2), -(1/2), -(1/2), 1/2, 0.9, 0.8] // N (* Out[3119]= AppellF1[-0.5, -0.5, -0.5, 0.5, 0.9, 0.8] *)  Here no numerical value is returned. Whereas, e.g. AppellF1[-(1/2), -(1/2), -(1/2), 1/2, 0.5, 0.6] // N (* 1.55065 *)  $$\\$$ Example 3a: Plot This plot takes about 5 minutes on my PC: Plot[Re[AppellF1[-(1/2), -(1/2), -(1/2), 1/2, E^(I θ), E^(-I θ)] /. {d -> 1/2}], {θ, 0, 2 π}]  Notice the numerical difficulties close to the ends of the interval and the broad gap in the middle. Example 3b: ListLinePlot The coarse-grained alternative is Table[{θ, AppellF1[-d, -(1/2), -(1/2), 1 - d, E^(I θ), E^(-I θ)] /. {d -> 1/2}}, {θ, 0, 2 π, π/10}]; Chop[N[%], 10^-5] (* Out[2993]= {{0, 2.}, {0.314159, 1.90825}, {0.628319, 1.7056}, {0.942478, 1.44035}, {1.25664, 1.14516}, {1.5708, 0.847213}, {1.88496, 0.569869}, {2.19911, 0.332818}, {2.51327, 0.151856}, {2.82743, AppellF1[-0.5, -0.5, -0.5, 0.5, -0.951057 + 0.309017 I, -0.951057 - 0.309017 I]}, {3.14159, 0}, {3.45575, AppellF1[-0.5, -0.5, -0.5, 0.5, -0.951057 - 0.309017 I, -0.951057 + 0.309017 I]}, {3.76991, 0.151856}, {4.08407, 0.332818}, {4.39823, 0.569869}, {4.71239, 0.847213}, {5.02655, 1.14516}, {5.34071, 1.44035}, {5.65487, 1.7056}, {5.96903, 1.90825}, {6.28319, 2.}} *) ListLinePlot[%]  The numerical chopping takes a few minutes before we can ListLinePlot the list. We can see that the gap in the middle of the graph is caused by non evaluated functions AppellF1. $$\\$$ Some definitions (cf. (3)) Series representation $$F1(a,b_1,b_2,c,x,y) = \sum _{m=0}^{\infty } \sum _{n=0}^{\infty }\frac{ (a)_{m+n} }{(c)_{m+n}} (b_1)_m (b_2)_n \frac {x^m}{m!} \frac {y^n}{n!}$$ Sum[ Pochhammer[a, m + n]/ Pochhammer[c, n + m] Pochhammer[b1, m] Pochhammer[b2, n] x^m/m! y^n/ n!, {m, 0, ∞}, {n, 0, ∞}] (* Out[2843]= AppellF1[a, b1, b2, c, x, y] *)  Integral representation ((3) formula (9)) $$F1(a,b_1,b_2,c,x,y) = \frac{\Gamma (c)}{\Gamma (a) \Gamma (c-a)} \int_0^1 t^{a-1} (1-t)^{-a+c-1} (1-t x)^{-b_1} (1-t y)^{-b_2} \, dt$$ fctF1[a_, b1_, b2_, c_, x_, y_] := Gamma[c]/(Gamma[a] Gamma[c - a]) Integrate[ t^(a - 1) (1 - t)^(c - a - 1) (1 - x t)^-b1 (1 - t y)^-b2, {t, 0, 1}]  References • I would most certainly report this to WRI. – ktm Mar 6, 2017 at 19:19 • Because AppellF1[-(1/2), -(1/2), -(1/2), 1/2, 0.9, -0.9] outputs$1.15196\$, this seems to be a graphic problem. Mar 6, 2017 at 19:28
• @user64494 No, try this AppellF1[-(1/2), -(1/2), -(1/2), 1/2, x, y] /. {x -> 0.9, y -> 0.91} // N (* Out[3023]= AppellF1[-0.5, -0.5, -0.5, 0.5, 0.9, 0.91] *) The function is not evaluated. That's the bug. Mar 7, 2017 at 15:19

The method I am about to propose is a bit finicky, and thus requires some manual intervention and checking. With that caveat, here is a method to evaluate the first Appell hypergeometric function for arguments where the built-in AppellF1[] fails.

This is a combination of two methods: a method, due to Cuyt, for transforming a double series into a sequence, and the Wynn $$\varepsilon$$ algorithm.

Normally, one would use Wynn through the built-in undocumented function NumericalMathNSequenceLimit[] (SequenceLimit[] in older versions), but I have found it unstable in this situation. Thus, I fell back to using my own implementation, implemented as the van den Broeck-Schwartz form:

wgvs[seq_?VectorQ, h_: 1] := Module[{n = Length[seq], ep, v, w},
Table[ep[k] = seq[[k]]; w = 0;
Do[v = w; w = ep[j];
ep[j] = v If[OddQ[k - j], h, 1] + 1/(ep[j + 1] - w),
{j, k - 1, 1, -1}];
ep[Mod[k, 2, 1]], {k, n}]]


Now, let's present the method. For reference, let's take a case where the built-in AppellF1[] evaluates numerically, for comparison purposes:

With[{a = -1/2, b1 = -1/2, b2 = -1/2, c = 1/2, x = 3/4, y = 2/3},
res = N[AppellF1[a, b1, b2, c, x, y], 30]]
1.70889493846451655825978155908


Here is Cuyt's method for summing the defining double series:

With[{a = -1/2, b1 = -1/2, b2 = -1/2, c = 1/2, x = 3/4, y = 2/3},
seq = N[Accumulate[Table[Sum[Function[{i, j},
Pochhammer[a, i + j]/Pochhammer[c, i + j]
Pochhammer[b1, i] Pochhammer[b2, j]
x^i/i! y^j/j!] @@ v,
{v, FrobeniusSolve[{1, 1}, k]}], {k, 0, 39}]], 60]];


Note the need to use 1. a large number of terms, and 2. very high precision. This is because there will inevitably some numerical cancellation during the application of Wynn $$\varepsilon$$.

You might be tempted to just take the last entry of seq, but that's not a very accurate result…

Last[seq] - res
-1.3407389500737610211*10^-10


…yet. This is where Wynn $$\varepsilon$$ comes in:

eps = wgvs[seq];

Take[eps, -10] - res
{-4.78916651602*10^-18, -1.42635137821*10^-18, -3.7802058433*10^-19,
-7.57389529*10^-21, 1.0837220906*10^-19, 1.4919244639*10^-19,
1.6200922471*10^-19, 1.6650702861*10^-19, 1.6792388370*10^-19,
1.6841995299*10^-19}


and we see that the results are much better after acceleration.

Now, let's use the OP's complex example, where AppellF1[] is not able to evaluate numerically:

With[{a = -1/2, b1 = -1/2, b2 = -1/2, c = 1/2, x = Exp[2 I π/3], y = Exp[-2 I π/3]},
N[AppellF1[a, b1, b2, c, x, y], 20] // Head] (* does not evaluate *)
AppellF1

With[{a = -1/2, b1 = -1/2, b2 = -1/2, c = 1/2, x = Exp[2 I π/3], y = Exp[-2 I π/3]},
seq = N[Accumulate[Table[Sum[Function[{i, j},
Pochhammer[a, i + j]/Pochhammer[c, i + j]
Pochhammer[b1, i] Pochhammer[b2, j]
x^i/i! y^j/j!] @@ v,
{v, FrobeniusSolve[{1, 1}, k]}], {k, 0, 39}]], 60]];


Since we don't know the result in advance, we will be a bit more conservative. First, take the differences of the transformed partial sums:

Differences[Take[seq, -10]] // Chop // N
{-0.0000182151, -0.0000474801, 0.0000594675, -0.000014441, -0.000037797,
0.0000476935, -0.0000116771, -0.0000306683, 0.0000389418}


and we see that they only agree to four or so digits. Now, apply wgvs[]:

eps = wgvs[seq];


and look at the differences again:

Differences[Take[eps, -10]] // N
{-4.86267*10^-17 + 0. I, 3.14393*10^-17 + 0. I, 6.35936*10^-18 + 0. I,
-1.57739*10^-18 + 0. I, 9.1468*10^-20 + 0. I, -1.62375*10^-19 + 0. I,
2.81443*10^-19 + 0. I, -1.45097*10^-20 + 0. I, 3.98562*10^-21 + 0. I}


and we now get agreement to twenty or so digits, so we can take the last one as a provisional result:

Last[eps] // Re
0.40629888645996024661252095193


As it turns out, this Appell value that cannot be evaluated by Mathematica has a relatively simple closed form in terms of complete elliptic integrals (derived after some amount of used paper):

N[2 EllipticE[1/4] - 3/2 EllipticK[1/4], 30]
0.406298886459960246612785047283


and we see that the double sum evaluated using the accelerated Cuyt transformation has indeed given twenty digits of accuracy.

More generally, we have the result

AppellF1[-1/2, -1/2, -1/2, 1/2, x, y] ==
y AppellF1[1/2, -1/2, 1/2, 3/2, x, y] +
x AppellF1[1/2, -1/2, 1/2, 3/2, y, x] + Sqrt[1 - x] Sqrt[1 - y]


which can be readily expressed in terms of the Carlson integrals (Mathematica package from here):

<<Carlson

4 CarlsonRG[1 - x, 1 - y, 1] + (x + y - 2) CarlsonRF[1 - x, 1 - y, 1] -
Sqrt[1 - x] Sqrt[1 - y]


or less elegantly in terms of more conventional special functions,

2 Sqrt[x] EllipticE[ArcCos[Sqrt[1 - x]], y/x] +
(y/Sqrt[x] - Sqrt[x]) InverseJacobiCN[Sqrt[1 - x], y/x] +
Sqrt[1 - x] Sqrt[1 - y]


As mentioned, I have found the current version of the accelerated Cuyt method to be finicky and a bit hard to automatize; there are argument ranges where you need to take more terms, higher precision, or both, and I usually make these judgments on convergence based on the behavior of the original and transformed sequences. Nevertheless, someone might be able to build a usable routine based on this someday.