Bernoulli polynomials.

We defined the Bernoulli polynomials by the relation

	((t((e^xt)-1))/((e^t)-1))= (∑_(j=1)^∞)(Subscript[B, j])(x) ((t^j)/(j!))
	
To actually compute them, we need to expand the left-hand side as a power series (in t) and then collect coefficients.

This is kind of a mess to do by hand, but Mathematica makes this kind of computation relatively easy.

BernoulliComputation[j_] := FullSimplify[j! SeriesCoefficient[t (((E^(x t)) - 1)/((E^t)- 1)),{t,0,j}]]

As you’d expect “SeriesCoefficient” gives you the coefficient of (t^j) in the series expansion (around 0) of the function.

We can then make a table of Bernoulli polynomials:

Table[{(Subscript[B, j])[x],BernoulliComputation[j]},{j,0,6}] // TableForm

We defined the Bernoulli numbers (Subscript[B, j]) by

	(t/((e^t)-1))=∑ (Subscript[B, j]) ((t^j)/(j!))
	
and it’s not hard to get lists of these using the same SeriesCoefficient function

BernoulliNumber[j_] := j! SeriesCoefficient[(t/((E^t)-1)),{t,0,j}]

We can make a table of Bernoulli numbers the same way:

Table[{(Subscript[B, j]),BernoulliNumber[j]},{j,0,11}]//TableForm

This displays the pattern that we mentioned above: all the odd Bernoulli numbers (Subscript[B, 2j+1]) for j ≥ 1 vanish.

The relationship between (Subscript[B, j])(x) and (Subscript[B, j])

We can see something interesting if we integrate the Bernoulli polynomials:

BernoulliIntegral[j_] := Integrate[BernoulliComputation[j],{x,0,1}]

Table[{j,BernoulliIntegral[j]},{j,0,11}] // TableForm

Except for the 0-th term, these are equal to minus the Bernoulli numbers!