ComplexInfinity/Indeterminate error when evaluating derivative

Question

Motivation: I want to find the coefficients for the polynomial that is obtained when one adds together the first $n$ natural numbers to the power of $a$; that is, when you consider $1^{a} + 2^{a} + \cdots + n^{a}$. The way I wanted to do this was: find the polynomial, then take the $j$-th derivative, evaluate at 0, and divide by $j!$. My code to just obtain the polynomials and to take the derivatives are:

f[n_, a_] = Sum[i^(a), {i, 1, n}]
g[d_, a_, n_] = Derivative[d, 0][f][n, a]

So, for example, $g[1,2,n]$ is the first derivative (with respect to $n$) of the polynomial corresponding to $1^2 + 2^2 + \cdots + n^2$.

The Problem: If I try $g[1,2,r]$ or $g[2,2,r]$, everything is fine. If I attempt to try $g[3,2,r]$ (the third derivative, which should just be 2) I get the error:

Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

I'm not sure what I'm doing wrong; this is true for $g[d,a,n]$ in general, if $d = a+1$. I'm using Mathematica 8.

Answer 1

Lets start step by step. First let us see that your sum is transformed into an build in function

Sum[i^a, {i, 1, r}]
(*  
  Out[36]= HarmonicNumber[r, -a]
*)

Now we can go two different ways: One way is to derivate this expression wrt r and keep the unknown a. Second is, we substitute a->2 and derive afterwards. If we replace a by the number 2 we get the expected result

D[HarmonicNumber[r, -2], {r, 3}]
(*    
  Out[35]= 2
*)

If we don't replace it, we get another expression

D[HarmonicNumber[r, -a], {r, 3}]
(*
  Out[42]= -(1 - a) (2 - a) a (-HarmonicNumber[r, 3 - a] + Zeta[3 - a])
*)

If we replace a->2 in this expression, we get in the last part Zeta[1] which is infinity. I first thought this is maybe a bug in D which uses wrong assumptions. Indeet, if we don't stupidly substitute a by 2 but we calculate the Limit we get

Limit[-(1 - a) (2 - a) a (-HarmonicNumber[r, 3 - a] + Zeta[3 - a]), a -> 2]
(*
  Out[57]= 2
*)

Considering this, you could have ensured that first the numbers are inserted into your sum and then the derivation is done. As example you could have used this here

g[d_, a_] := Block[{n}, D[Sum[i^a, {i, 1, n}], {n, d}]]

which is neither perfect nor foolproof but it works for your example when used with numbers. Note that it still gives the same expression for an unknown a

g[3, a]
(*
  Out[55]= -(1 - a) (2 - a) a (-HarmonicNumber[n, 3 - a] + Zeta[3 - a])
*)

Answer 2

Let's build the polynomial:

f[n_, a_] := Sum[i^(a), {i, 1, n}]

Now, build a function to take the j-th derivative of f:

g[j_, a_, n_] := D[f[n, a], {n, j}]

So with a=2:

f[n, 2] = 1/6 n (1 + n) (1 + 2 n)

Assuming the derivative is with respect to n:

g[1, 2, n]

(*1/3 n (1 + n) + 1/6 n (1 + 2 n) + 1/6 (1 + n) (1 + 2 n)*)

g[2, 2, n]

(*(2 n)/3 + 1/3 (1 + 2 n + 2 (1 + n))*)

g[3, 2, n]

(*2*)

To evaluate at zero and divide by j!:

g[j, 2, n]/j! /. j -> 2 /. n -> 0

Answer 3

The origin of the problem is a small bug appearing in your code here : Derivative[1, 0][HarmonicNumber] (even though this expression evaluates correctly in general, however you should have rather Derivative[1, 0][#1 &][r, 0]! ) for [r,0].

Since we've had :

f[n_, a_] = Sum[i^(a), {i, 1, n}]

HarmonicNumber[n, -a]

This result is a polynomial in n for any natural number a, e.g. :

And @@ (PolynomialQ[#, n] & /@ Sum[i^(Range[100]), {i, 1, n}])

True

Then if you evaluate g[d, a, n] setting numeric values d and a such that d == a + 1, subsequently your system will encounter the message evaluating :

Derivative[1, 0][HarmonicNumber][r, 0]

Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

Indeterminate

You can find out this issue checking the output of e.g. :

g[3, 2, r] // Trace

Zeta[1] alone (which yields ComplexInfinity) doesn't produce the message, however this derivative :

Derivative[1, 0][HarmonicNumber]

#2 (Zeta[#2 + 1] - HarmonicNumber[#1, #2 + 1]) &

does when evaluated on [r,0] i.e. it yields 0 Zeta[1] which is Indeterminate.

One can get rid of the problem simply defining e.g.

f[n_, a_] = Sum[ i^(a), {i, 1, n}]; 
g1[d_, a_, n_] := D[ f[n, a], {n, d}]

i.e. substituting Set (=) in g by SetDelayed (:=) and Derivative operator by D. See e.g. Assignments ( lhs = rhs — immediate assignment, with rhs evaluated at the time of assignment, and lhs := rhs — delayed assignment, with rhs reevaluated every time it is used ). Moreover Derivative is a functional operator acting on functions unlike D which acts on expressions. This post on Wolfram Blog you may find useful : Three Functions for Computing Derivatives.

Applying the above definition of g1 we get :

{ g1[5, 5, n], g1[3, 2, n], g1[4, 3, n], g1[3, 10, n], g1[21, 20, n]}

{ 60 + 120 n, 2, 6, -3 + 60 n^2 - 210 n^4 + 420 n^6 + 360 n^7 + 90 n^8, 2432902008176640000}