6
$\begingroup$

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}]

plot

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 π}]

plot

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[%]

enter image description here

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

(1) https://math.stackexchange.com/questions/2111566/summation-of-legendre-polynomials-sum-l-2-infty-frac2l1l1kl1/2116383#2116383

(2) Summation of Legendre polynomials related to the zeta function

(3) http://mathworld.wolfram.com/AppellHypergeometricFunction.html

$\endgroup$
  • $\begingroup$ I would most certainly report this to WRI. $\endgroup$ – user6014 Mar 6 '17 at 19:19
  • $\begingroup$ Because AppellF1[-(1/2), -(1/2), -(1/2), 1/2, 0.9, -0.9] outputs $1.15196$, this seems to be a graphic problem. $\endgroup$ – user64494 Mar 6 '17 at 19:28
  • $\begingroup$ @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. $\endgroup$ – Dr. Wolfgang Hintze Mar 7 '17 at 15:19
8
$\begingroup$

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 NumericalMath`NSequenceLimit[] (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

Addendum 7/25/2019

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.