I want to understand the Euler-Maclaurin sum formula. On page 136 of A First Course in Numerical Analysis by Anthony Ralston and Philip Rabinowitz there is a derivation of the Euler-Maclaurin sum formula. They leave out the details and define non-standard Bernoulli polynomials. My derivation will fill in the details and use standard Bernoulli polynomials.

Our interest is in approximating the sum

n å j=0  f (x0 + jh)

where n may be infinite.

We begin by introducing the polynomials S_k(x), defined to be the polynomials of degree k

Sk(x) = k x-1 å p=1  pk-1

Using DERIVE we can show that

Sk(1)=0

Sk(0)=0        k > 1

The DERIVE results are:

S(k,x):=k*SUM(p^(k-1),p,1,x-1)

k:epsilonInteger

;User=Simp(User) S(k,1)=0

k:epsilonInteger (1, inf)

;User=Simp(User) S(k,0)=0

Using DERIVE we can find the first few polynomials:

S0(x)=0

S1(x)=x-1

S2(x)=x2-x

S3(x)=x3-  3 x2 2 +  x 2

S(k,x):=k*SUM(p^(k-1),p,1,x-1)

;Expd(User') VECTOR(S(k,x),k,0,3)=[0,x-1,x^2-x,x^3-3*x^2/2+x/2]

The constants B_k are the Bernoulli numbers. The following two idenities are of particular interest to us:

B2 k+1=0        k > 0

d Sk(x) dx =k (Sk-1(x)+Bk-1)       k > 2

Using DERIVE we can verify these idenities for several values of k.

B(k):=IF(k=0,1,IF(k=1,-1/2,-k*ZETA(1-k)))

;User=Simp(User) VECTOR(B(2*k+1),k,1,12)=[0,0,0,0,0,0,0,0,0,0,0,0]

S(k,x):=k*SUM(p^(k-1),p,1,x-1)

;User=Simp(User)

VECTOR(DIF(S(k,x),x)-k*(S(k-1,x)+B(k-1)),k,3,12)=[0,0,0,0,0,0,0,0,0,0]

Now we define

Mk(x,h)=  1 (2 k)! ó õ x+h x  S2 k æ è  y-x h ö ø  f (2 k)(y) dy

Using DERIVE we find M₁(x,h).

[F(x):=, S(k,x):=k*SUM(p^(k-1),p,1,x-1)]

M(k,x,h):=INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)/(2*k)!

For k=1

M1(x,h)=  1 2 ó õ x+h x  S2 æ è  y-x h ö ø  f (2)(y) dy

;Simp(User') S(2,(y-x)/h)=(x-y)*(x-y+h)/h^2

Evaluate S.

S2 æ è  y-x h ö ø =  (x-y) (x-y+h) h2

M(1,x,h)=INT((x-y)*(x-y+h)*F''(y),y,x,x+h)/(2*h^2)

Substute that value of S.

M1(x,h)= ó õ x+h x  (x-y) (x-y+h) f (2)(y) dy 2 h2

F(x):=

INT((x-y)*(x-y+h)*F''(y),y,x,x+h)/(2*h^2)

[V(y):=(y-x)*(y-x-h), U(y):=DIF(F(y),y,1)]

INT(V(y)*DIF(U(y),y),y,x,x+h)=

V(x+h)*U(x+h)-V(x)*U(x)-INT(U(y)*DIF(V(y),y),y,x,x+h)

Integrate M₁ by parts.

V(y)=(y-x) (y-x-h)

U(y)=f (1)(y)

ó õ x+h x  V(y)   d U(y) dy dy = V(x+h)*U(x+h)-V(x)*U(x)- ó õ x+h x  U(y)   d V(y) dy dy

;Simp(User') INT((x-y)*(x-y+h)*F''(y),y,x,x+h)=

(2*x+h)*F(x+h)-(2*x+h)*F(x)-2*INT(y*F'(y),y,x,x+h)

Evaluate the integral.

ó õ x+h x  (x-y) (x-y+h) f (2)(y) dy = (2 x+h) (f (x+h)-f (x))-2  ó õ x+h x  y f (1)(y) dy

[F(x):=,S(k,x):=k*SUM(p^(k-1),p,1,x-1)]

[U(y):=F(y), V(y):=y]

INT(V(y)*DIF(U(y),y),y,x,x+h)=

V(x+h)*U(x+h)-V(x)*U(x)-INT(U(y)*DIF(V(y),y),y,x,x+h)

Integrate by parts again.

V(y)=y

U(y)=f (y)

;Simp(#3) INT(y*F'(y),y,x,x+h)=(x+h)*F(x+h)-x*F(x)-INT(F(y),y,x,x+h)

Evaluate the integral.

ó õ x+h x  y f (1)(y) dy = (x+h) f (x+h)-x f (x)- ó õ x+h x  f (y) dy

;Simp(#6') INT((x-y)*(x-y+h)*F''(y),y,x,x+h)=

-h*F(x+h)-h*F(x)+2*INT(F(y),y,x,x+h)

Subsitute the last integral.

ó õ x+h x  (x-y) (x-y+h) f (2)(y) dy = -h (f (x+h)-f (x))+2  ó õ x+h x  f (y) dy

;Expd(#10') M(1,x,h)=-F(x+h)/(2*h)-F(x)/(2*h)+INT(F(y),y,x,x+h)/h^2

Subsitute this result to get the value of M₁.

M1(x,h)=-  f (x+h)+f (x) 2 h + ó õ x+h x  f (y) dy h2

Now to find the other values of M_k we will find a recurrence relationship. Start again with the definition of M_k.

Mk(x,h)=  1 (2 k)! ó õ x+h x  S2 k æ è  y-x h ö ø  f (2 k)(y) dy

Do the same thing in DERIVE.

[F(x):=,S(k,x):=]

M(k,x,h):=INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)/(2*k)!

Integrate M_k by parts.

V(y)=S2*k æ è  y-x h ö ø

U(y)=f (2 k-1)(y)

ó õ x+h x  V(y)   d U(y) dy dy = V(x+h)*U(x+h)-V(x)*U(x)- ó õ x+h x  U(y)   d V(y) dy dy

[V(y):=, U(y):=]

INT(V(y)*DIF(U(y),y),y,x,x+h)=

LIM(V(y)*U(y),y,x+h,0)-LIM(V(y)*U(y),y,x,0)-

INT(U(y)*DIF(V(y),y),y,x,x+h)

[V(y):=S(2*k,(y-x)/h), U(y):=DIF(F(y),y,2*k-1)]

;Simp(#4) INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)=

-INT(DIF(F(y),y,2*k-1)*LIM(DIF(S(2*k,@),@),@,(y-x)/h,0),y,x,x+h)/h+

LIM(S(2*k,(y-x)/h)*DIF(F(y),y,2*k-1),y,x+h,0)-

LIM(S(2*k,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,0)

[U(y)=DIF(F(y),y,2*k-1),U(y):=]

;Simp(User) INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)=

-INT(U(y)*LIM(DIF(S(2*k,@),@),@,(y-x)/h,0),y,x,x+h)/h+

U(x+h)*S(2*k,1)-U(x)*S(2*k,0)

[S(k,x):=k*SUM(p^(k-1),p,1,x-1),k:epsilonInteger (1, inf)]

;User=Simp(User) [S(2*k,1),S(2*k,0)]=[0,0]

[S(k,x):=,B(k):=]

;Simp(User) INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)=

-INT(U(y)*LIM(DIF(S(2*k,@),@),@,(y-x)/h,0),y,x,x+h)/h

LIM(DIF(S(2*k,@),@),@,(y-x)/h,0)=LIM(DIF(S(2*k,t),t),t,(y-x)/h,0)

[DIF(S(k,x),x)=k*(S(k-1,x)+B(k-1)),k>2]

;Sub(User) DIF(S(2*k,t),t)=(2*k)*(S(2*k-1,t)+B(2*k-1))

;Sub(#13') LIM(DIF(S(2*k,@),@),@,(y-x)/h,0)=

LIM((2*k)*(S(2*k-1,t)+B(2*k-1)),t,(y-x)/h,0)

B(k):=IF(k=0,1,IF(k=1,-1/2,-k*ZETA(1-k)))

;Simp(#16') LIM(DIF(S(2*k,@),@),@,(y-x)/h,0)=2*k*S(2*k-1,(y-x)/h)

;Simp(User) INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)=

-2*k*INT(U(y)*S(2*k-1,(y-x)/h),y,x,x+h)/h

U(y):=DIF(F(y),y,2*k-1)

;Simp(#19) INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)=

-2*k*INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)/h

Evaluate the integral.

ó õ x+h x  S2 k æ è  y-x h ö ø  f (2 k)(y) dy=

S2 k(1) f (2 k-1)(x+h)-S2 k(0) f (2 k-1)(x)- ó õ x+h x  f (2 k-1)(y)  d S2 k æ è  y-x h ö ø dy dy

This can be simplified using the two idenities. For k > 1.

S2 k(1)=0

S2 k(0)=0

d Sk(x) dx =k (Sk-1(x)+Bk-1)       k > 2

d S2 k æ è  y-x h ö ø dy = 2 k  æ è S2 k-1 æ è  y-x h ö ø +B2 k-1 ö ø h

We also know B_2 k-1=0 so it simplifys even more.

d S2 k æ è  y-x h ö ø dy = 2 k S2 k-1 æ è  y-x h ö ø h

Substitute that result into the integral.

ó õ x+h x  S2 k æ è  y-x h ö ø  f (2 k)(y) dy = - 2 k  ó õ x+h x  f (2 k-1)(y) S2 k-1 æ è  y-x h ö ø dy h

;Simp(User') M(k,x,h)=

-INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)/(h*(2*k-1)!)

Substitute the integral in the definition of M_k.

Mk(x,h)=- ó õ x+h x  f (2 k-1)(y) S2 k-1 æ è  y-x h ö ø dy h (2 k-1)!

[F(x):=,S(k,x):=]

M(k,x,h):=INT(S(2*k,(y-x)/h)*DIF(F(y),y,2*k),y,x,x+h)/(2*k)!

M(k,x,h)=-INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)/(h*(2*k-1)!)

[V(y):=S(2*k-1,(y-x)/h), U(y):=DIF(F(y),y,2*k-2)]

;Simp(#5) INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)=

-INT(DIF(F(y),y,2*k-2)*LIM(DIF(S(2*k-1,@),@),@,(y-x)/h,0),y,x,x+h)/h+

LIM(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-2),y,x+h,0)-

LIM(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-2),y,x,0)

[U(y)=DIF(F(y),y,2*k-2),U(y):=]

;Simp(User) INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)=

-INT(U(y)*LIM(DIF(S(2*k-1,@),@),@,(y-x)/h,0),y,x,x+h)/h+

U(x+h)*S(2*k-1,1)-U(x)*S(2*k-1,0)

[S(k,x):=k*SUM(p^(k-1),p,1,x-1),k:epsilonInteger (1, inf)]

;User=Simp(User) [S(2*k-1,0),S(2*k-1,1)]=[0,0]

[S(k,x):=,B(k):=]

;Simp(User) INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)=

-INT(U(y)*LIM(DIF(S(2*k-1,@),@),@,(y-x)/h,0),y,x,x+h)/h

LIM(DIF(S(2*k-1,@),@),@,(y-x)/h,0)=LIM(DIF(S(2*k-1,t),t),t,(y-x)/h,0)

[DIF(S(k,x),x)=k*(S(k-1,x)+B(k-1)),k>2]

;Sub(User) DIF(S(2*k-1,t),t)=(2*k-1)*(B(2*k-2)+S(2*k-2,t))

;Simp(User) LIM(DIF(S(2*k-1,@),@),@,(y-x)/h,0)=

(2*k-1)*(B(2*k-2)+S(2*k-2,(y-x)/h))

;Simp(User) INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)=

(1-2*k)*(B(2*k-2)*INT(U(y),y,x,x+h)+

INT(U(y)*S(2*k-2,(y-x)/h),y,x,x+h))/h

U(y):=DIF(F(y),y,2*k-2)

;Simp(User) INT(S(2*k-1,(y-x)/h)*DIF(F(y),y,2*k-1),y,x,x+h)=

(1-2*k)*(B(2*k-2)*(LIM(DIF(F(y),y,2*k-3),y,x+h,0)-

LIM(DIF(F(y),y,2*k-3),y,x,0))+

INT(S(2*k-2,(y-x)/h)*DIF(F(y),y,2*k-2),y,x,x+h))/h

Integrate M_k by parts again.

V(y)=S2*k-1 æ è  y-x h ö ø

U(y)=f (2 k-2)(y)

ó õ x+h x  V(y)   d U(y) dy dy = V(x+h)*U(x+h)-V(x)*U(x)- ó õ x+h x  U(y)   d V(y) dy dy

Evaluate the integral.

ó õ x+h x  S2 k-1 æ è  y-x h ö ø  f (2 k-1)(y) dy=

S2 k-1(1) f (2 k-2)(x+h)-S2 k-1(0) f (2 k-2)(x)- ó õ x+h x  f (2 k-2)(y)  d S2 k-1 æ è  y-x h ö ø dy dy

This can be simplified using the two idenities. For k > 1.

S2 k-1(1)=0

S2 k-1(0)=0

d Sk(x) dx =k (Sk-1(x)+Bk-1)       k > 2

d S2 k-1 æ è  y-x h ö ø dy = (2 k-1)  æ è S2 k-2 æ è  y-x h ö ø +B2 k-2 ö ø h

Substitute that result into the integral.

ó õ x+h x  S2 k-1 æ è  y-x h ö ø  f (2 k-1)(y) dy=

- (2 k-1)  ó õ x+h x  f (2 k-2)(y)  æ è S2 k-2 æ è  y-x h ö ø +B2 k-2 ö ø dy h

;Simp(User') M(k,x,h)=

(INT(S(2*k-2,(y-x)/h)*DIF(F(y),y,2*k-2),y,x,x+h)-

B(2*k-2)*(LIM(DIF(F(y),y,2*k-3),y,x+h,0)-

LIM(DIF(F(y),y,2*k-3),y,x,0)))/(h^2*(2*k-2)!)

;User=Simp(User) M(k-1,x,h)=

INT(S(2*(k-1),(y-x)/h)*DIF(F(y),y,2*(k-1)),y,x,x+h)/(2*k-2)!

;User=Simp(User) M(k,x,h):=

M(k-1,x,h)*(2*k-2)!=

INT(S(2*(k-1),(y-x)/h)*DIF(F(y),y,2*(k-1)),y,x,x+h)

;Simp(User) M(k,x,h)=B(2*k-2)*(LIM(DIF(F(y),y,2*k-3),y,x+h,0)-

LIM(DIF(F(y),y,2*k-3),y,x,0))/(h^2*(2*k-2)!)+M(k-1,x,h)/h^2

Substitute the integral in the definition of M_k.

Mk(x,h)= ó õ x+h x  f (2 k-2)(y)  æ è S2 k-2 æ è  y-x h ö ø +B2 k-2 ö ø dy h2 (2 k-2)!

Expand the integral.

Mk(x,h)= ó õ x+h x  f (2 k-2)(y) S2 k-2 æ è  y-x h ö ø dy h2 (2 k-2)! +  B2 k-2 (f (2 k-3)(x+h)-f (2 k-3)(x)) h2 (2 k-2)!

Form the definition we can see M_k-1.

Mk-1(x,h)=  1 (2 (k-1))! ó õ x+h x  S2 (k-1) æ è  y-x h ö ø  f (2 (k-1))(y) dy

We substitute M_k-1 in the recurrence relation for M_k.

Mk(x,h)=  Mk-1(x,h) h2 +  B2 k-2 (f (2 k-3)(x+h)-f (2 k-3)(x)) h2 (2 k-2)!

Our interest is in approximating the sum

n å j=0  f (x0 + jh)

We really do not care what M_k is. We know what M₁ is.

M1(x,h)=-  f (x+h)+f (x) 2 h + ó õ x+h x  f (y) dy h2

Reorder the equation to be more useful.

f (x+h)+f (x) 2 = ó õ x+h x  f (y) dy h -h M1(x,h)

Reorder the recurrence relation for M_k.

Mk-1(x,h)=-  B2 k-2 (f (2 k-3)(x+h)-f (2 k-3)(x)) (2 k-2)! +h2 Mk(x,h)

We can get M₁ from M₂ from this recurrence relation with k=2.

M1(x,h)=-  B2 (f (1)(x+h)-f (1)(x)) 2! +h2 M2(x,h)

We can get M₁ from M₃ from the recurrence relation.

M2(x,h)=-  B4 (f (3)(x+h)-f (3)(x)) 4! +h2 M3(x,h)

M1(x,h)=-  B2 (f (1)(x+h)-f (1)(x)) 2! -h2   B4 (f (3)(x+h)-f (3)(x)) 4! +h4 M3(x,h)

We can get M₁ from M_l+1 from the recurrence relation.

M1(x,h)=- l å k=1   h2 k-2 B2 k (f (2 k-1)(x+h)-f (2 k-1)(x)) (2 k)! +h2 l Ml+1(x,h)

Substitute this for M₁ in a previous equation.

f (x+h)+f (x) 2 = ó õ x+h x  f (y) dy h +

l å k=1   h2 k-1 B2 k (f (2 k-1)(x+h)-f (2 k-1)(x)) (2 k)! -h2 l+1 Ml+1(x,h)

Now we can verify some examples using DERIVE.

[F(x):=, B(k):=]

;Simp(User') M(k,x,h):=

IF(k=1,-F(x+h)/(2*h)-F(x)/(2*h)+INT(F(y),y,x,x+h)/h^2,

B(2*k-2)*(LIM(DIF(F(y),y,2*k-3),y,x+h,0)-

LIM(DIF(F(y),y,2*k-3),y,x,0))/(h^2*(2*k-2)!)+M(k-1,x,h)/h^2)

(F(x+h)+F(x))/2=INT(F(y),y,x,x+h)/h+

SUM(h^(2*k-1)*B(2*k)*(LIM(DIF(F(y),y,2*k-1),y,x+h,0)-

LIM(DIF(F(y),y,2*k-1),y,x,0))/(2*k)!,k,1,l)-h^(2*l+1)*M(l+1,x,h)

l=12

;Simp(Sub(#3)) true

Now using derive we can find

n å j=0  f (x0 + jh)

F(x):=

B(k):=

M(k,x,h):=

;Sub(User) (F(x+h)+F(x))/2=INT(F(y),y,x,x+h)/h+

SUM(h^(2*k-1)*B(2*k)*(LIM(DIF(F(y),y,2*k-1),y,x+h,0)-

LIM(DIF(F(y),y,2*k-1),y,x,0))/(2*k)!,k,1,l)-h^(2*l+1)*M(l+1,x,h)

This is results on the previous page.

f (x+h)+f (x) 2 = ó õ x+h x  f (y) dy h +

l å k=1   h2 k-1 B2 k (f (2 k-1)(x+h)-f (2 k-1)(x)) (2 k)! -h2 l+1 Ml+1(x,h)

SUM(F(x0+j*h),j,0,n)=F(h*n+x0)/2+F(x0)/2+

SUM(F(h*j+x0),j,1,n)/2+SUM(F(h*j+x0),j,0,n-1)/2

Expand the sumation into four parts.

n å j=0  f (x0 + jh) =  f (x0) 2 + n-1 å j=0   f (x0 + jh) 2 + n å j=1   f (x0 + jh) 2 +  f (x0 + n h) 2

;Sub(User) SUM(F(h*j+x0),j,0,n)=F(h*n+x0)/2+F(x0)/2+

SUM((F(h*j+h+x0)+F(h*j+x0))/2,j,0,n-1)

Combine two of the parts.

n å j=0  f (x0 + jh) =  f (x0)+f (x0 + n h) 2 + n-1 å j=0   f (x0 + jh)+f (x0 + jh + h) 2

;Sub(#5) (F((x0+j*h)+h)+F(x0+j*h))/2=

INT(F(y),y,x0+j*h,(x0+j*h)+h)/h+

SUM(h^(2*k-1)*B(2*k)*(LIM(DIF(F(y),y,2*k-1),y,(x0+j*h)+h,0)-

LIM(DIF(F(y),y,2*k-1),y,x0+j*h,0))/(2*k)!,k,1,l)-

h^(2*l+1)*M(l+1,x0+j*h,h)

Substitute for the last part.

n å j=0  f (x0 + jh) =  f (x0)+f (x0 + n h) 2 + n-1 å j=0  ó õ x0+jh+h x0+jh  f (y) dy h +

n-1 å j=0  l å k=1   h2 k-1 B2 k (f (2 k-1)(x0+jh+h)-f (2 k-1)(x0+jh)) (2 k)! -h2 l+1 Ml+1(x0+jh,h)

;Sub(#7) SUM(F(h*j+x0),j,0,n)=F(h*n+x0)/2+

F(x0)/2+SUM(INT(F(y),y,x0+j*h,(x0+j*h)+h)/h+

SUM(h^(2*k-1)*B(2*k)*(LIM(DIF(F(y),y,2*k-1),y,(x0+j*h)+h,0)-

LIM(DIF(F(y),y,2*k-1),y,x0+j*h,0))/(2*k)!,k,1,l)-

h^(2*l+1)*M(l+1,x0+j*h,h),j,0,n-1)

;Simp(#9) SUM(F(h*j+x0),j,0,n)=F(h*n+x0)/2+F(x0)/2+

SUM(-h^(2*l+1)*M(l+1,h*j+x0,h)+INT(F(y),y,h*j+x0,h*j+h+x0)/h+

SUM(h^(2*k)*B(2*k)*(LIM(DIF(F(y),y,2*k-1),y,h*j+h+x0,0)-

LIM(DIF(F(y),y,2*k-1),y,h*j+x0,0))/(2*k)!,k,1,l)/h,j,0,n-1)

;Sub(#11) SUM(F(h*j+x0),j,0,n)=F(h*n+x0)/2+F(x0)/2+

(-h^(2*l+1)*SUM(M(l+1,h*j+x0,h),j,0,n-1)+ INT(F(y),y,x0,h*n+x0)/h+

SUM(h^(2*k)*B(2*k)*(LIM(DIF(F(y),y,2*k-1),y,h*n+x0,0)-

LIM(DIF(F(y),y,2*k-1),y,x0,0))/(2*k)!,k,1,l)/h)

Simplify the sums.

n å j=0  f (x0 + jh) =  f (x0)+f (x0 + n h) 2 + ó õ x0+n h x0  f (y) dy h +

l å k=1   h2 k-1 B2 k (f (2 k-1)(x0+n h)-f (2 k-1)(x0)) (2 k)! -h2 l+1  n-1 å j=0  Ml+1(x0+jh,h)