Introduction

A derivative measures how fast something changes. It is the cornerstone of physics, engineering, and economics — but it is also an ideal object. Derivatives are defined in the limit as the gap between two points shrinks to zero, and no computer can store an infinitely thin gap.

The finite difference method breaks that idealism: instead of letting the gap shrink to zero, keep it small but finite. Pick a set of evenly spaced grid points and replace each derivative with a simple ratio:

f′(x)≈f(x+h)−f(x)hf'(x) \approx \frac{f(x+h) - f(x)}{h}

That one substitution converts a differential equation — full of "how fast is this changing right now?" — into a system of ordinary arithmetic operations. Write the same approximation at every grid point, and a differential equation turns into a matrix equation or a stepping rule the computer can crunch in a loop.

The method is not new: Brook Taylor laid the groundwork with his series in 1715, and Lewis Fry Richardson used finite differences to compute the first ever numerical weather forecast by hand in 1922 (it took him months; a modern laptop does it in milliseconds). Today finite differences are baked into weather models, structural analysis codes, electronic circuit simulators, and image-processing pipelines.

Try It: Heat Along a Bar

The 1-D heat equation asks: given an initial temperature profile along a thin bar, how does heat spread over time? The exact answer involves an infinite series of sines. The finite difference answer is a loop.

<div class="controls">
  <label>Grid points (N): <strong id="nVal">30</strong>
    <input type="range" id="nSlider" min="10" max="60" value="30" step="5">
  </label>
  <label>Stability r: <strong id="rVal">0.40</strong>
    <input type="range" id="rSlider" min="0.05" max="0.60" value="0.40" step="0.05">
  </label>
  <div class="btn-row">
    <button id="btnStep" type="button">Step</button>
    <button id="btnRun"  type="button">Run</button>
    <button id="btnStop" type="button" disabled>Stop</button>
    <button id="btnReset" type="button" class="ghost">Reset</button>
  </div>
  <div id="info" class="info">Step 0 — r = 0.40 (stable)</div>
</div>
<canvas id="chart" width="460" height="220"></canvas>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; }
.controls { display: flex; flex-direction: column; gap: .45rem; margin-bottom: .6rem; }
label { font-size: .88rem; display: flex; align-items: center; gap: .5rem; flex-wrap: wrap; }
label strong { min-width: 2.8rem; }
input[type=range] { flex: 1; min-width: 120px; cursor: pointer; }
.btn-row { display: flex; gap: .45rem; flex-wrap: wrap; }
button { font: 600 13px system-ui; padding: .38rem .8rem; border-radius: 7px;
         border: 1px solid #1d3557; background: #1d3557; color: #fff; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.info { font-size: .82rem; color: #555; min-height: 1.2em; }
canvas { border: 1px solid #dde3ea; border-radius: 8px; width: 100%; max-width: 460px;
         display: block; background: #f9fbfc; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Explicit finite difference scheme for the 1-D heat equation:
//   u_t = alpha * u_xx
// Discretised: u[i]^{n+1} = u[i]^n + r*(u[i+1]^n - 2*u[i]^n + u[i-1]^n)
// where r = alpha*dt/dx^2. Stable when r <= 0.5.
// Boundary conditions: u[0] = u[N-1] = 0 (bar ends held at 0 °C).
// Initial condition: triangular spike — rises to 1.0 at the midpoint.
// Exact solution: Fourier sine series (first 20 terms shown in blue).

const canvas = document.getElementById('chart');
const ctx    = canvas.getContext('2d');
const CW = canvas.width, CH = canvas.height;
const PAD = { top: 18, right: 14, bottom: 30, left: 38 };
const PW = CW - PAD.left - PAD.right;
const PH = CH - PAD.top  - PAD.bottom;

// --- Parameters (updated from sliders) ---
let N  = 30;   // number of grid points (including endpoints)
let r  = 0.40; // stability parameter alpha*dt/dx^2

// --- State ---
let u    = [];   // current temperature array
let step = 0;
let timer = null;

// --- Initial condition: triangular spike ---
function init() {
  u = new Array(N).fill(0);
  const mid = Math.floor((N - 1) / 2);
  for (let i = 1; i < N - 1; i++) {
    u[i] = i <= mid ? i / mid : (N - 1 - i) / (N - 1 - mid);
  }
  step = 0;
}

// --- One explicit FD step ---
function fdStep() {
  const next = u.slice();
  for (let i = 1; i < N - 1; i++) {
    next[i] = u[i] + r * (u[i + 1] - 2 * u[i] + u[i - 1]);
  }
  next[0] = 0; next[N - 1] = 0;
  u = next;
  step++;
}

// --- Exact solution: Fourier sine series (first 20 terms) ---
// For the triangular initial condition on [0,1], the exact Fourier coefficients
// of the odd extension are: b_n = (8/(n^2 * pi^2)) * sin(n*pi/2) for n odd, 0 for n even.
// (This is the standard Fourier sine series of a triangle on [0,1] with peak at 0.5.)
function exactU(t) {
  const ex = new Array(N).fill(0);
  const TERMS = 20;
  for (let i = 1; i < N - 1; i++) {
    const x = i / (N - 1); // x in (0,1)
    let val = 0;
    for (let n = 1; n <= TERMS; n++) {
      // b_n for triangle on [0,1] with peak at 0.5:
      const bn = (8 / (n * n * Math.PI * Math.PI)) * Math.sin(n * Math.PI / 2);
      val += bn * Math.exp(-n * n * Math.PI * Math.PI * t) * Math.sin(n * Math.PI * x);
    }
    ex[i] = Math.max(0, val);
  }
  return ex;
}

// --- Drawing ---
function draw() {
  ctx.clearRect(0, 0, CW, CH);

  // Background
  ctx.fillStyle = '#f9fbfc';
  ctx.fillRect(0, 0, CW, CH);

  // Axes
  ctx.strokeStyle = '#b0bec8';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(PAD.left, PAD.top);
  ctx.lineTo(PAD.left, PAD.top + PH);
  ctx.lineTo(PAD.left + PW, PAD.top + PH);
  ctx.stroke();

  // Y axis labels
  ctx.fillStyle = '#888';
  ctx.font = '11px system-ui';
  ctx.textAlign = 'right';
  for (let v of [0, 0.5, 1.0]) {
    const py = PAD.top + PH - v * PH;
    ctx.fillText(v.toFixed(1), PAD.left - 4, py + 4);
    ctx.strokeStyle = '#e3e8ed';
    ctx.setLineDash([3, 3]);
    ctx.beginPath();
    ctx.moveTo(PAD.left, py);
    ctx.lineTo(PAD.left + PW, py);
    ctx.stroke();
    ctx.setLineDash([]);
  }

  // X axis labels
  ctx.textAlign = 'center';
  ctx.fillText('0', PAD.left, PAD.top + PH + 16);
  ctx.fillText('½', PAD.left + PW / 2, PAD.top + PH + 16);
  ctx.fillText('1', PAD.left + PW, PAD.top + PH + 16);

  // Axis labels
  ctx.fillStyle = '#555';
  ctx.font = '11px system-ui';
  ctx.textAlign = 'left';
  ctx.fillText('Temperature', PAD.left + 2, PAD.top - 4);
  ctx.textAlign = 'center';
  ctx.fillText('Position along bar', PAD.left + PW / 2, CH - 2);

  // Compute t (dimensionless time = step * r / (N-1)^2 * (N-1)^2 * ...)
  // In our scheme: t_physical = step * dt. We set dx = 1/(N-1), so
  // r = alpha*dt/dx^2 => dt = r*dx^2/alpha. With alpha=1 and dx=1/(N-1):
  const t = step * r / ((N - 1) * (N - 1));

  // Draw exact solution (blue)
  const ex = exactU(t);
  ctx.strokeStyle = '#1a6fc4';
  ctx.lineWidth = 1.8;
  ctx.setLineDash([5, 3]);
  ctx.beginPath();
  for (let i = 0; i < N; i++) {
    const px = PAD.left + (i / (N - 1)) * PW;
    const py = PAD.top  + PH - ex[i] * PH;
    i === 0 ? ctx.moveTo(px, py) : ctx.lineTo(px, py);
  }
  ctx.stroke();
  ctx.setLineDash([]);

  // Draw FD approximation (orange/red)
  ctx.strokeStyle = '#e05d2a';
  ctx.lineWidth = 2.2;
  ctx.beginPath();
  for (let i = 0; i < N; i++) {
    const px = PAD.left + (i / (N - 1)) * PW;
    const py = PAD.top  + PH - Math.max(0, Math.min(1, u[i])) * PH;
    i === 0 ? ctx.moveTo(px, py) : ctx.lineTo(px, py);
  }
  ctx.stroke();

  // Dots for FD
  ctx.fillStyle = '#e05d2a';
  for (let i = 0; i < N; i++) {
    const px = PAD.left + (i / (N - 1)) * PW;
    const py = PAD.top  + PH - Math.max(0, Math.min(1, u[i])) * PH;
    ctx.beginPath();
    ctx.arc(px, py, 2.4, 0, 2 * Math.PI);
    ctx.fill();
  }

  // Legend
  ctx.font = '11px system-ui';
  ctx.textAlign = 'left';
  ctx.fillStyle = '#e05d2a';
  ctx.fillText('▬ FD grid', PAD.left + 4, PAD.top + 14);
  ctx.fillStyle = '#1a6fc4';
  ctx.fillText('╌ Exact', PAD.left + 78, PAD.top + 14);

  // Info
  const stable = r <= 0.5;
  const infoEl = document.getElementById('info');
  infoEl.textContent = fmt('Step %step%  —  r = %r%  (%v%)', { step: step, r: r.toFixed(2), v: stable ? 'stable ✓' : 'UNSTABLE ✗ — expect blow-up!' });
  infoEl.style.color = stable ? '#555' : '#c92f3c';
}

// --- Controls ---
const nSlider = document.getElementById('nSlider');
const rSlider = document.getElementById('rSlider');

nSlider.addEventListener('input', () => {
  N = +nSlider.value;
  document.getElementById('nVal').textContent = N;
  stopRun(); init(); draw();
});
rSlider.addEventListener('input', () => {
  r = +rSlider.value;
  document.getElementById('rVal').textContent = r.toFixed(2);
  draw();
});

document.getElementById('btnStep').addEventListener('click', () => {
  fdStep(); draw();
});
document.getElementById('btnRun').addEventListener('click', startRun);
document.getElementById('btnStop').addEventListener('click', stopRun);
document.getElementById('btnReset').addEventListener('click', () => {
  stopRun(); init(); draw();
});

function startRun() {
  if (timer) return;
  document.getElementById('btnRun').disabled  = true;
  document.getElementById('btnStop').disabled = false;
  timer = setInterval(() => {
    for (let k = 0; k < 3; k++) fdStep();
    draw();
  }, 50);
}
function stopRun() {
  clearInterval(timer); timer = null;
  document.getElementById('btnRun').disabled  = false;
  document.getElementById('btnStop').disabled = true;
}

// --- Boot ---
init();
draw();
Result
Console

The bar is divided into N grid points. At each time step the temperature at an interior point is updated by the rule:

ui←ui+r (ui+1−2ui+ui−1)u_i \leftarrow u_i + r\,(u_{i+1} - 2u_i + u_{i-1})

where r=α dt/dx2r = \alpha \, dt/dx^2 is the stability parameter. Keep r≤0.5r \le 0.5 and the scheme is stable; push rr above 0.50.5 and errors explode. Try the Step button to advance one tick, or Run to watch diffusion unfold in real time. The exact solution (computed with the first 20 sine terms) is drawn in blue for comparison — notice how the grid approximation tracks it closely but not perfectly.

Error, Stability, and Cost

Finite differences are a bargain: you give up exactness and get a method that runs in time proportional to the number of grid points. Here is what that bargain looks like in detail.

Truncation error. Replacing f′(x)f'(x) with f(x+h)−f(x)h\frac{f(x+h)-f(x)}{h} is only exact if h=0h = 0. Taylor's theorem tells us the error is proportional to hh — first-order accuracy. The symmetric central difference f(x+h)−f(x−h)2h\frac{f(x+h)-f(x-h)}{2h} achieves second-order accuracy: halving hh cuts the error by four. Most production codes use second- or fourth-order stencils for this reason.

Stability. For the heat equation the explicit scheme is stable only when r=α dt/dx2≤12r = \alpha \, dt/dx^2 \le \tfrac{1}{2}. This is an instance of the CFL condition (Courant–Friedrichs–Lewy, 1928): information cannot travel more than one grid cell per time step. Halving dxdx forces dtdt to shrink by four to keep rr in bounds — one reason implicit schemes (which solve a linear system at each step but allow large dt) dominate in practice.

Cost. Updating N grid points for T time steps costs O(N⋅T)O(N \cdot T) work. For 2-D or 3-D problems the grid becomes N2N^{2} or N3N^{3} cells, and the cost grows accordingly. Techniques like adaptive mesh refinement and multigrid solvers tame these costs for real-world problems.

The method is firmly in P: solved (polynomial time) once you fix the grid and accuracy target. The hard questions are engineering ones — choosing the right stencil, time-stepper, and resolution — not fundamental computational barriers like those facing P vs NP.

Where It Matters

"Turn this differential equation into arithmetic" is one of the most useful things a computer ever does, and finite differences are the simplest way to do it:

  • Weather forecasting: Richardson's 1922 hand calculation used finite differences; every operational weather model today is a descendant, now running on petaflop machines.
  • Chip thermal analysis: engineers simulate heat spreading through a processor die to catch hot spots before fabrication — the same 1-D or 2-D heat equation from the demo above, just on a finer grid.
  • Computational fluid dynamics: the Navier–Stokes equations that govern airflow around a wing or water through a pipe are discretized on a grid. Finite differences and their close cousin finite volumes are two of the dominant approaches.
  • Image processing: the image gradient (edge detection) is a finite difference in 2-D; the anisotropic diffusion filter that smooths photos while preserving edges uses the same heat-equation stepping.
  • Options pricing: the Black–Scholes PDE that prices financial derivatives is solved numerically on a price–time grid using finite differences.

Anywhere a law of nature or a financial model is expressed as a differential equation, finite differences are often the first tool engineers and scientists reach for — simple to implement, easy to reason about, and fast enough for millions of grid points on modern hardware.

Conclusion

The finite difference method is an act of productive compromise: instead of demanding an infinitely thin gap, keep the gap small and do the arithmetic. The result is not exact — but it is computable, and often accurate enough for the task at hand.

Richardson dreamed of a "forecast factory" staffed by 64,000 human computers working in parallel; today a single chip runs the same discretized equations in real time. The idea is still the same: divide space into a grid, replace derivatives with differences, and march forward one step at a time.

If you want to go further, explore how linear programming solves large systems that arise from implicit finite difference schemes, or how matrix multiplication efficiency shapes the speed of every PDE solver. The finite difference method is where calculus meets computer science — and the meeting point is just subtraction.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/finite-difference/Content licensed under CC BY-NC 4.0.