The mystery of Bernoulli polynomials
Stefan Güttel, 8th February 2012
If there is another class of polynomials that is as fascinating and important to mathematics as orthogonal polynomials, then these are probably the Bernoulli polynomials B_j(x), deg(B_j) = j. These polynomials appear in the most different areas of mathematics, and have a variety of applications. In this example we have cited from the excellent Wikipedia articles  and , and a talk of Karl Dilcher .
Bernoulli polynomials are typically defined on the interval [0,1]. They can be generated recursively by integrating and adding a constant such that the definite integral equals zero. Let us build a quasimatrix whose (j+1)-st column is B_j(x), and plot the first 13 polynomials:
close all; clear all LW = 'linewidth'; lw = 2; format short x = chebfun('x',[0,1]); B(:,1) = 0*x + 1; for j = 1:100, B(:,j+1) = j*cumsum(B(:,j)); B(:,j+1) = B(:,j+1) - sum(B(:,j+1)); end plot(B(:,1:13),LW,lw) axis([0,1,-.3,.3])
The function values B_j(0) are called Bernoulli numbers. Here are the first 14 Bernoulli numbers:
ans = Columns 1 through 7 1.0000 -0.5000 0.1667 0.0000 -0.0333 -0.0000 0.0238 Columns 8 through 14 0.0000 -0.0333 0.0000 0.0758 -0.0000 -0.2531 -0.0000
These numbers turn out to be the Taylor coefficients of z/(exp(z)-1). Moreover, they are related to certain function values of the famous Riemann zeta function at integer arguments. In fact, the unresolved Riemann Hypothesis has an alternative reformulation due to Marcel Riesz (1916) in terms of Bernoulli numbers! Matlab doesn't come with a ZETA function, but if you have one available (see , for example), we can verify that the function values f(j), j = 0,...,13, coincide with the above Bernoulli numbers (for j = 1 the sign is switched):
if exist('zeta','file') f = chebfun(@(x) -x.*zeta(1-x),[0,13]); plot(f,LW,lw); hold on j = 0:13; f(j) plot(j,f(j),'ro',LW,lw); axis([0,13,-.4,1.1]); hold off end
ans = Columns 1 through 7 1.0000 0.5000 0.1667 0.0000 -0.0333 -0.0000 0.0238 Columns 8 through 14 -0.0000 -0.0333 0.0000 0.0758 -0.0000 -0.2531 0
Note that (except for j = 1) every second Bernoulli number is zero. These correspond to the trivial zeros of the Riemann zeta function. Using the above function f (and a generalization involving the so-called Hurwitz zeta function), one can define Bernoulli numbers (and polynomials) of non-integer index.
Bernoulli polynomials have the property that the number of (distinct) roots in the interval [0,1] is at most 3. We can easily verify this assertion numerically:
for j = 1:100, nrRoots(1,j) = length(roots(B(:,j))); end nrRoots(1:14) fprintf('The maximal number of roots is %d.\n',max(nrRoots))
ans = Columns 1 through 13 0 1 2 3 2 3 2 3 2 3 2 3 2 Column 14 3 The maximal number of roots is 3.
The multiplicity and location of complex roots has been of interest to many mathematicians over a long time. It is now known that all roots of the Bernoulli polynomials are distinct (Brillhart 1969, Dilcher 2008). It is also known that there exists a parabolic region above and below the interval [0,1] which is free of roots (Dilcher 1983/88):
figure for j = 1:20, r = roots(B(:,j),'all'); plot(r+1i*eps,'b*') hold on end hold off
Another interesting observation is the following: If one appropriately scales the even/odd Bernoulli polynomials, then these converge to cosine/sine functions, respectively. Let us visualize the first 50 rescaled Bernoulli polynomials of odd degree, and compute the distance to the expected limit function in the uniform norm:
err = ; limit = sin(2*pi*x); for j = 1:50, b = (-1)^(j)*(2*pi)^(2*j-1)/2/factorial(2*j-1)*B(:,2*j); plot(b) err(2*j) = norm(b - limit,inf); hold on end hold off
It is known that this convergence is geometric with rate 0.5:
semilogy(err,'b*','MarkerSize',10) hold on semilogy(0.5.^(0:99),'r--',LW,lw) hold off axis([0,100,1e-16,1])
The last property we like to mention and visualize here is the behavior of extrema of Bernoulli polynomials on [0,1]. D. H. Lehmer (1940) showed that the j-th degree Bernoulli polynomial is bounded by
for j > 1, except when j is 2 modulo 4, in which case the bound becomes
again with the Riemann zeta function.
fact = cumprod([1,1:99]); bound = 2*fact./(2*pi).^(0:99); for j = 1:100, M(j) = max(B(:,j)); if mod(j-1,4) == 2 && exist('zeta','file') bound(j) = bound(j)*zeta(j-1); end end semilogy(M,LW,lw); hold on semilogy(bound,'r--',LW,lw) axis([0,100,1e-5,1e80])
This bound looks quite sharp, and in fact, if one would remove the zeta(j) factor from the bound then it would be invalid every 4-th index.
 Wikipedia article on Bernoulli polynomials as of 08/01/2012, http://en.wikipedia.org/wiki/Bernoulli_polynomials
 Wikipedia article on Bernoulli numbers as of 08/01/2012, http://en.wikipedia.org/wiki/Bernoulli_numbers
 K. Dilcher, On Multiple Zeros of Bernoulli Polynomials, Talk at the 2011 "Special Functions in the 21st Century" conference in Washington, http://math.nist.gov/~DLozier/SF21/SF21slides/Dilcher.pdf
 Paul Godfrey, Special Functions math library, http://www.mathworks.com/matlabcentral/fileexchange/978