I've been thinking of using polynomial interpolation combined with least squares optimization to approximate various functions of a single variable, defined by linear differential or functional equations.

For differential equations the problem is rather simple and straighforward, since they are local. On the other hand, the cases such as Gamma function can be more complicated. Gamma function is defined by a non-local equation:

$$\Gamma(x+1)=x \Gamma(x)$$

Additionally, we know its values at positive integer points:


There's also the condition of log-convexity, which I will return to later.

I have applied the following approximation scheme. It can be used on any interval $[n,n+2]$, which is good enough for practical applications.

Let a polynomial $P_k(x)$ be an approximation for the function on the interval:

$$P_k(x)=a_k x^k+a_{k-1}x^{k-1}+\dots+a_1 x+a_0$$

First, it has to satisfy the constraints at integer points:

$$ \begin{cases} P_k(n)=(n-1)! \\ P_k(n+1)=n! \\ P_k(n+2)=(n+1)! \end{cases}$$

Second, it has to approximately satisfy the defining equation, in some sense. I chose the least squares method.

In other words, it has to minimize the following "Lagrangian":

$$L(a_j,\lambda_j)=\int_n^{n+1} \left(P_k(x+1)- x P_k(x) \right)^2 dx-\lambda_1 \left(P_k(n)-(n-1)! \right)- \\ - \lambda_2 \left(P_k(n+1)-n! \right)- \lambda_3 \left(P_k(n+2)-(n+1)! \right)$$

It's just the "norm" of the square of the defining equation, with constraints added with the use of Lagrange multipliers. To minimize it, we equate all the partial derivatives to zero.

Which, thanks to linearity of the defining equation, reduces to solving a system of linear equations, which has only one solution.

Thus, we find the coefficients $a_0,a_1,\dots,a_k$ and can approximate Gamma function. Or can we?

Does this algorithm allow us to approximate Gamma function on the interval $[n,n+2]$ for any positive integer $n$ with any degree of accuracy for large enough $k$? Does this approximation converge to Gamma function for $k \to \infty$?

As far as I understand, without the condition of log-convexity, Gamma function can't be uniquely defined, which is one of the reasons I have doubts. Obviously, each approximation is not log-convex, since the polynomials oscillate a lot.

Here's the example for $n=2$. I plot the error $\Gamma(x)-P_k(x)$ on the interval $[2,4]$ and list the coefficients in descending order (from $k$ to $0$).

  • $k=2$, in this case we have just enough variables to accomodate the constraints at integer points, so no optimization is involved.

$$\{a_j\}=\{3/2, -13/2, 8\}$$

enter image description here

  • $k=3$.

$$\{a_j\}=\{0.623249, -4.10924, 9.70447, -6.95797\}$$

enter image description here

  • $k=4$.

$$\{a_j\}=\{0.243426, -2.24618, 8.32714, -13.7811, 9.32815\}$$

enter image description here

  • $k=5$.

$$\{a_j\}=\{0.0762967, -0.874414, 4.20777, -10.0217, 11.8993, -4.82474\}$$

enter image description here

The method appears to converge to the Gamma function, and quite fast, since the maximum error falls by a factor of approximately $5$ with increasing $k$.


I'm not claiming this method is in any way optimal! I'm interested in it for other, more general, purposes.

  • $\begingroup$ Sorry, I'm an old man, could you try and concentrate your question a bit? Yes, $\Gamma(x+1)=x\Gamma(x)$ doesn't define it uniquely without the assumption of log-convexity, but who cares, there are other properties of that function. Whatever, your approximation scheme looks far too complicated, I use far simpler methods in my (arbitrary precision) library. $\endgroup$ – Professor Vector Jan 2 '18 at 21:06
  • $\begingroup$ @ProfessorVector, my question is highlighted. I'm not sure what kind of clarification you are asking for, I would be happy to provide any. You may also be interested in my edit at the end of the post. Thank you for your comment. $\endgroup$ – Yuriy S Jan 2 '18 at 21:10
  • $\begingroup$ Sorry, this may sound mean, but... there are well-known, established methods to compute the gamma function. Ok, in my own library, I use another one, solving a functional equation for the digamma function with polynomial approximation, integrating.... But I don't have to ask at MSE if that converges. You seem to use another scheme, and don't have a proof that converges, though your numerical evidence seems to indicate it, did I get that right? $\endgroup$ – Professor Vector Jan 2 '18 at 21:31
  • $\begingroup$ @ProfessorVector, I am not going to compute Gamma function with this method. I am interested in the method itself. This particular question deals with Gamma function, because it is well known. Sorry, for your last question, yes, this is what I'm asking $\endgroup$ – Yuriy S Jan 2 '18 at 21:39
  • $\begingroup$ It doesn't break the official rules. But there are others, violating of which guarantees zero answers. Look, people here volunteer... if they are interested. So you might want to add some motivation before the tedious staff. $\endgroup$ – Professor Vector Jan 2 '18 at 21:44

Your Answer

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

Browse other questions tagged or ask your own question.