%% The mystery of Bernoulli polynomials
% Stefan Güttel, 8th February 2012
%%
% (Chebfun example roots/BernoulliPolynomials.m)
% [Tags: #roots, #Bernoulli, #Reimann]
%%
% 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 [1] and [2], and a talk of Karl Dilcher [3].
%
% 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:
B(0,1:14)
%%
% 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 [4], 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
%%
% 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))
%%
% 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
%
% 2*factorial(j)/(2*pi)^j
%
% for j > 1, except when j is 2 modulo 4, in which case the bound becomes
%
% 2*zeta(j)*factorial(j)/(2*pi)^j,
%
% 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.
%%
% References:
%
% [1] Wikipedia article on Bernoulli polynomials as of 08/01/2012,
% http://en.wikipedia.org/wiki/Bernoulli_polynomials
%
% [2] Wikipedia article on Bernoulli numbers as of 08/01/2012,
% http://en.wikipedia.org/wiki/Bernoulli_numbers
%
% [3] 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
%
% [4] Paul Godfrey, Special Functions math library,
% http://www.mathworks.com/matlabcentral/fileexchange/978