infsup

by Zinan Huang 🌸

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:

  1. Ansatz: Assume $L$ has order $r$ with unknown polynomial coefficients of bounded degree.
  2. 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$.
  3. 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.
  4. 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:

  1. Telescope over $z$ → get an ODE in $(x, y, t)$
  2. Telescope over $y$ → get an ODE in $(x, t)$
  3. 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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
<< HolonomicFunctions`

(* Integrand for ζ(3): 1/(1-t·x·y·z) *)
F = 1/(1 - t*x*y*z);

(* Creative telescoping: eliminate x, y, z *)
ann = Annihilator[F, {Der[t], Der[x], Der[y], Der[z]}];

(* Find the ODE in t only *)
tel = CreativeTelescoping[F, Der[x], {0, 1}];

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}}$$
  1. Creative telescoping compresses the multi-dimensional integral into a single-variable ODE with polynomial coefficients.
  2. Time reparametrization $d\tau = dz/p(z)$ absorbs the ODE’s leading coefficient $p(z)$ into the time variable, making the system polynomial.
  3. 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


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.