Creative Telescoping: Turning Multi-Dimensional Integrals into ODEs
2026-04-17
A triple integral for $\zeta(3)$. A double integral for $\zeta(2)$. A single integral for $\log 2$. Each is a “period” in the sense of Kontsevich–Zagier — a number defined by integrating a rational function over a rational domain. Creative telescoping is the technique that compresses these multi-dimensional integrals into ordinary differential equations in one variable. Here is how it works, and why it matters for computing these numbers.
The Problem: From Volumes to Curves
The Riemann zeta values have beautiful integral representations:
$$\zeta(2) = \int_0^1\!\int_0^1 \frac{dx\,dy}{1-xy}, \qquad \zeta(3) = \int_0^1\!\int_0^1\!\int_0^1 \frac{dx\,dy\,dz}{1-xyz}.$$These are multi-dimensional integrals of rational functions over the unit cube — “periods” in the language of Kontsevich and Zagier. But a polynomial ODE has only one time variable $t$. How do you turn a triple integral into a single-variable differential equation?
The answer is creative telescoping, a technique introduced by Doron Zeilberger in 1991 for sums, and extended to integrals by Frédéric Chyzak in 2000. The idea is both simple and powerful: differentiate the integral with respect to a parameter, and discover that the result satisfies an ODE.
Warm-Up: Evaluating $\zeta(3)$ by Iterated Integration
Before introducing the general machinery, let us do the $\zeta(3)$ integral by hand. The parametric family is:
$$I(t) = \int_0^1\!\int_0^1\!\int_0^1 \frac{dx\,dy\,dz}{1-t\,xyz}, \qquad I(1) = \zeta(3).$$Step 1: Integrate over $z$. Using $\int_0^1 \frac{dz}{1-w z} = -\frac{\log(1-w)}{w}$ for $|w| < 1$:
$$\int_0^1 \frac{dz}{1-t\,xyz} = \frac{-\log(1-txy)}{txy}.$$Step 2: Integrate over $y$. Expand $-\log(1-w)/w = \sum_{n=1}^{\infty} w^{n-1}/n$ and integrate term by term:
$$\int_0^1 \frac{-\log(1-txy)}{txy}\,dy = \sum_{n=1}^{\infty} \frac{(tx)^{n-1}}{n^2} = \frac{\mathrm{Li}_2(tx)}{tx}.$$Step 3: Integrate over $x$. Similarly:
$$\int_0^1 \frac{\mathrm{Li}_2(tx)}{tx}\,dx = \sum_{n=1}^{\infty} \frac{t^{n-1}}{n^3} = \frac{\mathrm{Li}_3(t)}{t}.$$So $I(t) = \mathrm{Li}_3(t)/t$, and indeed $I(1) = \mathrm{Li}_3(1) = \zeta(3)$.
We verified this numerically by computing the triple integral at several values of $t$ and comparing with $\mathrm{Li}_3(t)/t$ — they match to machine precision ($\sim 10^{-16}$).
The Theta Operator and the ODE
Now that we know $I(t) = \mathrm{Li}_3(t)/t$, we need the ODE that $\mathrm{Li}_3$ satisfies. The key tool is the theta operator $\theta = t\frac{d}{dt}$, which acts cleanly on power series:
$$\theta\!\left(\sum a_n t^n\right) = \sum n\,a_n\,t^n.$$Applied to polylogarithms, it peels off one layer of the $1/n^k$ factor:
$$\theta\,\mathrm{Li}_3(t) = \sum \frac{t^n}{n^2} = \mathrm{Li}_2(t), \qquad \theta\,\mathrm{Li}_2(t) = \sum \frac{t^n}{n} = -\log(1-t).$$One more application gives:
$$\theta\bigl(-\log(1-t)\bigr) = t \cdot \frac{1}{1-t} = \frac{t}{1-t}.$$Chaining these together:
$$\theta^3\,\mathrm{Li}_3(t) = \frac{t}{1-t}, \qquad \text{hence} \qquad (1-t)\,\theta^3\,\mathrm{Li}_3(t) = t.$$Expanding $\theta^3 = t^3 D^3 + 3t^2 D^2 + tD$ (where $D = d/dt$), this becomes:
$$\boxed{t^2(1-t)\,f''' + 3t(1-t)\,f'' + (1-t)\,f' = 1}$$— a third-order linear ODE with polynomial coefficients and a rational forcing term. This is the Picard–Fuchs equation for the period $\zeta(3)$ (via the polylogarithmic route).
What Is Creative Telescoping?
The hand computation above — integrate one variable at a time, recognize polylogarithms, apply $\theta$ — works because we know the special functions involved. But what if the integrand is less familiar?
Creative telescoping automates this process. The core idea, for a single integral:
$$L \cdot F(x, t) = \frac{\partial}{\partial x}\bigl[G(x, t)\bigr].$$
Integrating both sides over $x \in [0, 1]$:
$$L \cdot \int_0^1 F(x, t)\,dx = \bigl[G(x, t)\bigr]_0^1.$$If the boundary terms are known (often zero or simple), this gives an ODE for $\int_0^1 F\,dx$. The operator $L$ is the telescoper, and $G$ is the certificate.
Why “Telescoping”?
The name comes from the discrete version. For a sum $S_n = \sum_k F(n, k)$, a telescoping identity takes the form:
$$L \cdot F(n, k) = G(n, k+1) - G(n, k).$$Summing over $k$ collapses the right side (like a telescope), leaving $L \cdot S_n$ equal to boundary terms. This is exactly Zeilberger’s method for proving hypergeometric identities and finding recurrences.
The continuous version replaces the discrete difference $\Delta_k$ with $\partial/\partial x$, and sums with integrals.
The Algorithm (Sketch)
For holonomic (D-finite) integrands — those satisfying polynomial-coefficient differential equations in each variable — the creative telescoping algorithm proceeds roughly as:
- Ansatz: Assume $L$ has order $r$ with unknown polynomial coefficients of bounded degree.
- Closure properties: Compute the left ideal of differential operators annihilating $F$ in the Weyl algebra $\mathbb{Q}\langle t, \partial_t, x, \partial_x \rangle$.
- Elimination: Find an operator in this ideal that involves only $t$ and $\partial_t$ (plus an $\partial_x$-exact remainder). This is a Gröbner basis computation in the Weyl algebra.
- Output: The $t$-only part is the telescoper $L$; the $\partial_x$-exact part gives the certificate $G$.
Chyzak (2000) proved that this always terminates for holonomic integrands, guaranteeing that a finite-order ODE exists. Koutschan’s HolonomicFunctions Mathematica package implements this algorithm.
Multi-Dimensional Integrals
For a double integral $\int_0^1\!\int_0^1 F(x, y, t)\,dx\,dy$, apply creative telescoping twice: first eliminate $y$ to get an ODE in $(x, t)$, then eliminate $x$ to get an ODE in $t$ alone.
For $\zeta(3) = \int_0^1\!\int_0^1\!\int_0^1 \frac{dx\,dy\,dz}{1-txyz}$, the algorithm would:
- Telescope over $z$ → get an ODE in $(x, y, t)$
- Telescope over $y$ → get an ODE in $(x, t)$
- Telescope over $x$ → get the ODE in $t$ alone
The result is the same ODE we derived by hand: $t^2(1-t)f''' + 3t(1-t)f'' + (1-t)f' = 1$.
The HolonomicFunctions Package
Christoph Koutschan at RISC (Johannes Kepler University Linz) maintains the HolonomicFunctions Mathematica package, which implements creative telescoping for both sums and integrals.
Installation
Download from RISC’s software page and place in Mathematica’s $UserBaseDirectory/Applications/ directory.
Basic Usage
|
|
The output is a differential operator in $\partial_t$ plus boundary terms. For simple rational integrands, the computation is fast (seconds). For more complex integrands (e.g., involving algebraic functions or high-dimensional integrals), it can take longer.
Classical Examples
Here are some well-known periods and what creative telescoping gives:
| Period | Integral | ODE order |
|---|---|---|
| $\log 2$ | $\int_0^1 \frac{dx}{1+x}$ | 1 |
| $\pi/4$ | $\int_0^1 \frac{dx}{1+x^2}$ | 1 |
| $\zeta(2)$ | $\int_0^1\!\int_0^1 \frac{dx\,dy}{1-xy}$ | 2 |
| $\zeta(3)$ | $\int_0^1\!\int_0^1\!\int_0^1 \frac{dx\,dy\,dz}{1-xyz}$ | 3 |
| $\mathrm{K}(1/\sqrt{2})$ | $\int_0^{\pi/2} \frac{d\theta}{\sqrt{1-\frac{1}{2}\sin^2\theta}}$ | 2 |
The ODE order equals the number of integration variables for the $\zeta(n)$ integrals — each layer of integration adds one order to the differential equation.
Why This Matters: The Pipeline from Periods to Polynomial ODEs
In a previous post, we explored computing $\zeta(3)$ via a polynomial ODE (PIVP). The challenge was: how do you go from the definition of a number (as an integral) to a computation of that number (as the limit of an ODE solution)?
Creative telescoping is the first step of a general pipeline:
$$\boxed{\text{Period (multi-dim integral)} \xrightarrow{\text{creative telescoping}} \text{ODE} \xrightarrow{\text{time reparam.}} \text{Polynomial PIVP}}$$- Creative telescoping compresses the multi-dimensional integral into a single-variable ODE with polynomial coefficients.
- Time reparametrization $d\tau = dz/p(z)$ absorbs the ODE’s leading coefficient $p(z)$ into the time variable, making the system polynomial.
- Normalization (ratio variables, if needed) keeps the solution bounded.
For first-order periods like $\log 2$ and $\pi$, the ODE is trivial and the pipeline gives a 3-variable polynomial PIVP with exact rational initial conditions. For higher-order periods like $\zeta(2)$ and $\zeta(3)$, the situation is more interesting: the ODE has a maximally unipotent monodromy (MUM) point at the origin, which freezes the PIVP. The workaround — starting from a nearby rational point with truncated Taylor series as initial conditions — works in practice but raises subtle theoretical questions.
We proved (see the working note) that every first-order period is in $R_{\mathrm{RTCRN}}$ (the class of CRN-computable numbers), and conjectured that all periods are CRN-computable. Creative telescoping is the mathematical bridge that makes this conjecture precise.
Further Reading
- D. Zeilberger, The method of creative telescoping, J. Symbolic Comput. 11 (1991), 195–204. The original paper.
- F. Chyzak, An extension of Zeilberger’s fast algorithm to general holonomic functions, Discrete Math. 217 (2000), 115–134. The generalization to continuous integrals.
- C. Koutschan, HolonomicFunctions package. The practical implementation.
- M. Kontsevich and D. Zagier, Periods, in Mathematics Unlimited — 2001 and Beyond, Springer, 2001, pp. 771–808. The definition of periods that motivates this direction.
- A. Bostan, F. Chyzak, M. van Hoeij, M. Kauers, and L. Pech, Hypergeometric expressions for generating functions of walks with small steps in the quarter plane, European J. Combin. 61 (2017), 242–275. A beautiful application of creative telescoping to combinatorics.
Creative telescoping doesn’t just evaluate integrals — it reveals the differential equation hiding inside them. And once you have the ODE, the entire machinery of polynomial dynamics becomes available. The distance from a triple integral to a polynomial ODE is shorter than it looks.
积分之中有微分,微分之中有机器。
Within the integral hides a differential equation; within the equation hides a machine.