[prev in list] [next in list] [prev in thread] [next in thread] 

List:       python-edu-sig
Subject:    [Edu-sig] Re: quadratic functions
From:       urner () alumni ! princeton ! edu (Kirby Urner)
Date:       2001-10-13 5:12:20
Message-ID: 18dfsto103oo1ghq01mcre325fk4me65r2 () 4ax ! com
[Download RAW message or body]

On Fri, 12 Oct 2001 12:08:13 -0700, I wrote:

>And is there a more general expression for
>
>               n
>   f(n,c) =  SIGMA k^c = some polynomial of degree c+1 ?
>             k = 1
>
>
>The answer to the 2nd question is yes, and the Bernoulli numbers
>give us the coefficients of such a c+1 degree polynomial.  

Just to spell that out at the command line, what we're 
looking for is a closed form expression for the sum of
consecutive integers all to the cth power.  What we 
expect is some polynomial of degree (c+1) with rational
coefficients in which the Bernoulli numbers figure (as
fractions).

So, for example, what is 1^4 + 2^4 + 3^4 + 4^4... n^4?

An easy function just does the arithmetic in a loop:

  >>> from operator import add
  >>> def sum(n,c):
	  return reduce(add,[i**c for i in range(1,n+1)])

Note that range(1,n) generates the list [1,2...n-1], hence
the parameter n+1.  [i**c for i in range(1,n+1)] is the 
sequence of terms in list form, and reduce(add, list) 
simply sums them up.

  >>> sum(10,4)
  25333

That's the sum of the first 10 natural numbers to the 4th power.

Bernoulli got his numbers from Johann Faulhaber, who wrote 
about this series in 1631.[1]  With this series, Bernoulli 
could brag about how quickly he could get the sum of the 10th
powers of the first 1000 natural numbers.  

Using our sum function:

  >>> sum(1000,10)
  91409924241424243424241924242500L

(the L means 'long integer' -- it'll be dropped in later 
versions of Python, as the integer types become more seamless).

But Bernoulli was using a polynomial:

  >>> summa(10)
  (1/11)*x**11 + (1/2)*x**10 + (5/6)*x**9 - x**7 + x**5 - (1/2)*x**3 + (5/66)*x

The coefficients above aren't naked Bernoulli numbers; they're 
folded into products which also involve choice(n+1,k), where
k is the kth term.  

  kth coefficient = b[k]*choice(n+1,k)   (b[k] = kth bernoulli)
                    ------------------
                         n+1

Given x=1000, the powers are easy to compute.  summa(10) 
returns a polynomial object into which we can substitute a 
value using p() syntax i.e.:

  >>> p = summa(10)
  >>> p(1000)
  91409924241424243424241924242500

Let's look at the first 10 bernoullis:

  >>> bernseries(10)
  [1, (1/2), (1/6), 0, (-1/30), 0, (1/42), 0, (-1/30), 0, (5/66)]

Not so bad.  The odd ones become 0 after the 0th, 1st, 2nd.

The Bernoulli numbers get pretty crazy (see below) and difficult
to compute by hand after the first few.  But our computer 
algorithm, using the Fraction object, is up to the task. 

Let's ask for the 100th Bernoulli number:

  >>> bern(100)
  (-94598037819122125295227433069493721872702841533066936
  133385696204311395415197247711/33330)

Wow, that took almost 10 seconds.  And it agrees with the result 
posted here (yay):

ftp://ftp.cdrom.com/pub/gutenberg/etext01/brnll10.txt

We calculate successive bernoullis using Pascal's Triangle.
The bernseries algorithm will show how this works if passed 
a non-zero "show" parameter (2nd argument):

  >>> bernseries(4,1)
  b[1] = (1/2) * (1)
  b[2] = (1/3) * (3*b[1]- 1)
  b[3] = (1/4) * (6*b[2] - 4*b[1]+ 1)
  b[4] = (1/5) * (10*b[3] - 10*b[2] + 5*b[1]- 1)
  [1, (1/2), (1/6), 0, (-1/30)]

Those'd be the bernoullis used to compute the sum of 
consecutive integers to the 4th power.  The corresponding
5th degree polynomial:

  >>> summa(4)
  (1/5)*x**5 + (1/2)*x**4 + (1/3)*x**3 - (1/30)*x

And the result for x = 10?

  >>> p = summa(4)
  >>> p(10)
  25333

That agrees with our earlier answer.
                                          n
And note we get the earlier formula for SIGMA k^2
                                        k = 1

>Well, you have metal cluster experiments, written up in 
>Scientific American (Dec 1998), about agglomerations of atoms
                          ^^^^
                          1997

>in icosahedral form transforming into cuboctahedra.[1]  How many
>atoms in all?  These formula will tell you.

So how do we do the total number of spheres in a cuboctahedron
again?  I posted quite a bit about this over on math-teach 
awhile ago.

For each layer, you've got the six squares and 8 triangles, 
but then you have shared edges and vertices, so adding all 
the squares and triangles involves double-counting edge 
spheres and quadruple counting the vertices.  The number of
edge spheres is 24 * (n-2) and the number of vertices is 12
(see http://www.inetarena.com/~pdx4d/ocn/urner.html for a
graphic).  So 8*n(n+1)/2 + 4*n^2 - 24*(n-2) - 3*12 should
be it.  That simplifies to:  10*n^2 - 20*n + 12.  This 
works for n>=2, i.e. a cuboctahedron with 2 spheres along
an edge has 12 spheres in that layer, with 3 spheres along
an edge has 42, then 92, then 162 and so on (icosahedral
shells the same, as per the jitterbug transformation:
http://www.inetarena.com/~pdx4d/ocn/numeracy0.html ).

What we want is a cumulative number of spheres in all 
layers.  In other words:

        n
  1 + SIGMA (10k^2 - 20k + 12)
       k=2

the 1 being the nuclear sphere.

       n
  10 SIGMA k^2 = 10 [(1/3)*n^3 + (1/2)*n^2 + (1/6)*n) - 1]
      k=2

       n
  20 SIGMA k =  20 [n(n+1)/2 - 1]
      k=2

       n
  12 SIGMA 1 = 12[n - 1]
      k=2

Combining and simplifying gives:

  Cubocta(n) = (10/3)*n^3 - 5*n^2 + (11/3)*n - 1

(remember to add 1 in front of all the SIGMA terms).

Making this a function:

  >>> def cubocta(n):
  	  return int((10/3)*n**3 - 5*n**2 + (11/3)*n - 1)

  >>> map(cubocta,range(1,11))
  [1, 13, 55, 147, 309, 561, 923, 1415, 2057, 2869]

That's the first 10 cuboctahedral numbers.

Citation:
http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=A005902

=========

ID Number: A005902 (Formerly M4898)
Sequence:  1,13,55,147,309,561,923,1415,2057,2869,3871,5083,6525,8217,
           10179,12431,14993,17885,21127,24739,28741,33153,37995,43287,
           49049,55301,62063,69355,77197,85609,94611,104223,114465,
           125357,136919,149171
Name:      Centered icosahedral (or cuboctahedral) numbers, also crystal ball
              sequence for f.c.c. lattice.
Comments:  Called "magic numbers" in some chemical contexts.
References S. Bjornholm, Clusters..., Contemp. Phys. 31 1990 pp. 309-324.
           T. P. Martin, Shells of atoms, Phys. Reports, 273 (1996), 199-241,
eq.
              (2).
           B. K. Teo and N. J. A. Sloane, Magic numbers in polygonal and
polyhedral
              clusters, Inorgan. Chem. 24 (1985), 4545-4558.
Links:     J. H. Conway and N. J. A. Sloane, Low-Dimensional Lattices VII:
Coordination Sequences, Proc. Royal Soc. London, A453 (1997), 2369-2389 (pdf,
ps).
Formula:   (2*n+1)*(5*n^2+5*n+3)/3.
Maple:     A005902:= n -> (2*n+1)*(5*n^2+5*n+3)/3;

=========

Kirby

PS:  Source code for bernoulli.py is at the Oregon Tips website,
along with the mathobjects package upon which it depends:
http://www.oregon-tips.org/

[1]  See John H. Conway, Richard K. Guy, 'The Book of Numbers' 
(Springer-Verlag, 1996) p. 106-109 re Faulhaber and Bernoulli 
numbers.



[prev in list] [next in list] [prev in thread] [next in thread]