Bernoulli Polynomials.

This demonstration shows us computing the first few Bernoulli polynomials with Mathematica. We recall that the Bernoulli polynomials were defined by

((t((e^xt)- 1))/((e^t)- 1))= (∑_(j=1)^∞)(Subscript[B, j])(x)  ((t^j)/(j!))

which is to say that we're supposed to take a series expansion for the left-hand side in t and copy down the coefficients (in x). Ouch. 

We try it with Mathematica, though:

f := t(Exp[x*t] - 1) / (Exp[t] - 1);
sForm = Series[f,{t,0,8}]

Table[Expand[SeriesCoefficient[sForm,i]*(i!)],{i,1,8}] // TraditionalForm

No problem! Now the Bernoulli Numbers are defined to be the coefficients 

(t/((e^t)- 1)) = (∑_(j=0)^∞)(Subscript[B, j])  ((t^j)/(j!))
again, this looks like a pain in the neck to compute. However, we get:

sNForm = Series[t/(Exp[t] - 1),{t,0,8}]

Table[SeriesCoefficient[sNForm,i]*(i!),{i,0,8}]

and we notice the repeating pattern that all of the odd Bernoulli numbers >= 3 vanish. 

To get a relationship between the Bernoulli numbers and polynomials, observe that when we integrate

Integrate[ t(Exp[x*t] - 1) / (Exp[t] - 1),{x,0,1}]

which means (by term-by-term integration) that the Bernoulli numbers are (-) the integrals of the Bernoulli polynomials.