# 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));
answer := 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
> readlib(eulermac):
> 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.
# What about using PartialSum(10)?
#
> 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.
