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

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})$

where $n$ may be infinite.

We begin by introducing the polynomials $S{}_k(x)$, defined to be the polynomials of degree $k$

$S{}_k(x)=k\sum _{p=1}^{x-1}{p}^{k-1}$

Using DERIVE we can show that

$S{}_k(1)=0$

$S{}_k(0)=0\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}k>1$

The DERIVE results are:

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

k:epsilonInteger

Using DERIVE we can find the first few polynomials:

$S{}_0(x)=0$

$S{}_1(x)=x-1$

$S{}_2(x)=x^2-x$

$S{}_3(x)=x^3-\frac{3\hspace{1em}{x}^2}{2}+\frac{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:

$B{}_{2\hspace{1em}k+1}=0\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}k>0$

$\frac{d\hspace{1em}S{}_k(x)}{\mathrm{dx}}=k\hspace{1em}(S{}_{k-1}(x)+B{}_{k-1})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}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]

Now we define

$M{}_k(x,h)=\frac{1}{(2\hspace{1em}k)!}{\int }_x^{x+h}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k)}(y)\hspace{1em}\mathrm{dy}$

Using DERIVE we find $M{}_1(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$

$M{}_1(x,h)=\frac{1}{2}{\int }_x^{x+h}S{}_2\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2)}(y)\hspace{1em}\mathrm{dy}$

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

Evaluate S.

$S{}_2\left(\frac{y-x}{h}\right)=\frac{(x-y)\hspace{1em}(x-y+h)}{h^2}$

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

Substute that value of S.

$M{}_1(x,h)=\frac{{\int }_x^{x+h}(x-y)\hspace{1em}(x-y+h)\hspace{1em}f\hspace{1em}{}^{(2)}(y)\hspace{1em}\mathrm{dy}}{2\hspace{1em}{h}^2}$

F(x):=

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

Integrate $M{}_1$ by parts.

$V(y)=(y-x)\hspace{1em}(y-x-h)$

$U(y)=f\hspace{1em}{}^{(1)}(y)$

${\int }_x^{x+h}V(y)\hspace{1em}\frac{d\hspace{1em}U(y)}{\mathrm{dy}}\mathrm{dy}=V(x+h)*U(x+h)-V(x)*U(x)-{\int }_x^{x+h}U(y)\hspace{1em}\frac{d\hspace{1em}V(y)}{\mathrm{dy}}\mathrm{dy}$

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

Evaluate the integral.

${\int }_x^{x+h}(x-y)\hspace{1em}(x-y+h)\hspace{1em}f\hspace{1em}{}^{(2)}(y)\hspace{1em}\mathrm{dy}=(2\hspace{1em}x+h)\hspace{1em}(f\hspace{1em}(x+h)-f\hspace{1em}(x))-2\hspace{1em}{\int }_x^{x+h}y\hspace{1em}f\hspace{1em}{}^{(1)}(y)\hspace{1em}\mathrm{dy}$

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

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

Integrate by parts again.

$V(y)=y$

$U(y)=f\hspace{1em}(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.

${\int }_x^{x+h}y\hspace{1em}f\hspace{1em}{}^{(1)}(y)\hspace{1em}\mathrm{dy}=(x+h)\hspace{1em}f\hspace{1em}(x+h)-x\hspace{1em}f\hspace{1em}(x)-{\int }_x^{x+h}f\hspace{1em}(y)\hspace{1em}\mathrm{dy}$

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

Subsitute the last integral.

${\int }_x^{x+h}(x-y)\hspace{1em}(x-y+h)\hspace{1em}f\hspace{1em}{}^{(2)}(y)\hspace{1em}\mathrm{dy}=-h\hspace{1em}(f\hspace{1em}(x+h)-f\hspace{1em}(x))+2\hspace{1em}{\int }_x^{x+h}f\hspace{1em}(y)\hspace{1em}\mathrm{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{}_1$.

$M{}_1(x,h)=-\frac{f\hspace{1em}(x+h)+f\hspace{1em}(x)}{2\hspace{1em}h}+\frac{{\int }_x^{x+h}f\hspace{1em}(y)\mathrm{dy}}{h^2}$

Now to find the other values of $M{}_k$ we will find a recurrence relationship. Start again with the definition of $M{}_k$.

$M{}_k(x,h)=\frac{1}{(2\hspace{1em}k)!}{\int }_x^{x+h}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k)}(y)\hspace{1em}\mathrm{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)=S{}_{2*k}\left(\frac{y-x}{h}\right)$

$U(y)=f\hspace{1em}{}^{(2\hspace{1em}k-1)}(y)$

${\int }_x^{x+h}V(y)\hspace{1em}\frac{d\hspace{1em}U(y)}{\mathrm{dy}}\mathrm{dy}=V(x+h)*U(x+h)-V(x)*U(x)-{\int }_x^{x+h}U(y)\hspace{1em}\frac{d\hspace{1em}V(y)}{\mathrm{dy}}\mathrm{dy}$

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

Evaluate the integral.

${\int }_x^{x+h}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k)}(y)\hspace{1em}\mathrm{dy}=$

$S{}_{2\hspace{1em}k}(1)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x+h)-S{}_{2\hspace{1em}k}(0)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x)-{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(y)\hspace{1em}\frac{d\hspace{1em}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)}{\mathrm{dy}}\mathrm{dy}$

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

$S{}_{2\hspace{1em}k}(1)=0$

$S{}_{2\hspace{1em}k}(0)=0$

$\frac{d\hspace{1em}S{}_k(x)}{\mathrm{dx}}=k\hspace{1em}(S{}_{k-1}(x)+B{}_{k-1})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}k>2$

$\frac{d\hspace{1em}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)}{\mathrm{dy}}=\frac{2\hspace{1em}k\hspace{1em}(S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)+B{}_{2\hspace{1em}k-1})}{h}$

We also know $B{}_{2\hspace{1em}k-1}=0$ so it simplifys even more.

$\frac{d\hspace{1em}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)}{\mathrm{dy}}=\frac{2\hspace{1em}k\hspace{1em}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)}{h}$

Substitute that result into the integral.

${\int }_x^{x+h}S{}_{2\hspace{1em}k}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k)}(y)\hspace{1em}\mathrm{dy}=-\frac{2\hspace{1em}k\hspace{1em}{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(y)\hspace{1em}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)\mathrm{dy}}{h}$

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

Substitute the integral in the definition of $M{}_k$.

$M{}_k(x,h)=-\frac{{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(y)\hspace{1em}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)\mathrm{dy}}{h\hspace{1em}(2\hspace{1em}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)!

Integrate $M{}_k$ by parts again.

$V(y)=S{}_{2*k-1}\left(\frac{y-x}{h}\right)$

$U(y)=f\hspace{1em}{}^{(2\hspace{1em}k-2)}(y)$

${\int }_x^{x+h}V(y)\hspace{1em}\frac{d\hspace{1em}U(y)}{\mathrm{dy}}\mathrm{dy}=V(x+h)*U(x+h)-V(x)*U(x)-{\int }_x^{x+h}U(y)\hspace{1em}\frac{d\hspace{1em}V(y)}{\mathrm{dy}}\mathrm{dy}$

Evaluate the integral.

${\int }_x^{x+h}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(y)\hspace{1em}\mathrm{dy}=$

$S{}_{2\hspace{1em}k-1}(1)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k-2)}(x+h)-S{}_{2\hspace{1em}k-1}(0)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k-2)}(x)-{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-2)}(y)\hspace{1em}\frac{d\hspace{1em}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)}{\mathrm{dy}}\mathrm{dy}$

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

$S{}_{2\hspace{1em}k-1}(1)=0$

$S{}_{2\hspace{1em}k-1}(0)=0$

$\frac{d\hspace{1em}S{}_k(x)}{\mathrm{dx}}=k\hspace{1em}(S{}_{k-1}(x)+B{}_{k-1})\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}k>2$

$\frac{d\hspace{1em}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)}{\mathrm{dy}}=\frac{(2\hspace{1em}k-1)\hspace{1em}(S{}_{2\hspace{1em}k-2}\left(\frac{y-x}{h}\right)+B{}_{2\hspace{1em}k-2})}{h}$

Substitute that result into the integral.

${\int }_x^{x+h}S{}_{2\hspace{1em}k-1}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}k-1)}(y)\hspace{1em}\mathrm{dy}=$

$-\frac{(2\hspace{1em}k-1)\hspace{1em}{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-2)}(y)\hspace{1em}(S{}_{2\hspace{1em}k-2}\left(\frac{y-x}{h}\right)+B{}_{2\hspace{1em}k-2})\mathrm{dy}}{h}$

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

Substitute the integral in the definition of $M{}_k$.

$M{}_k(x,h)=\frac{{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-2)}(y)\hspace{1em}(S{}_{2\hspace{1em}k-2}\left(\frac{y-x}{h}\right)+B{}_{2\hspace{1em}k-2})\mathrm{dy}}{h^2\hspace{1em}(2\hspace{1em}k-2)!}$

Expand the integral.

$M{}_k(x,h)=\frac{{\int }_x^{x+h}f\hspace{1em}{}^{(2\hspace{1em}k-2)}(y)\hspace{1em}S{}_{2\hspace{1em}k-2}\left(\frac{y-x}{h}\right)\mathrm{dy}}{h^2\hspace{1em}(2\hspace{1em}k-2)!}+\frac{B{}_{2\hspace{1em}k-2}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-3)}(x+h)-f\hspace{1em}{}^{(2\hspace{1em}k-3)}(x))}{h^2\hspace{1em}(2\hspace{1em}k-2)!}$

Form the definition we can see $M{}_{k-1}$.

$M{}_{k-1}(x,h)=\frac{1}{(2\hspace{1em}(k-1))!}{\int }_x^{x+h}S{}_{2\hspace{1em}(k-1)}\left(\frac{y-x}{h}\right)\hspace{1em}f\hspace{1em}{}^{(2\hspace{1em}(k-1))}(y)\hspace{1em}\mathrm{dy}$

We substitute $M{}_{k-1}$ in the recurrence relation for $M{}_k$.

$M{}_k(x,h)=\frac{M{}_{k-1}(x,h)}{h^2}+\frac{B{}_{2\hspace{1em}k-2}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-3)}(x+h)-f\hspace{1em}{}^{(2\hspace{1em}k-3)}(x))}{h^2\hspace{1em}(2\hspace{1em}k-2)!}$

Our interest is in approximating the sum

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})$

We really do not care what $M{}_k$ is. We know what $M{}_1$ is.

$M{}_1(x,h)=-\frac{f\hspace{1em}(x+h)+f\hspace{1em}(x)}{2\hspace{1em}h}+\frac{{\int }_x^{x+h}f\hspace{1em}(y)\mathrm{dy}}{h^2}$

Reorder the equation to be more useful.

$\frac{f\hspace{1em}(x+h)+f\hspace{1em}(x)}{2}=\frac{{\int }_x^{x+h}f\hspace{1em}(y)\mathrm{dy}}{h}-h\hspace{1em}M{}_1(x,h)$

Reorder the recurrence relation for $M{}_k$.

$M{}_{k-1}(x,h)=-\frac{B{}_{2\hspace{1em}k-2}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-3)}(x+h)-f\hspace{1em}{}^{(2\hspace{1em}k-3)}(x))}{(2\hspace{1em}k-2)!}+h^2\hspace{1em}M{}_k(x,h)$

We can get $M{}_1$ from $M{}_2$ from this recurrence relation with $k=2$.

$M{}_1(x,h)=-\frac{B{}_2\hspace{1em}(f\hspace{1em}{}^{(1)}(x+h)-f\hspace{1em}{}^{(1)}(x))}{2!}+h^2\hspace{1em}M{}_2(x,h)$

We can get $M{}_1$ from $M{}_3$ from the recurrence relation.

$M{}_2(x,h)=-\frac{B{}_4\hspace{1em}(f\hspace{1em}{}^{(3)}(x+h)-f\hspace{1em}{}^{(3)}(x))}{4!}+h^2\hspace{1em}M{}_3(x,h)$

$M{}_1(x,h)=-\frac{B{}_2\hspace{1em}(f\hspace{1em}{}^{(1)}(x+h)-f\hspace{1em}{}^{(1)}(x))}{2!}-h^2\hspace{1em}\frac{B{}_4\hspace{1em}(f\hspace{1em}{}^{(3)}(x+h)-f\hspace{1em}{}^{(3)}(x))}{4!}+h^4\hspace{1em}M{}_3(x,h)$

We can get $M{}_1$ from $M{}_{l+1}$ from the recurrence relation.

$M{}_1(x,h)=-\sum _{k=1}^l\frac{h^{2\hspace{1em}k-2}\hspace{1em}B{}_{2\hspace{1em}k}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x+h)-f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x))}{(2\hspace{1em}k)!}+h^{2\hspace{1em}l}\hspace{1em}M{}_{l+1}(x,h)$

Substitute this for $M{}_1$ in a previous equation.

$\frac{f\hspace{1em}(x+h)+f\hspace{1em}(x)}{2}=\frac{{\int }_x^{x+h}f\hspace{1em}(y)\mathrm{dy}}{h}+$

$\sum _{k=1}^l\frac{h^{2\hspace{1em}k-1}\hspace{1em}B{}_{2\hspace{1em}k}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x+h)-f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x))}{(2\hspace{1em}k)!}-h^{2\hspace{1em}l+1}\hspace{1em}M{}_{l+1}(x,h)$

Now we can verify some examples using DERIVE.

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

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

Now using derive we can find

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})$

F(x):=

This is results on the previous page.

$\frac{f\hspace{1em}(x+h)+f\hspace{1em}(x)}{2}=\frac{{\int }_x^{x+h}f\hspace{1em}(y)\mathrm{dy}}{h}+$

$\sum _{k=1}^l\frac{h^{2\hspace{1em}k-1}\hspace{1em}B{}_{2\hspace{1em}k}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x+h)-f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x))}{(2\hspace{1em}k)!}-h^{2\hspace{1em}l+1}\hspace{1em}M{}_{l+1}(x,h)$

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

Expand the sumation into four parts.

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})=\frac{f\hspace{1em}\left(x_0\right)}{2}+\sum _{j=0}^{n-1}\frac{f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})}{2}+\sum _{j=1}^n\frac{f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})}{2}+\frac{f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}n\hspace{1em}h)}{2}$

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

Combine two of the parts.

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})=\frac{f\hspace{1em}\left(x_0\right)+f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}n\hspace{1em}h)}{2}+\sum _{j=0}^{n-1}\frac{f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})+f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh}\hspace{1em}+\hspace{1em}h)}{2}$

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

Substitute for the last part.

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})=\frac{f\hspace{1em}\left(x_0\right)+f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}n\hspace{1em}h)}{2}+\sum _{j=0}^{n-1}\frac{{\int }_{x_0+\mathrm{jh}}^{x_0+\mathrm{jh}+h}f\hspace{1em}(y)\mathrm{dy}}{h}+$

$\sum _{j=0}^{n-1}\sum _{k=1}^l\frac{h^{2\hspace{1em}k-1}\hspace{1em}B{}_{2\hspace{1em}k}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x_0+\mathrm{jh}+h)-f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x_0+\mathrm{jh}))}{(2\hspace{1em}k)!}-h^{2\hspace{1em}l+1}\hspace{1em}M{}_{l+1}(x_0+\mathrm{jh},h)$

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

Simplify the sums.

$\sum _{j=0}^{n}f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}\mathrm{jh})=\frac{f\hspace{1em}\left(x_0\right)+f\hspace{1em}(x_0\hspace{1em}+\hspace{1em}n\hspace{1em}h)}{2}+\frac{{\int }_{x_0}^{x_0+n\hspace{1em}h}f\hspace{1em}(y)\mathrm{dy}}{h}+$

$\sum _{k=1}^l\frac{h^{2\hspace{1em}k-1}\hspace{1em}B{}_{2\hspace{1em}k}\hspace{1em}(f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x_0+n\hspace{1em}h)-f\hspace{1em}{}^{(2\hspace{1em}k-1)}(x_0))}{(2\hspace{1em}k)!}-h^{2\hspace{1em}l+1}\hspace{1em}\sum _{j=0}^{n-1}M{}_{l+1}(x_0+\mathrm{jh},h)$