Central difference interpolation

Newton's forward and backward interpolation formulae are lopsided: each leans on differences taken from only one side of a chosen starting point, so they converge fast only when that point sits near an end of the table. What happens when the value you actually need lies somewhere in the middle? Pushing a one-sided formula through a mid-table point wastes accuracy — the series converges slowly because $p=(x-x_0)/h$ is forced to be large. Central-difference formulae fix this by centering the difference table on the tabulated point nearest $x$ and drawing symmetrically from both sides. Gauss's forward and backward formulae are the building blocks; averaging them gives Stirling's formula, the workhorse for mid-table interpolation. This article builds the idea, derives Stirling's formula from the Gauss pair, and reworks the syllabus example in full.

The core insight

For a point near the middle of an equally spaced table, averaging Gauss's forward and backward formulae into Stirling's formula uses differences symmetric about the nearest tabulated point $x_0$ , keeping $p=(x-x_0)/h$ small and the series fast-converging.

The core idea

Why one-sided formulae struggle mid-table. Newton's forward formula fixes $x_0$ at the first row and marches forward through $\Delta f_0,\Delta^2f_0,\ldots$ ; Newton's backward formula fixes $x_0$ at the last row and marches backward. Either way, if the interpolation point $x$ sits in the interior, $p=(x-x_0)/h$ is far from zero for whichever endpoint you anchor at, so the truncated polynomial needs many terms to be accurate. The fix is to anchor $x_0$ at the tabulated point nearest $x$ , so $p$ stays small (typically $|p|\le 0.5$ ), regardless of where in the table $x$ falls.

Gauss's forward formula (for $0<p<0.5$ , i.e. $x$ just ahead of $x_0$ ) uses one diagonal of differences that zig-zags outward from $x_0$ :

f(x) = f_0 + p\,\Delta f_0 + \frac{p(p-1)}{2!}\Delta^2 f_{-1} + \frac{(p+1)p(p-1)}{3!}\Delta^3 f_{-1} + \frac{(p+1)p(p-1)(p-2)}{4!}\Delta^4 f_{-2} + \cdots

where: $f_0=f(x_0)$ is the table value at the origin row, $\Delta^k f_j$ denotes the $k$ -th forward difference starting at the row $j$ steps from $x_0$ (negative $j$ means $j$ rows above $x_0$ ), and $p=(x-x_0)/h$ .

Gauss's backward formula (for $-0.5<p<0$ , i.e. $x$ just behind $x_0$ ) uses the mirror-image diagonal:

f(x) = f_0 + p\,\Delta f_{-1} + \frac{(p+1)p}{2!}\Delta^2 f_{-1} + \frac{(p+1)p(p-1)}{3!}\Delta^3 f_{-2} + \cdots

Averaging the two gives Stirling's formula, which is symmetric about $x_0$ and is the standard choice whenever $-0.25\le p\le 0.25$ :

f(x) = f_0 + p\,\frac{\Delta f_{-1}+\Delta f_0}{2} + \frac{p^2}{2!}\Delta^2 f_{-1} + \frac{p(p^2-1)}{3!}\cdot\frac{\Delta^3 f_{-2}+\Delta^3 f_{-1}}{2} + \frac{p^2(p^2-1)}{4!}\Delta^4 f_{-2} + \cdots

The differences needed — $\Delta f_{-1}$ , $\Delta f_0$ , $\Delta^2 f_{-1}$ , $\Delta^3 f_{-2}$ , $\Delta^3 f_{-1}$ , $\Delta^4 f_{-2}$ — sit around $x_0$ in a diamond shape when written into a difference table, which is why picking the right origin row matters as much as picking the right formula. Like Newton's formulae, all three of these are still truncated polynomial fits through the tabulated points; nothing about the underlying mathematics changes. The only thing that changes is bookkeeping — which diagonal of the difference table is used, and where it is centred — chosen so that $p$ stays small and the series converges in as few terms as possible for a point that happens to sit mid-table rather than near an end.

Figure

A difference table diamond centered on the origin row x zero showing which forward differences feed Stirling formula

Derivation Aof Stirling's formula from Gauss's forward and backward formulae

The goal is to show that Stirling's formula is nothing new — it is simply the term-by-term average of Gauss's forward and Gauss's backward formulae, and to confirm that averaging is what makes each odd-order difference term appear as a symmetric average of two neighbouring differences rather than a single lopsided one.

Step 1. Write both Gauss formulae with matching term order. Gauss forward: $f_0+p\,\Delta f_0+\dfrac{p(p-1)}{2!}\Delta^2f_{-1}+\dfrac{(p+1)p(p-1)}{3!}\Delta^3f_{-1}+\cdots$ . Gauss backward: $f_0+p\,\Delta f_{-1}+\dfrac{(p+1)p}{2!}\Delta^2f_{-1}+\dfrac{(p+1)p(p-1)}{3!}\Delta^3f_{-2}+\cdots$ . Both share the same constant term $f_0$ and the same second-order coefficient's difference $\Delta^2f_{-1}$ .

Step 2. Average the two first-order terms. Forward contributes $p\,\Delta f_0$ , backward contributes $p\,\Delta f_{-1}$ ; their average is $p\cdot\dfrac{\Delta f_{-1}+\Delta f_0}{2}$ — exactly the first-order term seen in Stirling's formula above.

Step 3. Average the second-order terms. Forward's coefficient is $\dfrac{p(p-1)}{2!}$ and backward's is $\dfrac{(p+1)p}{2!}$ ; both multiply the same difference $\Delta^2f_{-1}$ , so averaging the coefficients gives $\dfrac{1}{2}\left[\dfrac{p(p-1)+(p+1)p}{2!}\right]=\dfrac{1}{2}\cdot\dfrac{2p^2}{2!}=\dfrac{p^2}{2!}$ , matching Stirling's second-order term.

Step 4. Average the third-order terms. Both share the coefficient $\dfrac{(p+1)p(p-1)}{3!}$ , but forward attaches it to $\Delta^3f_{-1}$ while backward attaches it to $\Delta^3f_{-2}$ — different differences. Averaging gives $\dfrac{(p+1)p(p-1)}{3!}\cdot\dfrac{\Delta^3f_{-2}+\Delta^3f_{-1}}{2}$ , and writing $(p+1)p(p-1)=p(p^2-1)$ reproduces the third-order term of Stirling's formula exactly. Higher-order terms follow the same pattern: even orders share one difference and their coefficients simply average; odd orders keep the same coefficient but the two differences themselves get averaged.

Worked example

Given the table $x=20,25,30,35,40$ with $f(x)=512,552,585,610,628$ (spacing $h=5$ ), Stirling's formula centered at $x_0=30$ was used earlier to estimate $f(29)$ . Here the same difference table is reused to estimate $f(31)$ , confirming the near-symmetric behaviour predicted by the derivation above.

Step 1. The difference table centered at $x_0=30$ (built once, reused for any $x$ near 30) gives $\Delta f_{-2}=40$ , $\Delta f_{-1}=33$ , $\Delta f_0=25$ , $\Delta f_1=18$ ; $\Delta^2f_{-2}=-7$ , $\Delta^2f_{-1}=-8$ , $\Delta^2f_0=-7$ ; $\Delta^3f_{-2}=-1$ , $\Delta^3f_{-1}=1$ ; $\Delta^4f_{-2}=2$ .

Step 2. Compute $p=(31-30)/5=0.2$ — small, as required for Stirling's formula to be accurate.

Step 3. Substitute term by term into Stirling's formula with $f_0=585$ : Term 1 is $585$ . Term 2 is $(0.2)\cdot\dfrac{33+25}{2}=(0.2)(29)=5.8$ . Term 3 is $\dfrac{0.04}{2}(-8)=-0.16$ . Term 4 is $\dfrac{(0.2)(0.04-1)}{6}\cdot\dfrac{-1+1}{2}=0$ (the two third differences cancel here, so this term vanishes regardless of the coefficient). Term 5 is $\dfrac{(0.04)(0.04-1)}{24}(2)\approx-0.0032$ .

Step 4. Sum the terms: $f(31)\approx 585+5.8-0.16+0-0.0032\approx \mathbf{590.64}$ .

Competitive stretch

Compare with $f(29)\approx579.04$ from the same table (obtained with $p=-0.2$ , same origin). Because $p=+0.2$ and $p=-0.2$ are mirror displacements, the odd-order terms (Term 2, Term 4) flip sign while the even-order terms (Term 3, Term 5) stay identical — a direct consequence of Step 2 and Step 4 of the derivation above, and a fast sanity check: $585+(-5.8-0.16+0-0.0032)=579.04$ and $585+(5.8-0.16+0-0.0032)=590.64$ differ only in the sign of the odd terms, exactly as the symmetric structure predicts.

Common traps

Choosing $x_0$ far from $x$ instead of the nearest tabulated point — this defeats the whole purpose of a central formula by making $|p|$ large again.
Mixing up which diagonal a difference like $\Delta^3f_{-2}$ belongs to — central formulae draw from both sides of $x_0$ , unlike Newton's one-sided formulae.
Forgetting to average both the coefficients (even orders) and the differences themselves (odd orders) when deriving or applying Stirling's formula — skipping either average breaks the symmetry that makes it accurate.
Applying Stirling's formula too close to a table edge, where there are not enough points on one side to build the required differences — fall back to Newton's forward or backward formula instead.

Ask the tutor: If building the difference diamond for a new table feels confusing, ask the AI tutor to walk through constructing it row by row for a fresh set of data points.

Why does Stirling formula need p between minus 0.25 and 0.25 while Gauss formulae only need p between minus 0.5 and 0.5?

Stirling's formula is built by averaging Gauss forward and Gauss backward, and that averaging cancels some error but does not remove the underlying truncation error of a polynomial fit. In practice the averaged series converges fastest, and stays most accurate, closest to the origin $x_0$ , so the tighter window $-0.25\le p\le0.25$ is where Stirling's formula outperforms using either Gauss formula alone; outside that window Gauss forward or backward (whichever fits) is still fine up to $|p|=0.5$ .

If Term 4 in the worked example is zero does that mean the third-difference correction never matters?

No — it was zero here only because $\Delta^3f_{-2}=-1$ and $\Delta^3f_{-1}=1$ happened to average to zero for this particular table. In general the two third differences will not cancel, and Term 4 contributes a real correction. Always compute both third differences and average them rather than assuming the term vanishes; this table's cancellation is a coincidence of the data, not a property of Stirling's formula.

How do I decide whether to center the difference table at x0 equal to 30 or some other row?

Pick the tabulated $x_0$ nearest to the interpolation point $x$ , so that $p=(x-x_0)/h$ has the smallest possible magnitude. For $x=29$ or $x=31$ in this table, $x_0=30$ gives $|p|=0.2$ , well inside Stirling's accurate range. If you instead anchored at $x_0=25$ , you would get $p=0.8$ or $p=1.2$ , far too large for a central formula to converge quickly — you would need many more terms for the same accuracy.

Central difference interpolation

The core idea

Worked example

Competitive stretch

Common traps

Before this

Where this leads