4407 days ago by SC2009

Computational and Experimental Mathematics

Introduction: Some Sums

Neil Calkin and Dan Warner

November 14, 2009

# Consider an arbitrary arithmetic progression alfa = 5 beta = 2 # Work through these examples and then experiment # with your own values for alfa and beta # Choose a convenient number of terms for your experiment n = 50 ap = range(alfa, n, beta) print 'The Arithmetic Progression' print ap 
# Here's an alternate way to define the arithmetic progression ap2 = [(alfa+beta*j) for j in range(n)] print ap2 


The variable $n$ played a slightly different role here. In the first case it was the upper limit, and in this case it is was the number of terms, while $x$ takes on the values 0 through $n-1$.

# Consider the first differences np = len(ap) var('y') d1 = differences(ap) print 'The sequence of first differences' print d1 


We can conclude that this defines an arithmetic progression of the form

$\alpha + \beta \,x$

with   $\alpha$ = ap[0]   and   $\beta$ = d1[0] = d1[1] = ....

This is the exact discrete counterpart to solving the first order differential equation

$\frac{d}{dt}x(t) = \beta$,    $x(0) = \alpha$

where $\beta$ is a constant.

# Now consider the sequence of sums print 'The Arithmetic Progression' print ap np = len(ap) s1 = copy(ap) for k in range(1,np): s1[k] += s1[k-1] print 'The sequence of sums' print s1 
d1 = differences(s1) print 'The sequence of first differences' print d1 d2 = differences(d1) print 'The sequence of second differences' print d2 


By analogy with the second order differential equation

$\frac{d^2}{dt^2}x(t) = \gamma$   $x(0) = \alpha$   $\frac{d}{dt}x(0) = \beta$

with   $\alpha$ = s1[0],   $\beta$ = d1[0],   and  $\gamma$ a constant.

To get the solution to the differential equation we integrate twice.  The first integration gives us

$\frac{d}{dt}x(t) = \gamma x + \beta$

and the second integration gives us

$x(t) = \frac{1}{2} \gamma x^2 + \beta x + \alpha$

So, is there a similar formula in $n$ for the discrete case?  If so, should it be of the form 

$a n^2 + b n + c$ ?

Three Experimental Math ideas:

First, paste the sequence into the Online Encyclopedia of Integer Sequences 


For $\alpha =5$ and $\beta=2$ this is sequence A028347, which has the formula $n^2 - 4$ with an offset of 2.  That is, $(n+2)^2-4$ for the sequence as displayed. 

Since our sequence starts with the second term of A028347 our formula is $(n+3)^2 - 4 = n^2 + 6n + 5$

We also see that this sequence also occurs in a number of other situations including the Balmer series of the Hydrogen atom, and there is even a Sage formula for generating the sequence based on the Lucas numbers.


print 'Adjusted OEIS formula' ap3 = [(s^2+6*s+5) for s in range(np)] print ap3 print print 'OEIS connection with the Lucas numbers' [lucas_number2(2, n, 2) for n in xrange(2, 44)] 

Second, recall that you probably saw this problem a long time ago as

$S_n = \alpha + (\alpha + \beta) + (\alpha + 2\beta) + \cdots + (\alpha + n\beta) = \sum_{k=0}^{n} \alpha + k \beta$

which has a nice geometric interpretation as a staircase that covers half the rectangle of height $2 \alpha + n \beta$ and width $n+1$.

n = 5 staircase = [(0,0)] for k in range(n+1): staircase.append((k,ap[k])) staircase.append((k+1,ap[k])) staircase.append((n+1,0)) g = Graphics() g += polygon(staircase) top = ap[n]+ap[0] g += line([(0,0),(0,top),(n+1,top),(n+1,0)]) show(g) 

In other words

$S_n = \frac{1}{2}(2\alpha + n \beta) (n+1)) = \frac{1}{2}(\alpha + [\alpha+n \beta]) (n+1)$

or "one half of the first plus the last times the number of terms".  For the case $\alpha = 5$ and $\beta = 2$ we get

$S_n = \frac{5 + (5 + n 2)}{2}(n+1) = (5+n)(n+1) = n^2 + 6n + 5$


Third, you can assume that the sum is a quadratic polynomial of the form

$S_n = a n^2 + b n + c$ 

and solve for the coefficients using three different values of ($n$, $S_n$).

var('x,a,b,c') print s1[:3] eq1 = (a*x^2+b*x+c)(x=0) == s1[0] eq2 = (a*x^2+b*x+c)(x=1) == s1[1] eq3 = (a*x^2+b*x+c)(x=2) == s1[2] sol = solve([eq1,eq2,eq3],a,b,c,solution_dict=True) coef = sol[0] print coef 
poly = coef[a]*x^2+coef[b]*x+coef[c] ap4 = [poly(x=k) for k in range(np)] print ap4 print s1 

The next step would then be to show that the resulting formula works for all values of $n$.  This provides an excellent opportunity to use mathematical induction.

Base step: $S_0 = \frac{1}{2}(2\alpha + 0 \beta) (0+1)) = \alpha $

Inductive step:  

$S_{n+1} = S_n + \alpha + (n+1) \beta = \frac{1}{2}(\alpha + [\alpha + n \beta]) (n+1)  + \alpha + (n+1) \beta = \frac{1}{2}\left( (\alpha + [\alpha + n \beta]) (n+1)  + 2\alpha +  2(n+1) \beta \right)$

$S_{n+1} = \frac{1}{2} \left( \alpha(n+2) + \alpha(n+2)+(n+2)\beta(n+1) \right) = \frac{1}{2} (\alpha + [\alpha+(n+1)\beta])(n+2) $


More Sums

Let's now explore binomial expansions, Pascal's triangle, and more interesting sequences of sums.

var('x,y') p = [1] n = 7 for k in range(1,n): p.append(expand((x+y)*p[k-1])) print p[k] 
# That leads directly to Pascal's triangle n = 13 Pt = matrix(ZZ,n) Pt[0,0] = 1 for i in range(1,n): Pt[i,0] = 1 for j in range(i): Pt[i,j+1] = Pt[i-1,j] + Pt[i-1,j+1] print Pt 
pc = Pt.columns() # print pc s0 = copy(pc[0]) for k in range(1,len(s0)): s0[k] += s0[k-1] print s0 s1 = copy(pc[1]) for k in range(1,len(s1)): s1[k] += s1[k-1] print s1 s2 = copy(pc[2]) for k in range(1,len(s2)): s2[k] += s2[k-1] print s2 s3 = copy(pc[3]) for k in range(1,len(s3)): s3[k] += s3[k-1] print s3 

Observation 1:

If we sum the first $N$ elements in column $k$ that sum will be the value of the element in the next row and column, namely, row $N+1$ and column $k+1$.

$\sum\limits_{n=k}^N P_{n,k} = P_{N+1,k+1}$

Since $P_{n,k} = \left( {\begin{array}{*{20}c} n \\ k \\ \end{array}} \right)$ this is the same as the statement

$\sum\limits_{n=k}^N \left( {\begin{array}{*{20}c} n \\ k \\ \end{array}} \right) = \left( {\begin{array}{*{20}c} N+1 \\ k+1 \\ \end{array}} \right)$

or, if we use the factorial definition of $n$ choose $k$

$\sum\limits_{n=k}^N \frac{n!}{(n-k)! k!} = \frac{(N+1)!}{(N-k)!(k+1)!}$

Actually, this observation holds for all cases, and the Theorem is easy to prove by mathematical induction. 


Observation 2:

Look at each row of Pascal's Triangle.  In all the cases in the Table, if $n$ is a prime (2, 3, 5, 7, 11), then $n$ divides all the entries in the row between the first and last, which are ones.  On the other hand, when $n$ is not a prime (4, 6, 8, 9, 10), then $n$ does not divide every entry in the row.  In particular, 4 doesn't divide 6; 6 doesn't divide 15 or 20; 8 doesn't divide 28; 9 doesn't divide 84; 10 doesn't divide 45; and 12 doesn't divide 66.

This observation also holds for all cases, and the result along with Eisenstein's criterion can be used to show that polynomials of the form

$x^{p-1} + x^{p-2} + \cdots + x + 1$

cannot be factored over the rational numbers when $p$ is a prime.

# Now let's explore some partial row sums nc = 6 Pr = matrix(ZZ,n,nc) for j in range(nc): for i in range(n): if j == 0: Pr[i,j] = Pt[i,j] else: Pr[i,j] = Pt[i,j] + Pr[i,j-1] print Pr 
pc = Pr.columns() 
# A000124 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[2] 
# ?? Numbers g = Graphics() g += circle((2,2),2) g += line([(0,0),(4,4)]) g += line([(0,2),(4,2)]) g += line([(0,4),(4,0)]) # g += line([(0,3.5),(3.5,0)], rgbcolor=(1,0,0)) show(g,aspect_ratio=1) 
# A000125 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[3] 
# A000127 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[4] 
np = 3 g = Graphics() g += circle((2,2),2) ck = cos(2.0*pi/np) sk = sin(2.0*pi/np) xk = 0 yk = 2 points = [(xk+2,yk+2)] for k in range(1,np): xt = ck*xk - sk*yk yk = sk*xk + ck*yk xk = xt points.append((xk+2,yk+2)) g += point(points,pointsize=20) for i in range(np): for j in range(i,np): g += line([points[i],points[j]]) show(g,aspect_ratio=1) 
# A006261 numbers -- Check with the Online Encyclopedia of Integer Sequences print pc[5] 
# A000045 numbers diag_sums = [] for i in range(n): s = 0 for j in range(i+1): s += Pt[i-j,j] diag_sums.append(s) print diag_sums 
F = transpose(Matrix(diag_sums)) print F 
R = Matrix(ZZ,n) for i in range(n): R[n-1-i,i] = 1 
F2 = Pt*R*F print F2 
F3 = Pt*R*F2 print F3 
print (F[n-1]/F[n-2]).n(digits=30) print (F2[n-1]/F2[n-2]).n(digits=30) print (F3[n-1]/F3[n-2]).n(digits=30) 
# Check the inverse symbolic calculator # Original: http://oldweb.cecm.sfu.ca/projects/ISC/ISCmain.html # New version: http://glooscap.cs.dal.ca:8087/ # ISC's Best Guess is fn = (F3[n-1]/F3[n-2]).n(digits=30) phi = (1 + sqrt(5))/2 phin = phi.n(digits=30) relerr = (fn-phin)/phin print phin,relerr 
# Golden Ratio is defined by GR = x/1 == 1/(x-1) show(GR) sol = solve(GR,x,solution_dict=True) show(sol) 
phi = sol[1][x] show(phi) print phi.n(digits=1000) 
F1 = 1 F2 = 1 Fnm1 = F1 Fn = F2 for k in range(75): Fnp1 = Fn + Fnm1 print (Fnp1/Fn).n(digits=30),' ',(Fn/Fnp1).n(digits=29) Fnm1 = Fn Fn = Fnp1 
var('phi') Eq2 = phi^2 == phi + 1 print Eq2 print expand(phi*Eq2) 
Eq3 = phi^3 == (phi+1) + phi print Eq3 print expand(phi*Eq3) 
Eq4 = phi^4 == expand(2*(phi+1) + phi) print Eq4 print expand(phi*Eq4) 
Eq5 = phi^5 == expand(3*(phi+1) + 2*phi) print Eq5 print expand(phi*Eq5) 
Eq6 = phi^6 == expand(5*(phi+1) + 3*phi) print Eq6 
# Here's a cute matrix identity PP = Pt * transpose(Pt) print PP