# 1.28 Lesson 28

```# Lesson 28.  Euler-Maclaurin Summation Formula
# ---------------------------------------------
#
> restart:
# Suppose f is a smooth function, and F an anti-difference of f, i.e.
> f(x) = F(x+1)-F(x):
# We expanded F(x+1) (formally, i.e. we don't worry about convergence) using a Taylor
# series involving derivatives of F at x.
> F(x+1) = F(x) + Sum((D@@k)(F)(x)/k!, k=1..infinity);

/infinity           \
| -----    (k)      |
|  \      D   (F)(x)|
F(x + 1) = F(x) + |   )     ----------|
|  /          k!    |
| -----             |
\ k = 1             /
# So our formula for f becomes
> subs(", "");

infinity
-----    (k)
\      D   (F)(x)
f(x) =    )     ----------
/          k!
-----
k = 1
# Basically what we want to do is "invert" this formula to write F(x) in terms of
# f(x), its antiderivative and its derivatives.  Let's work with a finite partial sum.
> f(x) = Sum((D@@k)(F)(x)/k!, k=1 .. 7) + error;

/  7             \
|-----  (k)      |
| \    D   (F)(x)|
f(x) = |  )   ----------| + error
| /        k!    |
|-----           |
\k = 1           /
# where "error" will involve higher order derivatives of F.  In fact, according to
# Taylor's Theorem, "error" in this case is (D@@8)(F)(s)/8! for some s between x
# and x+1.  We suppose (D@@8)(F)(s) -> 0 (hopefully, fairly rapidly) as s -> infinity,
# so we get some control on the error.
# For convenience, make the following definition:
> Q:= m -> sum((D@@(k+m))(F)(x)/k!, k = 1 .. 7-m);

7 - m
-----  (k + m)
\    D       (F)(x)
Q := m ->   )   --------------
/          k!
-----
k = 1
> eqn[0]:= f(x) = Q(0);

(2)              (3)               (4)
eqn[0] := f(x) = D(F)(x) + 1/2 D   (F)(x) + 1/6 D   (F)(x) + 1/24 D   (F)(x)

(5)                (6)                 (7)
+ 1/120 D   (F)(x) + 1/720 D   (F)(x) + 1/5040 D   (F)(x)
# We're neglecting the error here, but remember that it involves (D@@8)(F).
> eqn[1]:= D(f)(x) = Q(1);

(2)              (3)              (4)
eqn[1] := D(f)(x) = D   (F)(x) + 1/2 D   (F)(x) + 1/6 D   (F)(x)

(5)                (6)                (7)
+ 1/24 D   (F)(x) + 1/120 D   (F)(x) + 1/720 D   (F)(x)
# This just came from differentiating eqn[0].  We can do this as many times as we want.
# But since we're only including terms up to (D@@7)(F), it is only useful up to eqn[6].
> for kk from 2 to 6 do eqn[kk]:= (D@@kk)(f)(x) = Q(kk) od:
> eqn[6];

(6)          (7)
D   (f)(x) = D   (F)(x)
# None of these equations yet involved F(x) itself.  But we can integrate eqn[0].
> eqn[-1]:= int(f(x), x) = Q(-1);

/
|                                      (2)               (3)
eqn[-1] :=  |  f(x) dx = F(x) + 1/2 D(F)(x) + 1/6 D   (F)(x) + 1/24 D   (F)(x)
|
/

(4)                (5)                 (6)
+ 1/120 D   (F)(x) + 1/720 D   (F)(x) + 1/5040 D   (F)(x)

(7)
+ 1/40320 D   (F)(x)
# Note that we aren't specifying which antiderivative of f to take, i.e. there is an unknown
# constant here.
# Now the idea is to solve for F(x).
> solve({ seq(eqn[kk], kk=-1 .. 6) }, { seq((D@@kk)(F)(x), kk = 0 .. 7) }):
# I didn't show the whole result, because I only want the F(x).
> answer:= F(x) = subs(", F(x));

/
|                     (5)                                          (3)
|  f(x) dx + 1/30240 D   (f)(x) - 1/2 f(x) + 1/12 D(f)(x) - 1/720 D   (f)(x)
|
/
# Again, this is not exact, but the error will involve (D@@8)(F), i.e. if
# |(D@@8)(F)(s)| < K(x) for x <= s <= x+1, then |error| < const * K(x).
#
# The same result could be obtained as follows:
# formally  f(x) = (exp(D)-1) F(x) (if we think of "D" as an operator that can be
# manipulated algebraically, with multiplication corresponding to composition).
# "exp(D)-1" just means: expand exp(u)-1 in a series in powers of u, and replace
#  u^n by (D@@n).
# So F(x) = (1/(exp(D)-1)) f(x) = (D/(exp(D)-1)) D^(-1) f(x)
# where "D/(exp(D)-1)" means  expand g(u) = u/(exp(u)-1) in a series in powers of
# u, and replace u^n by (D@@n).
# D^(-1) f(x) means an antiderivative of f: we don't know which one, so we include an
# arbitrary constant.
> taylor(u/(exp(u)-1), u, 9);

2          4            6      8
1 - 1/2 u + 1/12 u  - 1/720 u  + 1/30240 u  + O(u )
> int(f(x), x) - 1/2 * f(x) + 1/12 * D(f)(x) - 1/720 * (D@@3)(f)(x) + 1/30240 * (D@@5)(f)(x);

/
|                     (5)                                          (3)
|  f(x) dx + 1/30240 D   (f)(x) - 1/2 f(x) + 1/12 D(f)(x) - 1/720 D   (f)(x)
|
/
# We just kept the derivatives of F up to (D@@7)(F), but of course this could be done
# with as many terms as you want.  The result is the Euler-Maclaurin summation formula.
# In general,
> F(x) = int(f(x), x) + Sum(c[k+1]*(D@@k)(f)(x), k = 0 .. N) + error;

/  N                      \
/           |-----                    |
|            | \              (k)      |
F(x) =  |  f(x) dx + |  )   c[k + 1] D   (f)(x)| + error
|            | /                       |
/             |-----                    |
\k = 0                    /
# where "error" involve (D@@(N+2))(F), and
> u/(exp(u)-1) = 1 + Sum(c[k]*u^k, k=1..infinity);

/infinity        \
| -----          |
u           |  \            k|
---------- = 1 + |   )     c[k] u |
exp(u) - 1       |  /             |
| -----          |
\ k = 1          /
> taylor(u/(exp(u)-1), u, 15);

2          4            6              8               10
1 - 1/2 u + 1/12 u  - 1/720 u  + 1/30240 u  - 1/1209600 u  + 1/47900160 u

691       12      14
- ------------- u   + O(u  )
1307674368000
# It looks like (except for the term in u), there are only even powers of u here, and so
# (except for the f(x) term) only odd derivatives of f will appear in the Euler-Maclaurin
# formula.  In fact, u/(exp(u)-1) + u/2 is an even function of u.
> u/(exp(u)-1) + u/2 - subs(u=-u, u/(exp(u)-1) + u/2);

u                 u
---------- + u + ------------
exp(u) - 1       exp(- u) - 1
> normal(");

u (- 1 + exp(u) exp(- u))
---------------------------
(exp(u) - 1) (exp(- u) - 1)
> simplify(");

0
> eulermac(f(x),x,4);

/                                               /   3      \        5
|                            /  d      \         |  d       |       d
|  f(x) dx - 1/2 f(x) + 1/12 |---- f(x)| - 1/720 |----- f(x)| + O(----- f(x))
|                            \ dx      /         |   3      |        5
/                                                 \ dx       /      dx
> eulermac(f(x),x,5);

/                                               /   3      \
|                            /  d      \         |  d       |
|  f(x) dx - 1/2 f(x) + 1/12 |---- f(x)| - 1/720 |----- f(x)|
|                            \ dx      /         |   3      |
/                                                 \ dx       /

/   5      \        7
|  d       |       d
+ 1/30240 |----- f(x)| + O(----- f(x))
|   5      |        7
\ dx       /      dx
> eulermac(f(x),x,6);

/                                               /   3      \
|                            /  d      \         |  d       |
|  f(x) dx - 1/2 f(x) + 1/12 |---- f(x)| - 1/720 |----- f(x)|
|                            \ dx      /         |   3      |
/                                                 \ dx       /

/   5      \        7
|  d       |       d
+ 1/30240 |----- f(x)| + O(----- f(x))
|   5      |        7
\ dx       /      dx
# Example: in Lessons 21 and 22 we wanted to calculate the sum of the following series:
> f:= k -> ln(k+1)/k^3:
> S:= Sum(f(k), k=1..infinity);

infinity
-----
\      ln(k + 1)
S :=    )     ---------
/           3
-----       k
k = 1
> PartialSum:= N -> Sum(f(k), k=1..N);

N
-----
\
PartialSum := N ->   )   f(k)
/
-----
k = 1
> F:= N ->  PartialSum(N-1) - S;

F := N -> PartialSum(N - 1) - S
> F(N+1) - F(N) = f(N);

/  N            \   /N - 1          \
|-----          |   |-----          |
| \    ln(k + 1)|   | \    ln(k + 1)|   ln(N + 1)
|  )   ---------| - |  )   ---------| = ---------
| /         3   |   | /         3   |        3
|-----     k    |   |-----     k    |       N
\k = 1          /   \k = 1          /
# So the sum of the series will be
> PartialSum(N-1) - F(N);

infinity
-----
\      ln(k + 1)
)     ---------
/           3
-----       k
k = 1
> em:= eulermac(f(k), k,5);

1                    ln(k + 1) (k + 1) (k - 1)       ln(k + 1)
em := - --- - 1/2 ln(k) + 1/2 ------------------------- - 1/2 ---------
2 k                              2                         3
(k + 1)  - 1 - 2 k              k

1             ln(k + 1)          1                 1
+ ------------- - 1/4 --------- - --------------- - --------------
3            4                 3  3             2  4
12 (k + 1) k            k       360 (k + 1)  k    80 (k + 1)  k

1              ln(k + 1)           1                 1
- ------------- + 1/12 --------- + ---------------- + ---------------
5             6                  5  3              4  4
20 (k + 1) k             k       1260 (k + 1)  k    336 (k + 1)  k

1                 5                5              ln(k + 1)
+ --------------- + --------------- + ------------- - 1/12 --------- + O(
3  5              2  6               7             8
126 (k + 1)  k    252 (k + 1)  k    84 (k + 1) k             k

720           2520          6048         12600         25200
----------- + ----------- + ----------- + ----------- + -----------
7  3          6  4          5  5          4  6          3  7
(k + 1)  k    (k + 1)  k    (k + 1)  k    (k + 1)  k    (k + 1)  k

52920        141120            ln(k + 1)
+ ----------- + ---------- - 181440 ---------)
2  8            9              10
(k + 1)  k    (k + 1) k              k
> oterm:= op(nops(em), em);

720           2520          6048         12600         25200
oterm := O(----------- + ----------- + ----------- + ----------- + -----------
7  3          6  4          5  5          4  6          3  7
(k + 1)  k    (k + 1)  k    (k + 1)  k    (k + 1)  k    (k + 1)  k

52920        141120            ln(k + 1)
+ ----------- + ---------- - 181440 ---------)
2  8            9              10
(k + 1)  k    (k + 1) k              k
# So the "O" term is really O(ln(k+1)/k^10) (or O(ln(k)/k^10)).
# How well can we do using PartialSum(3)) and the "eulermac" approximation for F(4)?
> Digits:= 30:\
p:=evalf( PartialSum(3));

p := .881817952240492006724395710996
> S:=  [seq(eulermac(f(N),N,2*dd-1), dd=1..20)]:
> S[1];

1                    ln(N + 1) (N + 1) (N - 1)       ln(N + 1)
- --- - 1/2 ln(N) + 1/2 ------------------------- - 1/2 ---------
2 N                              2                         3
(N + 1)  - 1 - 2 N              N

1             ln(N + 1)
+ ------------- - 1/4 ---------
3            4
12 (N + 1) N            N

2             9            36          ln(N + 1)
+ O(----------- + ----------- + ---------- - 60 ---------)
3  3          2  4            5           6
(N + 1)  N    (N + 1)  N    (N + 1) N           N
> S:= subs(N=4, O=0, S):
> S:= evalf(S);

S := [-.077608192842046927480099003669, -.077587514716115786657590578933,

-.077588556277433903593917990399, -.077588458201894897432628775656,

-.077588472828733987882426461297, -.077588469667696606810961643052,

-.077588470601626899014664065564, -.077588470239897919377653248088,

-.077588470417793155396485886173, -.077588470309472422408564338691,

-.077588470389485285716290381652, -.077588470318986538207610671826,

-.077588470392036307061297951691, -.077588470304088314740403192465,

-.077588470425841417342257207269, -.077588470233778893208110489890,

-.077588470576273214502727965277, -.077588469890708723502297478972,

-.077588471421428500160147553478, -.077588467630506371144331916901]
> seq(S[nn]-S[nn-1], nn=2..20);

-5
.000020678125931140822508424736, -.1041561318116936327411466*10  ,

-7                              -7
.98075539006161289214743*10  , -.14626839090449797685641*10  ,

-8                            -9
.3161037381071464818245*10  , -.933930292203702422512*10  ,

-9                            -9
.361728979637010817476*10  , -.177895236018832638085*10  ,

-9                           -10
.108320732987921547482*10  , -.80012863307726042961*10   ,

-10                           -10
.70498747508679709826*10   , -.73049768853687279865*10   ,

-10                            -9
.87947992320894759226*10   , -.121753102601854014804*10  ,

-9                            -9
.192062524134146717379*10  , -.342494321294617475387*10  ,

-9                             -8
.685564491000430486305*10  , -.1530719776657850074506*10  ,

-8
.3790922129015815636577*10
# The terms, which are indicative of the error, reach a minimum at about nn = 12, then get
# bigger again.
# So S[11] or S[12] probably gives us our best estimate, which should be within 10^(-10).
> a11:= p - S[11]; a12:= p - S[12];

a11 := .959406422629977292440686092648

a12 := .959406422559478544932006382822
# For comparison, our last approximation from Lesson 22, with an estimated error less than
# .5*10^(-10), was
> A155:= .95940642259480509054;

A155 := .95940642259480509054
> a11 - A155, a12 - A155;

-10                           -10
.35172201900686092648*10   , -.35326545607993617178*10
# Of course, 3 is an almost ridiculously small number of terms to take in a partial sum.
#
> p:= evalf(PartialSum(10));

p := .946218752907129242715574153661
> S:= [seq(eulermac(f(N),N,2*dd-1), dd=11 .. 12)]:
> S:= subs(N=11, O=0, S):
> S:= evalf(S);

S := [-.01318766968687974052279721719, -.0131876696868797405221710290318]
> S[2]-S[1];

-21
.6261881582*10
> [p-S[1], p-S[2]];

[.959406422594008983238371370851, .959406422594008983237745182693]
# The true answer is (probably) between these numbers.

```