SC11-I-a

4227 days ago by SC2011

Computational and Experimental Mathematics

Introduction: Some Sums

Neil Calkin, Holly Hirst, and Dan Warner

November, 2011
 
       

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

Note

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 $j$ takes on the values 0 through $n-1$.
# Let np denote the length of the arithmetic progression, ap. np = len(ap) 
       
# Consider the first differences d1 = [(ap[j]-ap[j-1]) for j in range(1,np)] print 'The sequence of first differences' print d1 
       

Observation

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

$\alpha + \beta \,n$

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 s1 = copy(ap) for k in range(1,np): s1[k] += s1[k-1] print 'The sequence of sums' print s1 
       
d1 = [(s1[j]-s1[j-1]) for j in range(1,np)] print 'The sequence of first differences' print d1 d2 = [(d1[j]-d1[j-1]) for j in range(1,np-1)] print 'The sequence of second differences' print d2 
       

Observation

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 $s_1$, the sequence of sums, into the

       Online Encyclopedia of Integer Sequences 

http://oeis.org/classic/

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 with $n$ starting at 0. 

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

We see that this sequence also occurs in a number of other situations including the Balmer series of the Hydrogen atom.

 
print 'Adjusted OEIS formula' ap3 = [(s^2+6*s+5) for s in range(np)] print ap3 
       

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) print sol 
       
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] 
       
# Lazy Caterer 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 print R print R*F 
       
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/ # One of ISC's best guesses is A001622 the decimal expansion of the Golden Ratio phi = (1 + sqrt(5))/2 phin = phi.n(digits=30) fn = (F3[n-1]/F3[n-2]).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 
       
(1/phi).n(digits=100) 
       

How about the powers of $\phi$?
var('phi') Eq = phi^2 == phi + 1 print Eq Eq1 = Eq for k in range(10): Eq1 = expand(phi*Eq1) # print Eq1 Eq1 = Eq1.subs_expr(Eq) print Eq1 
       
# Here's a cute matrix identity PP = Pt * transpose(Pt) print PP