Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
(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

3 Answers 3

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
share|improve this answer
    
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

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])
*)
share|improve this answer
    
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

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}
share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.