5
$\begingroup$

It was originally asked here.

According to About the confluent versions of Appell Hypergeometric Function and Lauricella Functions the integral

$$\int_0^1t^{a-1}(1-t)^{c-a-1}(1-xt)^{-b}e^{yt}~dt $$

can be expressed as

$$\int_0^1t^{a-1}(1-t)^{c-a-1}(1-xt)^{-b}e^{yt}~dt=\dfrac{\Gamma(a)\Gamma(c-a)}{\Gamma(c)}\lim\limits_{k\to\infty}F_1\left(a,b,k,c;x,\dfrac{y}{k}\right).$$

It look like true, but gives numeric error. Look that, calculated by Mathematica:

F[a_?NumericQ, b_?NumericQ, c_?NumericQ, x_?NumericQ, y_?NumericQ] := 
 NIntegrate[t^(a - 1) (1 - t)^(c - a - 1) (1 - x t)^(-b) Exp[ y t], {t, 0, 1}]

N[F[3/2, 1, 2, .4, .3], 20]
{2.8964403550198865`}

G[a_?NumericQ, b_?NumericQ, c_?NumericQ, x_?NumericQ, y_?NumericQ] := 
 Gamma[a] Gamma[c - a]/Gamma[c] Limit[AppellF1[a, b, k, c, x, y/k],k -> Infinity]

N[G[3/2, 1, 2, .4, .3], 20] 
{2.2854650559595466`}

where I was careful to ensure that $|x|<1,|y|<1$, and $\text{Re}(c)>\text{Re}(a)>0$, which are the condition of $F_1$.Note that the results are different. Furthermore, according to Wolfram Alpha

$$\lim\limits_{k\rightarrow \infty}F_1[3/2,1,k;2;0.4,\text{any/k}]=1.4549$$

So, is It a math issue or Mathematica issue?

$\endgroup$
6
  • $\begingroup$ Question to other users: should we include the tag bugs? $\endgroup$ Commented May 2, 2018 at 19:19
  • $\begingroup$ @AccidentalFourierTransform Did you report this bug to Mathematica? Do you know if they fixed it? $\endgroup$ Commented Aug 17, 2018 at 20:53
  • $\begingroup$ In my 13.1 you now need to Rationalize to get to your N[,20] result. N[G[3/2, 1, 2, 4/10, 3/10], 20]. And F[3/2, 1, 2, 4/10, 3/10] does not give normal symbolic form (anymore?), and thus does not allow N manipulation, it is small. WorlframAlpha gives -5+5 Sqrt[5/3] for your x = 4/10 and so does 13.1. $\endgroup$ Commented Jul 4, 2022 at 16:11
  • $\begingroup$ Oh, yes! I got why it is so "cool" now for N[F[3/2, 1, 2, 0.4, .3], 20], that it cannot even listen to N[]! See, if you will change the first NIntegrate to Integrate it will actually produce an integrand from which it is calculated, but will fail to even find the antiderivative! Wow! And then you can just change Integrate to NIntegrate back again. Another bug in there (though 2.8964403550198865 is the correct result, it is G that has a problem). NIntegrate[(E^(3 t/10) Sqrt[t])/( Sqrt[1 - t] (1 - (2 t)/5)), {t, 0, 1}, PrecisionGoal -> 20, PrecisionGoal -> 50] $\endgroup$ Commented Jul 4, 2022 at 17:15
  • $\begingroup$ @ Валерий Заподовников Replacing Infinity in G by a large number, say 100 or 200, leads to a good approximation $G \simeq F$. Hance the limit is the problem. Sorry, I saw afterwards that this defect was already noticed by others. $\endgroup$ Commented Jul 1, 2023 at 13:33

1 Answer 1

7
$\begingroup$

There is a bug in the calculation of the large $k$ limit of AppellF1. This is easy to illustrate:

Needs["NumericalCalculus`"]
wrongLimit = Limit[AppellF1[3/2, 1, k, 2, 4/10, (3/10)/k], k -> Infinity];
correctLimit = NLimit[AppellF1[3/2, 1, k, 2, 4/10, (3/10)/k], k -> Infinity];
Plot[{AppellF1[3/2, 1, k, 2, 4/10, (3/10)/k], wrongLimit, correctLimit}, {k, .5, 5}
, PlotStyle -> {Blue, Directive[Red, Dashed], Directive[Green, Dashed]}, PlotRange -> {Automatic, {0, 3}}]

enter image description here

Perhaps the bug should be reported. In the meantime, you can use NLimit to get the correct limit. If you use NLimit in your definition of G, you get results that consistent with those of F (up to numerical inaccuracy).

$\endgroup$
9
  • 1
    $\begingroup$ Thank you a lot! I didn't know 'NLimit'. $\endgroup$ Commented May 2, 2018 at 19:16
  • 1
    $\begingroup$ @user3321 I'm glad I could help :-) $\endgroup$ Commented May 2, 2018 at 19:17
  • $\begingroup$ Is it normal this code NLimit[AppellF1[3/2, -1, k,1, 16(1)^2 Sin[Pi/2]^2, -(8 (1)^2 Sin[Pi/2]^2)/k], k -> Infinity] take a long time? $\endgroup$ Commented May 9, 2018 at 20:00
  • $\begingroup$ I don't know how NLimit works internally, so I can't tell. But if you want to accelerate the process, you can use SequenceLimit[AppellF1[3/2, -1, #, 1, 16, -(8/#)] & /@ N[Range[17, 35, 2]]], which yields the approximate limit 0.065 in a few seconds on my PC. It is less reliable than NLimitthough, so use at your own risk (I plotted the function and it seems that the result is more or less good, but don't take my word for it). $\endgroup$ Commented May 9, 2018 at 20:15
  • 1
    $\begingroup$ @AccidentalFourierTransform SequenceLimit is now called NumericalMath`NSequenceLimit, still being used. Why, oh, why, LOL... $\endgroup$ Commented Jul 15, 2023 at 14:31

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.