ComplexInfinity/Indeterminate error when evaluating derivative

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.

• (1) It's not clear what "polynomial" this question refers to, because $1^a+\cdots+n^a$ is neither a polynomial in $n$ nor in $a$. It sounds rather like you are trying to develop a MacLaurin series in terms of $a$. (2) Have you looked at Series? Try, for example, Table[Series[Sum[i^(a), {i, 1, n}], {a, 0, 4}] // FullSimplify, {n, 1, 3}] // TableForm – whuber Sep 3 '12 at 16:35
• @whuber It seems he meant the result of $1^a +...+n^a$. It is a polynomial in n for any natural number a. – Artes Sep 3 '12 at 17:44
• Yes, @Artes, that makes sense: but taking the derivative with respect to $a$ requires $a$ to be a real number, not a natural number. I think this is goes to the heart of the difficulties, as analyzed in the answer by halirutan. – whuber Sep 3 '12 at 20:19
• @whuber In e.g. modern theoretical physics many problems require e.g. regularization techniques as a usual practice. I find this is a really nice question showing a common problem. I am a bit confused why it deserved (according to the community) only a few upvotes (far too few). – Artes Sep 3 '12 at 21:08

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= 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= 2
*)


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

D[HarmonicNumber[r, -a], {r, 3}]
(*
Out= -(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 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= 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= -(1 - a) (2 - a) a (-HarmonicNumber[n, 3 - a] + Zeta[3 - a])
*)

• Whew, thank you for your analysis on that. Yeah, I was wondering why there was an error; that explains it! – james Sep 3 '12 at 4:37

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

• Ah, thank you! I think you mean $g[j_, a_, n_\,] := D[f[n, a], \{n, j\}]$, but it works perfectly! – james Sep 3 '12 at 2:54
• @james Oops. Yep. Right in code, wrong in Answer. Fixing now. – kale Sep 3 '12 at 2:59

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].

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), {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 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 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}

• Thank you for your reply, I did just that! – james Sep 3 '12 at 4:37
• @james I extended my answer to explain the issue more extensively. – Artes Sep 3 '12 at 4:44