Recently was asked on the site the question : how to compute
$$f(t)=\int_0^t\frac{\mu^x}{\Gamma(x+1)}\ dx$$ for which I assumed $0 \leq t \leq 1$. There is no antiderivative available from any CAS I tried.
My first idea was to approximate $\Gamma(x+1)$ by a more or less constrained polynomial in $x$, then to use partial fraction decomposition. This leads to a linear combination of elliptic integrals.
Later, I found that $\Big[\Gamma(1+x)\Big]^{-1}$ looks to be very similar to a Gibbs excess energy diagram for a binary system. Some models of this physical properties are very simple (Van Laar, Margules), but they only have two parameters which is not enough for an "accurate" representation of the function. Much more accurate are Scatchard-Hildebrand, Wilson, NRTL or Uniquac models but their complexity would not allow the required integration.
So, with a limited choice, I decided to try with a model looking like the one proposed by Redlich-Kister which contains a pure polynomial component. So my idea was to write
$$\frac{1}{\Gamma(x+1)}\sim 1+x(x-1) \sum_{k=0}^p d_k\, x^k$$
$$f(t)=\frac{\mu ^t-1}{L}+$$ $$\sum_{k=0}^p (-1)^k\,d_k L^{-(k+3)} (L\, \Gamma (k+2,-L t)+\Gamma (k+3,-L t)-(k+L+2) \Gamma (k+2))$$ where $L=\log(\mu)$ (provided that $\Re(\log (\mu ))<0\land \Re(k)>-2$).
Willing to stay non-empirical, in my answer, I used $p=3$ and the $d_k$'s were computed in order to match the function and first derivative values at $x=0$, $x=\frac 12$ and $x=1$. The result were no too bad.
For the art for art's sake, I decided to use $p=6$, the $d_k$'s were computed in order to match the function, first and second derivative values at $x=0$, $x=\frac 12$ and $x=1$. This leads to
$$d_0=-\gamma \qquad d_1=-\gamma -\frac{\gamma ^2}{2}+\frac{\pi ^2}{12}$$
$$4 \sqrt \pi\,d_2=6 (178+3 \gamma (8+\gamma )) \sqrt{\pi }+32 \pi ^2-3 \pi ^{5/2}-64 \left(P^2+4 P+36\right)$$
$$3 \sqrt \pi\,d_3=-3 (1356+\gamma (111+16 \gamma )) \sqrt{\pi }-144 \pi ^2+8 \pi ^{5/2}+96
\left(3 P^2+8 P+92\right)$$
$$3 \sqrt \pi\,d_4=2 \left(3 (1376+\gamma (61+14 \gamma )) \sqrt{\pi }+156 \pi ^2-7 \pi ^{5/2}-24
\left(13 P^2+20 P+372\right)\right)$$
$$\frac{\sqrt{\pi }}{4}\,d_5=-\left(628+9 \gamma +6 \gamma ^2\right) \sqrt{\pi }-24 \pi ^2+\pi ^{5/2}+16
\left(3 P^2+2 P+84\right)$$
$$3 \sqrt \pi\,d_6=-4 \left(-6 (106+(\gamma -1) \gamma ) \sqrt{\pi }-24 \pi ^2+\pi ^{5/2}+48
\left(P^2+28\right)\right)$$ where
$P=\psi \left(\frac{3}{2}\right)=2-\gamma -2\log (2)$
The maximum absolute error is $1.5 \times 10^{-8}$ which seems to be decent.
Making the numbers rational, this would give
$$\frac{1}{\Gamma(x+1)}\sim 1+x(1-x) P_6(x)$$
$$P_6(x)=\frac{2807}{4863}-\frac{247}{3140}x-\frac{461 }{3820}x^2+\frac{66
}{1435}x^3+\frac{11 }{3303}x^4-\frac{15 }{2726}x^5+\frac{3
}{2750}x^6$$
A few values of the definite integral (just for comparison)
$$\left(
\begin{array}{cccc}
\mu & t & \text{approximation} & \text{exact} \\
0.25 & 0.2 & 0.182857067200 & 0.182857068268 \\
0.25 & 0.4 & 0.329965795201 & 0.329965797571 \\
0.25 & 0.6 & 0.443007596935 & 0.443007599310 \\
0.25 & 0.8 & 0.526656210875 & 0.526656212492 \\
0.25 & 1.0 & 0.586607844050 & 0.586607845209 \\
& & & \\
0.50 & 0.2 & 0.195704099746 & 0.195704100926 \\
0.50 & 0.4 & 0.376453536397 & 0.376453539149 \\
0.50 & 0.6 & 0.535922898037 & 0.535922900791 \\
0.50 & 0.8 & 0.671420159082 & 0.671420160584 \\
0.50 & 1.0 & 0.782934567076 & 0.782934567750 \\
& & & \\
0.75 & 0.2 & 0.203784473259 & 0.203784474510 \\
0.75 & 0.4 & 0.407825165336 & 0.407825168343 \\
0.75 & 0.6 & 0.602995831714 & 0.602995834721 \\
0.75 & 0.8 & 0.782793915893 & 0.782793917219 \\
0.75 & 1.0 & 0.943235581611 & 0.943235581767
\end{array}
\right)$$
I am more than happy with these results but,again, just for the art for art's sake, I would like to be even better. For sure, I could use $p=9$ and obtain the parameters in order to match the function, first, second and third derivative values at $x=0$, $x=\frac 12$ and $x=1$; this would lead to monstreous expressions.
So, my question is : without any curve fit and avoiding so many terms, is there a way to obtain a better (mathematically based) approximation of $\Big[\Gamma(1+x)\Big]^{-1}$ for $0 \leq x \leq 1$
Best Answer
$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} &\bbox[5px,#ffd]{1 \over \Gamma\pars{x + 1}} = {1 \over x!} = {x + 1 \over \pars{x + 1}!} = {\pars{x + 2}\pars{x + 1} \over \pars{x + 2}!} \\[5mm] = &\ \cdots = {\prod_{k = 1}^{n}\pars{x + k} \over \pars{x + n}!} \\[5mm] \approx &\ {\prod_{k = 1}^{n}\pars{x + k} \over \root{2\pi}\pars{x + n}^{x + n + 1/2}\,\,\,\expo{-n - x}} = {\tt mySol}\pars{x,n} \end{align}