How to test a recursive formula using Mathematica?

Question

Through taking the derivative of $\binom{n}{k}$ w.r.t $n$ repeatedly, I found the recursive formula:

$$\frac{\partial^a}{\partial n^a}\binom{n}{k}=\sum_{j=0}^{a-1}\binom{a-1}{i}\frac{\partial^j}{\partial n^j}\binom{n}{k}\left[\psi^{(a-j-1)}(n+1)-\psi^{(a-j-1)}(n-k+1)\right]$$

and I want to test it as I am not sure if its correct. First I defined $\frac{\partial^a}{\partial n^a}\binom{n}{k}$ as $f(a,n,k):$

 f[a_,n_,k_]=D[Binomial[n,k],{n,a}];
f[a_,n_,k_]=Sum[Binomial[a-1,j] fj,n,k,{j,0,a-1}]

But this command gives zero for any value for $a,n$ and $k$.

What I did wrong? Sorry for asking a lot and thank you.

Answer 1

Without trying to solve the specific problem, here is some guidance on using recursion.

First, make sure to use the := symbol, which is SetDelayed. This will make your expressions act like a function that you can substitute values into.

Second, recursion requires that you identify the first (or possibly multiple) solution(s) and define them using the = symbol, which is Set. Using Set ensures that the value won't change unless reset. So for instance, you'll need to set

g[0, n_, k_] = g0[n,k];
g[1, n_, k_] = g1[n,k];

Now you are ready to define g[a_, n_, k_]. Since your recursive problem involves addition, you'll start with the previous solution, g[a-1, n, k], and add the incremental argument, gIncrement. gIncrement is the argument inside the summation. Here each solution is defined using := and saved using =; so that higher values of "a" can recall previous answers.

g[a_, n_, k_] := g[a, n, k] = g[a-1, n, k] + gIncrement[a-1, n, k];

Here's what the output will be:

g[#, n, k] & /@ {0, 1, 2, 3, 4} // TableForm
g0[n, k]
g1[n, k]
g1[n, k] + gIncrement[1, n, k]
g1[n, k] + gIncrement[1, n, k] + gIncrement[2, n, k]
g1[n, k] + gIncrement[1, n, k] + gIncrement[2, n, k] + 
  gIncrement[3, n, k]

Presto, you've got an efficient recursive function.

Answer 2

A first check can be done as follows:

f[a_, n_, k_] = D[Binomial[n, k], {n, a}];
g[a_, n_, k_] = 
  Sum[Binomial[a - 1, j] f[j, n, 
     k] (PolyGamma[a - j - 1, n + 1] - 
      PolyGamma[a - j - 1, n - k + 1]), {j, 0, a - 1}];
FullSimplify@Table[f[a, n, k] - g[a, n, k], {a, 1, 5}, {n, 1, 5}]

And a second one - perhaps more convincing

FullSimplify@Table[f[a, n, k] - g[a, n, k], {a, 1, 13}]

Answer 3

For a fully recursive implementation, for which @kenric gives an outline, the definition must be inductive. In particular, you need a base case (or more in some cases, like for the Fibonacci recursion).

ff[0, n_, k_] := Binomial[n, k]; (* base case *)
ff[a_?Positive, n_, k_] := 
 Sum[Binomial[a - 1, j] ff[j, n, 
    k] (PolyGamma[a - j - 1, n + 1] - 
     PolyGamma[a - j - 1, n - k + 1]), {j, 0, a - 1}];

We can speed it up, as @kenric points out, with memoization:

ff[0, n_, k_] := Binomial[n, k];
mem : ff[a_?Positive, n_, k_] := mem =
 Sum[Binomial[a - 1, j] ff[j, n, 
    k] (PolyGamma[a - j - 1, n + 1] - 
     PolyGamma[a - j - 1, n - k + 1]), {j, 0, a - 1}];

Memoization stores a computed value for reuse. The storage has a cost, but given that for large a, ff is defined in terms summing many values of ff, the savings will be substantial.

The main reasons for using := instead of = are to protect the variables in the definition from being evaluated, and in this case, to keep Sum from evaluating during the definition (which wastes a second or so only). (Sometimes you want something like Sum to be evaluated with symbolic parameters and replaced by its result during the definition of a function, but not here.) For instance, if n and k happen to have values and you use =, the function won't depend on n and k:

n = 4; k = 3;
fff[0, n_, k_] = Binomial[n, k];
?fff

Definitions fff[0, n_, k_] = 4

The value of f[0, n, k] will always be 4, no matter what n and k are.