Numerical analysis

A time-stepping method approximates the evolution of an ordinary differential equation by successive updates at discrete times. It can be applied after a spatial finite element method has converted a time-dependent partial differential equation into an ordinary differential equation system. The implicit time-stepping method is one class of such methods.

Implicit time-stepping method

An implicit time-stepping method determines a new solution value through an equation involving that same value. A small-step contraction mapping argument often establishes local solvability when the vector field is Lipschitz continuous. Examples include the Backward Euler method, implicit Runge-Kutta methods and implicit linear multistep methods. Solvability of the update is a separate issue from stability of a numerical method.

Numerical differentiation ( $f^{'} (x) \approx \frac{f ( x + h ) - f ( x )}{h}$ )

Numerical differentiation estimates derivatives from finitely many function values, often by differentiating an interpolation polynomial or using a finite difference. Its error combines truncation error from the formula with errors already present in the sampled values and additional roundoff error. Decreasing the sample spacing can decrease truncation error while increasing the effect of uncertain values. This tradeoff is explained by the ill-posedness of numerical differentiation when only an absolute value-error bound is available.

Ill-posedness of numerical differentiation ( $∥ f_{k} ∥_{\infty} \to 0 need not imply ∥ f_{k}^{'} ∥_{\infty} \to 0$ )

Recovering a derivative from values controlled only in the uniform norm is an ill-posed problem: the differentiation map is not continuous in those norms. On

[0, 1]

, take

f_{k} (x) = sin (k x) / k

. Then

∥ f_{k} ∥_{\infty} \leq 1/ k \to 0

, while

f_{k}^{'} (x) = cos (k x)

has uniform norm one for every

k

. Arbitrarily small value perturbations can therefore conceal nonvanishing derivative changes.

Numerical differentiation can still be useful when additional smoothness bounds, a suitable sampling scale, or a regularization method control the admissible changes. The basic noise amplification by differentiation of a linear interpolant gives a finite-data version: endpoint uncertainty of size

δ

produces derivative uncertainty of order

δ / h

at spacing

h

Noise amplification by differentiation of a linear interpolant ( $Δ p (x) = δ x / h, Δ p^{'} (x) = δ / h$ )

For linear interpolation at zero and

h > 0

, an error

δ

in the second sampled value, with the first exact, changes the linear interpolant by

Δ p (x) = δ x / h

. On

[0, h]

, its uniform norm is exactly

∣ δ ∣

, whereas the derivative change is the constant

δ / h

. With independent absolute bounds

∣ δ_{0} ∣, ∣ δ_{1} ∣ \leq δ

on both samples, the value change remains at most

δ

and the derivative change is at most

2 δ / h

, attained for opposite endpoint errors.

Thus bounded absolute sample noise can remain bounded under value interpolation while becoming unbounded under numerical differentiation as

h \to 0

. Fixed nonzero sample noise also prevents exact convergence of the value approximation; stability against perturbations and convergence to exact data are different assertions. For a second-sample error

a ϵ (h)

with fixed

a > 0

, value convergence requires

ϵ (h) \to 0

, while derivative convergence requires

ϵ (h) = o (h)

in little o notation, in addition to vanishing truncation error.

Discretization

 2  0

Replacing continuous space, time or fields by finite degrees of freedom permits a numerical forecast. Approximation and truncation errors depend on resolution, algorithms and unresolved scales; they are distinct from observation errors and from roundoff. Refinement reduces a consistent truncation error only when stability and the relevant smoothness assumptions permit convergence.

Implicit numerical method

A method defining the next state through equations involving that unknown state. Solving these equations can be costly, while enlarged stability regions are useful for stiff problems.

Explicit numerical method

A method whose next state is obtained from known stages without solving simultaneous equations involving the unknown next state. Stability may restrict its admissible step sizes.

One-step method ( $y_{n + 1} = y_{n} + h Φ (x_{n}, y_{n}, h)$ )

A one-step method approximates an ordinary differential equation using the current approximate solution and evaluations within the next step. The function

Φ

is its increment function. A Runge-Kutta method may use several intermediate stages while remaining a one-step method: it does not require additional previous solution values as a linear multistep method does. Its consistency error is the exact-start update defect divided by

h

Runge-Kutta method

A Runge-Kutta method advances an initial value problem

y^{'} = f (t, y)

y (t_{0}) = y_{0}

, by sampling its vector field at several trial states within one time step. These samples are Runge-Kutta stages; a weighted combination of their slopes produces the next state. The construction improves on the single tangent line used by the Euler method while requiring evaluations of

f

, rather than explicit higher derivatives of the solution. It applies to scalar equations and systems of ordinary differential equations alike. The coefficient choices distinguish methods such as Heun's method, the classical fourth-order Runge-Kutta method and implicit Runge-Kutta methods.

For

s

stages, write

A = (a_{ij})

for the stage matrix,

c_{i}

for relative time nodes and

b_{i}

for output weights. A time step of length

h

computes

Y_{i} = y_{n} + h \sum_{j = 1}^{s} a_{ij} k_{j}, k_{i} = f (t_{n} + c_{i} h, Y_{i}), y_{n + 1} = y_{n} + h \sum_{i = 1}^{s} b_{i} k_{i} .

Here

k_{i}

is a slope, with units of state per time;

h k_{i}

is a state increment. Another common convention writes

K_{i} = h k_{i}

, so its stage states are

Y_{i} = y_{n} + \sum_{j} a_{ij} K_{j}

and its output is

y_{n} + \sum_{i} b_{i} K_{i}

. These are equivalent conventions: including

h

both inside

K_{i}

and again in the output would count the time step twice. The textbook formulation uses the increment convention.^[ref]

Usually the Butcher tableau obeys

c_{i} = \sum_{j} a_{ij}

, called internal consistency. Stage

i

then combines slopes over an effective time

c_{i} h

. The second-order conditions with independent Runge-Kutta abscissae explain what changes when time nodes are specified independently. A Runge-Kutta stage is generally an auxiliary approximation, not a separate accepted solution value at its node. Several stages can even have the same time node but different state arguments, as in the classical fourth-order Runge-Kutta method.

The useful idea is to approximate the integral

y (t_{n} + h) - y (t_{n}) = \int_{t_{n}}^{t_{n} + h} f (t, y (t)) d t

without knowing the intermediate exact trajectory. Trial states provide places to sample

f

, and Butcher order conditions make the resulting combination agree with the exact Taylor expansion. More stages alone do not guarantee a higher order of a numerical method: their coefficients must satisfy the appropriate conditions. Embedded Runge-Kutta pairs reuse these same samples for a second output and a practical local error estimator.

Runge-Kutta stage ( $k_{i} = f (t_{n} + c_{i} h, Y_{i})$ )

A Runge-Kutta stage consists of a trial state

Y_{i}

at the time node

t_{n} + c_{i} h

and the corresponding vector field evaluation

k_{i}

. In the slope convention,

Y_{i} = y_{n} + h \sum_{j} a_{ij} k_{j}

. Stages belong to the construction of one Runge-Kutta method step; they need not have the same accuracy as its final output. Stage states can leave the physically allowed region even when the endpoint lies inside it, so

f

must be defined at every stage used.

Explicit Runge-Kutta method

An explicit Runge-Kutta method has

a_{ij} = 0

whenever

j \geq i

. Its stage matrix is strictly lower triangular: compute each Runge-Kutta stage using already known slopes, starting with

Y_{1} = y_{n}

. In contrast, an implicit Runge-Kutta method requires solving coupled stage equations, whose existence is addressed by stage solvability of an implicit Runge-Kutta method.

Heun's method ( $y_{n + 1} = y_{n} + h (k_{1} + k_{2}) /2$ )

Here Heun's method means the explicit endpoint-averaging Runge-Kutta method: take an Euler predictor

Y_{2} = y_{n} + h f (t_{n}, y_{n})

, sample

k_{2} = f (t_{n} + h, Y_{2})

there, and average it with

k_{1} = f (t_{n}, y_{n})

for the output. The second Runge-Kutta stage is evaluated at the predictor, not at the corrected endpoint. In particular it requires no implicit solve as the ODE trapezoidal rule does. Naming conventions vary between books; the formulas specify the method unambiguously.^[ref]

For smooth

f

, a Taylor expansion about

(t_{n}, y_{n})

gives

k_{2} = f + h (f_{t} + f_{y} f) + O (h^{2})

, where

f_{y}

is the state Jacobian matrix. Consequently the output is

y_{n} + h f + h^{2} (f_{t} + f_{y} f) /2 + O (h^{3})

, matching the exact solution through degree two by the chain rule. This gives local truncation error

O (h^{3})

and, with the usual bounded smoothness and Lipschitz continuity hypotheses on a fixed finite interval, global error

O (h^{2})

. The Euler-Heun embedded pair compares this output with its Euler predictor.

Figure 1.
Euler-Heun stage geometry for y prime equals y: sample at the Euler predictor, then average slopes for the corrected endpoint
. **Details**.

For

y^{'} = y

y_{n} = 1

and

h = 0.2

, the sampled slopes are

1

and

1.2

. Transporting the second slope to the starting point would end at

1.24

; the average slope

1.1

instead ends at

1.22

. The exact exponential function reaches

e^{0.2} \approx 1.22140275816

. The magnified endpoint comparison separates the pair difference

0.02

from the actual Heun local defect

0.00140275816

Maximal-order Runge-Kutta methods are A-stable ( $p = 2 s ⟹ A -stability$ )

s

-stage Runge-Kutta method has stability function

R (z) = 1 + z b^{T} (I - z A)^{- 1} 1

, a rational function with numerator and denominator degrees at most

s

. Order

2 s

forces

R (z) - e^{z} = O (z^{2 s + 1})

, so uniqueness of the diagonal Padé approximant identifies

R

with

[s / s]_{e^{z}}

. That approximant has no poles in the closed left half-plane and modulus at most one there. Hence every such maximal-order method is A-stable. Since the Padé approximant numerator and denominator are coprime and both have degree

s

, the original stage determinant cannot have hidden canceled poles in that half-plane.

Theta method ( $y_{n + 1} - y_{n} = h [(1 - a) f (y_{n}) + a f (y_{n + 1})]$ )

This one-step family interpolates between the Forward Euler method at

a = 0

, the trapezoidal rule at

a = 1/2

, and the Backward Euler method at

a = 1

. Its stability function is

R (z) = [1 + (1 - a) z] / [1 - a z]

. It is first order except at one half, where it is second order. Applying it to spatial diffusion gives the theta diffusion stability criterion.

Theta diffusion stability criterion ( $G (s) = [1 - (1 - a) s] / [1 + a s], s \geq 0$ )

On an orthogonal diffusion eigenmode, the theta method has the displayed amplification with

s = - Δ t λ

. The identity

(1 + a s)^{2} - [1 - (1 - a) s]^{2} = 2 s + (2 a - 1) s^{2}

proves contraction at all nonnegative

s

exactly when

a \geq 1/2

. For

a < 1/2

, contraction holds up to

s \leq 2/ (1 - 2 a)

, with a positive denominator in that range. Orthogonality is what transfers these modal bounds to a full L2 norm estimate.

Rough-data damping distinction for the theta method ( $G (s) \to - (1 - a) / a (s \to \infty, a > 0)$ )

a = 1/2

, arbitrarily stiff heat modes approach amplification

- 1

: the Crank-Nicolson method is contractive but does not rapidly damp unresolved rough components. The Backward Euler method instead has limit zero. Thus unconditional linear stability and strong stiff decay are different requirements. High smooth-data consistency order alone does not imply a uniform high-order estimate down to the initial instant for arbitrary L2 initial conditions.

Strong stability preserving Runge-Kutta method

A strong stability preserving Runge-Kutta method writes its stages as convex combinations of previously computed states and suitable Forward Euler method steps. It transfers any convex-functional nonincrease property of the forward step, such as a norm bound, under a proportionally scaled time-step restriction.

For example, if

E_{k} (U) = U + k F (U)

is nonexpansive in a norm, then

Y = E_{k} (U)

and

U_{new} = (U + E_{k} (Y)) /2

is also nonexpansive: the triangle inequality bounds the new difference by one half of the initial difference plus one half of the twice-advanced difference. Both are at most the initial difference. A Taylor expansion yields

U_{new} = U + k F (U) + k^{2} F^{'} (U) F (U) /2 + O (k^{3})

, so this is a second-order Runge-Kutta method. The forward-step hypothesis must hold on the stage states; preservation of one chosen convex bound is not the same as unconditional B-stability.

Classical fourth-order Runge-Kutta method ( $R (z) = 1 + z + z^{2} /2 + z^{3} /6 + z^{4} /24$ )

The classical four-stage explicit Runge-Kutta method uses

A = 0 1/2 00 00 1/2 0 00010000, b = (1/6, 1/3, 1/3, 1/6)^{T}, c = (0, 1/2, 1/2, 1)^{T} .

It satisfies the fourth-order conditions for a Runge-Kutta method. Its polynomial stability function excludes A-stability, while

∣ R (i q) ∣^{2} = 1 - q^{6} /72 + q^{8} /576

gives the imaginary-axis stability interval

∣ q ∣ \leq 22

Butcher tableau

A Butcher tableau records a Runge-Kutta method's stage matrix

A

, time nodes

c

, and output weights

b

in the block array

c A b^{T} .

Read row

i

above the line as the recipe for Runge-Kutta stage

i

: the left entry gives

t_{n} + c_{i} h

, and the right entries give the coefficients in

Y_{i} = y_{n} + h \sum_{j} a_{ij} k_{j}

. Below the line,

b_{i}

weights slope

k_{i}

in the final update. Zeros on and above the diagonal identify an Explicit Runge-Kutta method. For the internally consistent convention, row sums of

A

equal

c

; first-order consistency of a numerical method also requires

\sum_{i} b_{i} = 1

. An embedded Runge-Kutta pair adds another weight row while retaining

A, c

Butcher order condition

A Butcher order condition equates the coefficient assigned by a Runge-Kutta method to a rooted tree with the corresponding exact Taylor coefficient. Satisfying every condition through tree order

p

is equivalent to Runge–Kutta order at least

p

Second-order conditions with independent Runge-Kutta abscissae ( $b^{T} e = 1, b^{T} c = b^{T} A e = 1/2$ )

For a general nonautonomous Runge-Kutta method, Taylor expansion of stage

i

gives

k_{i} = f + h [c_{i} f_{t} + (A e)_{i} f_{y} f] + O (h^{2})

. Matching the exact increment

h f + h^{2} (f_{t} + f_{y} f) /2

gives the three displayed conditions. When internal consistency

c = A e

is imposed, the last two coincide. Internal consistency is a usual tableau convention, but is not logically necessary for order two of a generalized tableau.

Fourth-order conditions for a Runge-Kutta method

With

c = A e

and

C = diag (c)

, a Runge-Kutta method has order at least four exactly when the eight Butcher order conditions

b^{T} e = 1, b^{T} c = \frac{1}{2}, b^{T} c^{2} = \frac{1}{3}, b^{T} A c = \frac{1}{6}, b^{T} c^{3} = \frac{1}{4}, b^{T} C A c = \frac{1}{8}, b^{T} A c^{2} = \frac{1}{12}, b^{T} A^{2} c = \frac{1}{24}

hold; powers of

c

are componentwise. These conditions match the numerical and exact coefficients indexed by rooted trees through order four. Matching the scalar Dahlquist test equation alone is insufficient to check the nonlinear Butcher order conditions.

Implicit Runge-Kutta method

An implicit Runge-Kutta method defines one or more stages through equations involving those same stages.

Two-stage Runge-Kutta family with an implicit second stage ( $R (z) = \frac{1 + ( 1 - a ) z + ( 1/2 - a ) z ^{2}}{1 - a z}$ )

For stages

k_{1} = f (y)

and

k_{2} = f (y + (1 - a) h k_{1} + ah k_{2})

with update

y + h (k_{1} + k_{2}) /2

, the method has order two for every fixed real

a

. The third-order coefficient of

f^{''} f^{2}

is always

1/4

, not

1/6

, so no parameter gives third order on general nonlinear equations. The displayed stability function is bounded throughout the left half-plane only when

a = 1/2

, when it is the A-stable trapezoidal function

(1 + z /2) / (1 - z /2)

. It is not L-stable.

Stage solvability of an implicit Runge-Kutta method

For a vector field uniformly Lipschitz continuous in the state with constant

L

, the stage fixed-point map of an implicit Runge-Kutta method is a contraction in the maximum stage norm if

h L max_{i} \sum_{j} ∣ a_{ij} ∣ < 1

. The contraction mapping theorem then gives unique stages. This is a sufficient small-step condition; a dissipative vector field may permit solvability for larger steps. An algebraic stability or B-stability estimate compares existing stage solutions and should not be mistaken for an unqualified stage-existence theorem.

Implicit midpoint rule

The implicit midpoint rule advances an autonomous ordinary differential equation by

y_{n + 1} = y_{n} + h f (\frac{y _{n} + y _{n + 1}}{2}) .

It is a second-order, A-stable, time-symmetric numerical method.

The implicit midpoint rule is the midpoint method with

y_{n + 1} = y_{n} + h f (t_{n} + h /2, (y_{n} + y_{n + 1}) /2)

. It is a one-stage implicit Runge-Kutta method.

Implicit midpoint Lipschitz estimate ( $Lip (Ψ_{h}) \leq (1 + h L /2) / (1 - h L /2)$ )

For a uniformly Lipschitz continuous state field

f

with constant

L

, an implicit midpoint rule step is uniquely solvable when

h L /2 < 1

, by the contraction mapping theorem. Write its step map as

Ψ_{h} (u) = u + h Φ_{h} (u)

, where

Φ_{h} (u) = f (t + h /2, u + h Φ_{h} (u) /2)

. Subtracting the increment equations gives

(1 - h L /2) ∥ Φ_{h} (u) - Φ_{h} (v) ∥ \leq L ∥ u - v ∥.

Therefore

Lip (Φ_{h}) \leq L / (1 - h L /2)

, and the triangle inequality gives the bound for

Ψ_{h}

. Combining it with the implicit midpoint consistency bound yields

∥ e_{n} ∥ \leq R_{h}^{n} ∥ e_{0} ∥ + \frac{K h ^{2}}{L} (R_{h}^{n} - 1), R_{h} = \frac{1 + h L /2}{1 - h L /2}

for

L > 0

. This follows by summing the geometric error recurrence with forcing

K h^{3} / (1 - h L /2)

. For

h L \leq 1

and

nh \leq T

R_{h}^{n} \leq e^{2 L T}

, proving second-order global error with exact starting data. If

L = 0

, the direct bound is

∥ e_{n} ∥ \leq ∥ e_{0} ∥ + K T h^{2}

Implicit midpoint consistency bound ( $∥ τ_{n} ∥ \leq (M_{3} /24 + L M_{2} /8) h^{2}$ )

For an exact solution

y^{'} = f (t, y)

, define the normalized local truncation error of the implicit midpoint rule by

τ_{n} = \frac{y ( t _{n} + h ) - y ( t _{n} )}{h} - f (t_{n} + \frac{h}{2}, \frac{y ( t _{n} + h ) + y ( t _{n} )}{2}) .

∥ y^{''} ∥ \leq M_{2}

∥ y^{'''} ∥ \leq M_{3}

and the state Jacobian matrix of

f

has operator norm at most

L

on the intervening state segments, then the displayed bound holds uniformly on a fixed finite interval. To prove it, put

t = t_{n} + h /2

. Symmetric Taylor expansions give

\frac{y ( t + h /2 ) - y ( t - h /2 )}{h} - y^{'} (t) \leq \frac{M _{3} h ^{2}}{24}, \frac{y ( t + h /2 ) + y ( t - h /2 )}{2} - y (t) \leq \frac{M _{2} h ^{2}}{8} .

For the first bound, integrate the second-order remainder of

y^{'} (t + s)

over the symmetric step; for the second, average the second-order remainders of

y

at the endpoints. The odd terms cancel. The Lipschitz condition for

f

and

y^{'} (t) = f (t, y (t))

finish the proof by the triangle inequality. The unnormalized step residual is

h τ_{n} = O (h^{3})

Crank-Nicolson method

The Crank--Nicolson method applies the trapezoidal rule to a semidiscrete evolution equation. For a linear system

u^{'} = Lu

its update is

(I - \frac{h}{2} L) u^{n + 1} = (I + \frac{h}{2} L) u^{n} .

Contractivity of split Crank-Nicolson diffusion ( $∥ C_{A} (h) C_{B} (h) ∥_{2} \leq 1$ )

For symmetric negative semidefinite matrices

A, B

, each factor

C_{A} (h) = (I - h A /2)^{- 1} (I + h A /2)

has norm at most one, since its eigenvalues are

(1 + hλ /2) / (1 - hλ /2)

with

λ \leq 0

. The same holds for

B

, so their ordered product is a contraction for every

h \geq 0

, whether or not

A, B

commute. This proves unconditional discrete L2 norm stability of the split Crank-Nicolson method for directional diffusion. Its splitting error remains first order in general, because the product's quadratic defect is

h^{2} [A, B] /2

Collocation Runge-Kutta method

Given distinct nodes

c_{i}

, a collocation Runge-Kutta method uses the Lagrange basis

ℓ_{j}

with

a_{ij} = \int_{0}^{c_{i}} ℓ_{j} (τ) d τ

and

b_{j} = \int_{0}^{1} ℓ_{j} (τ) d τ

. Its stages enforce the differential equation at the collocation nodes.

This applies the collocation method to the time-step polynomial of an ordinary differential equation.

Two-node collocation A-stability criterion ( $c_{1} + c_{2} \geq 1$ )

For distinct nodes

c_{1}, c_{2} \in [0, 1]

, the two-stage collocation Runge-Kutta method has stability function

R (z) = \frac{1 + ( 1 - a ) z + r z ^{2}}{1 - a z + d z ^{2}}, a = (c_{1} + c_{2}) /2, d = c_{1} c_{2} /2, r = (1 - c_{1}) (1 - c_{2}) /2.

The denominator has no zeros in the closed left half-plane. On the imaginary axis the denominator's squared modulus minus the numerator's is

(d^{2} - r^{2}) y^{4}

. Thus

d \geq r

, equivalently

c_{1} + c_{2} \geq 1

, is necessary and sufficient for A-stability, using the maximum modulus principle and the bound at infinity.

Collocation order theorem ( $p = s + k$ )

For

s

distinct collocation nodes, let

π_{s} (t) = \prod_{i = 1}^{s} (t - c_{i})

. If

\int_{0}^{1} t^{j} π_{s} (t) d t = 0

for

0 \leq j < k

, and the next moment is nonzero, the associated interpolatory quadrature has order

s + k

, and the collocation Runge-Kutta method has the same classical order for sufficiently smooth differential equations. Here

0 \leq k \leq s

. Equivalently, if the quadrature is exact through degree

p - 1

but not degree

p

, the method has order

p

. The generic order is at least

s

, Gauss collocation methods reach

2 s

, and Radau IIA methods reach

2 s - 1

Collocation endpoint superconvergence by quadrature orthogonality ( $p = q$ )

For

s

collocation nodes, the residual of the collocation polynomial vanishes at each node and contains the factor

ω (τ) = \prod_{i} (τ - c_{i})

. Quadrature exactness through degree

q - 1

implies

\int_{0}^{1} ω (τ) p (τ) d τ = 0

for every polynomial of degree at most

q - s - 1

. Expand the smooth vector-field interpolation error and the derivative of the exact flow in powers of the step. Every endpoint-error term of degree below

q

has this vanishing weighted moment, so the endpoint local error is

O (h^{q + 1})

, although the interior stage accuracy is generally only

s

. A step Lipschitz estimate and discrete Gronwall inequality then give global order

q

. Gaussian nodes attain

q = 2 s

, while the first nonexact quadrature moment bounds the order from above.

Collocation tableau row-sum identity ( $A 1 = c$ )

For a collocation Runge-Kutta method,

a_{ij} = \int_{0}^{c_{i}} ℓ_{j} (t) d t

, where the Lagrange interpolation basis obeys

\sum_{j} ℓ_{j} = 1

. Hence each row sum is its collocation node. Failure of this identity disproves a claimed collocation tableau with those nodes. The identity alone is necessary rather than sufficient: the individual integrated basis coefficients must also agree.

Lobatto IIIA method

A Lobatto IIIA method collocates at the Lobatto quadrature nodes, including both endpoints of the time interval. Its three-stage version has

A = 0 5/24 1/6 0 1/3 2/3 0 - 1/24 1/6, b = (1/6, 2/3, 1/6)^{T}, c = (0, 1/2, 1)^{T} .

The Butcher order conditions give order four. Its stability function is

R (z) = \frac{1 + z /2 + z ^{2} /12}{1 - z /2 + z ^{2} /12},

which is A-stable but not L-stable. It is not algebraically stable: the first diagonal entry of the defining matrix is

2 b_{1} a_{11} - b_{1}^{2} = - 1/36

. This distinguishes linear A-stability from nonlinear algebraic stability.

Gauss-Legendre method

The

s

-stage Gauss--Legendre Runge--Kutta method collocates at the Gauss--Legendre quadrature nodes on

[0, 1]

. It is symmetric, A-stable, and has order

2 s

Gauss collocation coefficient construction ( $a_{ij} = \int_{0}^{c_{i}} ℓ_{j} (t) d t$ )

Choose

c_{1}, \dots, c_{s}

as the polynomial roots of

P_{s} (2 c - 1)

(0, 1)

, and form the Lagrange interpolation polynomials

ℓ_{j} (t) = \prod_{m \neq = j} (t - c_{m}) / (c_{j} - c_{m})

. The Gauss collocation method has

a_{ij} = \int_{0}^{c_{i}} ℓ_{j} (t) d t

and

b_{j} = \int_{0}^{1} ℓ_{j} (t) d t

. This gives explicit coefficient formulas after finding polynomial roots: if

ℓ_{j} (t) = \sum_{r} d_{j r} t^{r}

, these integrals are

\sum_{r} d_{j r} c_{i}^{r + 1} / (r + 1)

and

\sum_{r} d_{j r} / (r + 1)

. Gaussian quadrature gives order

2 s

Fourth-order two-stage Gauss collocation method ( $R (z) = (1 + z /2 + z^{2} /12) / (1 - z /2 + z^{2} /12)$ )

The two-stage Gauss--Legendre Runge-Kutta method uses nodes

c_{1, 2} = 1/2 \mp 3 /6

, weights

b_{1} = b_{2} = 1/2

and coefficients

a_{11} = a_{22} = 1/4

a_{12} = 1/4 - 3 /6

a_{21} = 1/4 + 3 /6

. The Butcher order conditions give order four. Its stability function is the displayed rational function, which is A-stable. The algebraic stability matrix is zero, so positive weights also give nonlinear contractivity for dissipative vector fields.

Radau IIA method

s

-stage Radau IIA method collocates at the right-endpoint Radau nodes. It has order

2 s - 1

and is both algebraically stable and L-stable.

Padé identification of a Radau stability function ( $R (z) = [s - 1/ s]_{e^{z}}$ )

Suppose an

s

-stage Runge-Kutta method has order

2 s - 1

, an invertible matrix

A

, and

b^{T} A^{- 1} 1 = 1

. Its stability function is

R (z) = 1 + z b^{T} (I - z A)^{- 1} 1

. The adjugate matrix formula gives a rational function with denominator degree

s

and numerator degree at most

s

; its zero limit at infinity lowers the latter degree to at most

s - 1

. The order of a numerical method gives

R (z) - e^{z} = O (z^{2 s})

, identifying the

[s - 1/ s]

Padé approximant. Uniqueness follows because the cross-multiplied difference of two candidates has degree at most

2 s - 1

but a zero of order at least

2 s

, so is identically zero. The A-stability of near-diagonal exponential Padé approximants then proves A-stability, and the zero limit proves L-stability. These conditions identify the stability function; they do not alone claim uniqueness of all entries of the Runge-Kutta method tableau.

Algebraic stability of a Runge-Kutta method

A Runge-Kutta method is algebraically stable when

b_{i} \geq 0

and the matrix

M

with

M_{ij} = b_{i} a_{ij} + b_{j} a_{ji} - b_{i} b_{j}

is positive semidefinite. Algebraic stability implies contractivity for dissipative differential equations.

Runge-Kutta conservation of quadratic invariants ( $b_{i} a_{ij} + b_{j} a_{ji} = b_{i} b_{j}$ )

For stage states

Y_{i} = y_{n} + h \sum_{j} a_{ij} k_{j}

and a vector field with

Y_{i}^{T} k_{i} = 0

, expansion of the new squared norm gives

∥ y_{n + 1} ∥^{2} - ∥ y_{n} ∥^{2} = - h^{2} \sum_{i, j} M_{ij} k_{i}^{T} k_{j}

, where

M_{ij} = b_{i} a_{ij} + b_{j} a_{ji} - b_{i} b_{j}

. Hence

M = 0

preserves the quadratic form exactly for every well-defined step. The same argument with a fixed symmetric matrix in the inner product preserves any quadratic invariant of the vector field. In particular it preserves the Euclidean norm for a skew-symmetric matrix differential equation.

Butcher contractivity theorem

An algebraically stable Runge-Kutta method is B-stable: for a dissipative vector field, a well-defined step cannot increase the distance between two starting values. This assumes that the compared stage equations have solutions, as required by stage solvability of an implicit Runge-Kutta method. The Runge-Kutta contractivity identity expresses the distance change as a nonpositive weighted term coming from the dissipative vector field minus a nonnegative quadratic form involving the matrix defining algebraic stability. Positive weights also rule out stage poles in the open left half-plane on the scalar test equation, giving A-stability there; zero-weight redundant stages need separate solvability attention.

Runge-Kutta contractivity identity

For differences

D_{i}

between corresponding stages and differences

F_{i}

between their vector fields, a Runge--Kutta step satisfies

∥ d_{n + 1} ∥^{2} = ∥ d_{n} ∥^{2} + 2 h \sum_{i} b_{i} Re ⟨ D_{i}, F_{i} ⟩ - h^{2} \sum_{i, j} m_{ij} Re ⟨ F_{i}, F_{j} ⟩ .

To derive the identity, expand the squared norm of

d_{n + 1} = d_{n} + h \sum_{i} b_{i} F_{i}

and substitute

d_{n} = D_{i} - h \sum_{j} a_{ij} F_{j}

in the linear term. Symmetrizing the double sum gives the displayed coefficients

m_{ij}

. For a dissipative vector field, the first sum is nonpositive when

b_{i} \geq 0

, and the second sum is nonnegative when the matrix

(m_{ij})

is a positive semidefinite matrix. Thus algebraic stability implies B-stability.

Trapezoidal rule ( $T (f) = (b - a) (f (a) + f (b)) /2$ )

As a quadrature rule, the trapezoidal rule replaces an integrand by the linear interpolation joining its two endpoint values:

I (f) = \int_{a}^{b} f (x) d x \approx T (f) = \frac{b - a}{2} (f (a) + f (b)), a < b .

For nonnegative endpoint values, this is the area of a trapezoid; in general it is a signed area. The quadrature nodes are

a, b

, each with quadrature weight

(b - a) /2

. Unlike the midpoint quadrature rule, it uses endpoint samples, which can be reused between adjacent panels.

Every affine function

f (x) = αx + β

agrees with its linear interpolation, so the trapezoidal rule integrates it exactly. On

[0, 2]

f (x) = 2 x + 1

gives

T (f) = (2/2) (1 + 5) = 6 = I (f)

. On

[0, 1]

, the polynomial

f (x) = x^{2}

instead gives

T (f) = 1/2

, exceeding

I (f) = 1/3

1/6

. Its degree of precision of a quadrature rule is exactly one.

f \in C^{2} [a, b]

, its signed approximation error satisfies

T (f) - I (f) = \frac{( b - a ) ^{3}}{12} f^{''} (ξ_{T}), a < ξ_{T} < b .

A convex function lies below its endpoint chord, so the trapezoidal rule overestimates its integral; a concave function gives the reverse inequality. The error formula assumes that the second derivative exists and is a continuous function throughout the closed interval.^[ref] A short derivation is the identity

T (f) - I (f) = \frac{1}{2} \int_{a}^{b} (t - a) (b - t) f^{''} (t) d t,

obtained by two integrations by parts. The nonnegative Peano kernel has integral

(b - a)^{3} /12

, and the weighted mean value theorem for integrals gives the formula.

The trapezoidal rule for an ordinary differential equation averages vector fields at the old and new states and is A-stable.

This time-stepping use applies the same trapezoidal rule to the integral of the derivative in an ordinary differential equation. That connection explains its place under implicit Runge-Kutta methods; its quadrature behavior and the nonlinear issue in trapezoidal rule fails B-stability concern different uses of the same endpoint averaging.

Sharp one-panel trapezoidal error bound ( $∣ I (f) - T (f) ∣ \leq \frac{( b - a ) ^{3}}{12} max_{[a, b]} ∣ f^{''} ∣$ )

For

f \in C^{2} [a, b]

with

a < b

, the trapezoidal rule satisfies the displayed bound, and the constant

1/12

is best possible. The polynomial interpolation error for the endpoint linear interpolant satisfies

∣ f (x) - p (x) ∣ \leq \frac{1}{2} M (x - a) (b - x)

, where

M = max ∣ f^{''} ∣

. Integrating this inequality and using the triangle inequality for an integral proves the bound because

\int_{a}^{b} (x - a) (b - x) d x = (b - a)^{3} /6

. Equality holds for

f (x) = (x - a)^{2}

: the exact integral is

(b - a)^{3} /3

, its trapezoidal rule estimate is

(b - a)^{3} /2

, and

M = 2

Composite trapezoidal rule ( $T_{n} (f) = h (f (a) /2 + \sum_{j = 1}^{n - 1} f (a + jh) + f (b) /2)$ )

On an equal-panel partition of an interval, with

h = (b - a) / n

, sum the trapezoidal rule over adjacent panels:

T_{n} (f) = h [\frac{f ( a )}{2} + \sum_{j = 1}^{n - 1} f (a + jh) + \frac{f ( b )}{2}] .

There are

n + 1

distinct quadrature nodes. Each interior sample belongs to two panels and therefore has quadrature weight

h

, while each endpoint has quadrature weight

h /2

. Equivalently, the composite trapezoidal rule integrates the piecewise linear function obtained by linear interpolation exactly.^[ref]

For

f \in C^{2} [a, b]

and

K = max_{x \in [a, b]} ∣ f^{''} (x) ∣

, summing panel quadrature errors gives

T_{n} (f) - I (f) = \frac{b - a}{12} h^{2} f^{''} (ξ), ∣ T_{n} (f) - I (f) ∣ \leq \frac{b - a}{12} h^{2} K

for some

ξ \in (a, b)

. Thus panel quadrature errors of order

O (h^{3})

accumulate to order

O (h^{2})

, using big O notation.^[ref] On a nonuniform partition of an interval with widths

h_{j}

and largest width

H

, the bound is

K \sum_{j} h_{j}^{3} /12 \leq K (b - a) H^{2} /12

For the polynomial

f (x) = x^{2}

[0, 1]

, direct summation yields

T_{n} (f) = \frac{1}{3} + \frac{1}{6 n ^{2}}, T_{n} (f) - I (f) = - 2 (M_{n} (f) - I (f)) .

Both quadrature errors shrink by exactly four when

n

doubles, with opposite signs. This exact relation belongs to this polynomial example; the trapezoidal rule is not always twice as inaccurate as the midpoint quadrature rule.

The left plot uses exact polynomial sums; the right adds a common sample bias

δ = 0.002

to the composite trapezoidal rule. That mathematical model converges to

1/3 + δ

, showing why smaller panels alone do not ensure accuracy for noisy samples. It illustrates noise sensitivity of a quadrature rule, rather than measured data.

Symmetric trapezoidal exactness for odd functions ( $T_{N} (g) = \int_{- L}^{L} g (x) d x = 0 (g odd)$ )

On a partition of an interval symmetric about zero, the composite trapezoidal rule gives opposite quadrature nodes equal quadrature weights. Samples of an odd function cancel pairwise, and a zero quadrature node contributes zero. For a continuous odd function, its exact integral over

[- L, L]

is also zero. Adding any affine function preserves exactness by linearity of the integral and of the quadrature rule. For example,

1 + sin x

is nonpolynomial but both the exact integral and every uniform-grid composite trapezoidal rule value on

[- 1, 1]

are two.

Trapezoidal endpoint correction for a right Riemann sum ( $R_{N} - T_{N} = \frac{h}{2} (f (b) - f (a))$ )

On a uniform partition of

[a, b]

with width

h

, the right Riemann sum and the composite trapezoidal rule have equal weights at every interior mesh point. They differ only by adding half a final-endpoint weight and removing half an initial-endpoint weight. This proves the displayed identity. For continuous

f

, both converge to its definite integral; their finite endpoint correction need not be confused with the smaller trapezoidal discretization error for a smooth integrand.

Composite trapezoidal error bound ( $∣ T_{n} - I ∣ \leq \frac{( b - a ) ^{3}}{12 n ^{2}} max_{[a, b]} ∣ f^{''} ∣$ )

f

is twice continuously differentiable on

[a, b]

, take a uniform bound

∣ f^{''} ∣ \leq M

. On one panel

[u, v]

, two integrations by parts give the signed quadrature error

\frac{v - u}{2} (f (u) + f (v)) - \int_{u}^{v} f (t) d t = \frac{1}{2} \int_{u}^{v} (t - u) (v - t) f^{''} (t) d t .

The nonnegative kernel has integral

(v - u)^{3} /6

, so the absolute error is at most

M (v - u)^{3} /12

. Summing

n

equal panels proves the displayed bound. The maximum exists on the closed interval; a supremum or any known uniform bound also suffices. The estimate certifies a sufficient panel count, which need not be the least count achieving a particular actual error.

Trapezoidal rule fails B-stability

The trapezoidal rule is A-stable but is not B-stable. For the dissipative vector field

f (y) = - y^{3}

at step size one, its uniquely defined step map

T

obeys

T + T^{3} /2 = y - y^{3} /2

. At

y = 2

T = 0

and implicit differentiation gives

T^{'} = - 2

, so nearby states expand rather than contract. This is a concrete nonlinear obstruction, rather than merely failure of the sufficient algebraic stability criterion.

Time-symmetric numerical method

A one-step numerical method with step map

Φ_{h}

is time symmetric when

Φ_{- h} = Φ_{h}^{- 1}

. Its leading local-error power is odd, so every nontrivial consistent time-symmetric method has even order.

Order of a Runge-Kutta method

The order is the highest power through which the one-step expansion matches the exact Taylor expansion.

Third-order conditions for an explicit Runge-Kutta method

For an explicit Runge-Kutta method with coefficients

A, b

, put

c = A 1

and let

c^{2}

mean componentwise squaring. The conditions

b^{T} 1 = 1

b^{T} c = 1/2

b^{T} c^{2} = 1/3

and

b^{T} A c = 1/6

match the exact flow expansion through degree three: the two cubic differentials are

f^{''} [f, f]

and

(f^{'})^{2} f

. For sufficiently smooth vector fields these give local error

O (h^{4})

and, under the usual finite-interval Lipschitz continuity assumptions, global error

O (h^{3})

Energy error

The difference between a computed system energy and its exact value. Long-time bounded oscillation of this error is distinct from secular numerical energy drift.

Manufactured solution

A chosen smooth function inserted into a differential equation to define a matching source and boundary data. Its known exact values permit direct verification of an implementation's convergence order.

Error norm

A norm used to quantify numerical error on a grid, such as the maximum norm or a mesh-weighted root-mean-square norm. Convergence studies must compare the same norm at a fixed physical time.

Vectorized numerical computation

Vectorized numerical computation applies numerical operations to entire arrays of values instead of explicitly iterating over scalar entries. For a power series, constructing a vector of term indices, evaluating all corresponding terms and summing that vector is a vectorized evaluation. A loop over bounded blocks can control memory while preserving vectorized computation within each block.

Backward stability ( $f (d) = f (d + δ d)$ )

A numerical algorithm is backward stable when its computed output is the exact output for input data perturbed by an amount comparable to the arithmetic's round-off error. The allowed input perturbations and their norm must be specified. A bound on backward error gives a forward error bound only after accounting for the problem's condition number. For example, a singular value decomposition with orthonormal computed outer factors is backward stable in the spectral norm if

∥ A - U Σ V^{T} ∥_{2} \leq c u ∥ A ∥_{2}

, with machine precision

u

and a modest dimension-dependent constant

c

Root finding

Find a value

x

satisfying

f (x) = 0

for a given function. Methods include the bisection method, which retains a bracket around a root, and derivative-based iteration. Finding a root numerically is distinct from the root as a mathematical object.

Residual root error bound ( $∣ x - x_{*} ∣ \leq ∣ f (x) ∣/ m$ )

Suppose

f

is continuous on an interval, differentiable in its interior, and

∣ f^{'} ∣ \geq m > 0

throughout the interior. If

x_{*}

is a root in the interval and

x

is another point in it, the mean value theorem gives

f (x) - f (x_{*}) = f^{'} (ξ) (x - x_{*})

for an intermediate point

ξ

. Taking absolute values proves the stated bound from the residual of a root equation. This converts an observed residual into a stopping certificate, provided the evaluation error in that residual is also controlled.

Residual of a root equation ( $r = f (x)$ )

The residual of a root equation

f (x) = 0

at a proposed solution

x

is the value

f (x)

. It measures violation of the equation, while the forward error

∣ x - x_{*} ∣

measures distance from a chosen exact root. Rescaling

f

rescales the residual of a root equation without moving its roots, so a small residual needs a scale or a derivative bound before it certifies a small forward error.

Secant method ( $x_{k + 1} = x_{k} - f (x_{k}) (x_{k} - x_{k - 1}) / (f (x_{k}) - f (x_{k - 1}))$ )

The secant method approximates a zero of a function by repeatedly taking the zero of a secant line through the two most recent sampled points. The formula follows by setting the line

f (x_{k}) + (x - x_{k}) (f (x_{k}) - f (x_{k - 1})) / (x_{k} - x_{k - 1})

equal to zero. It requires distinct sampled positions and a nonzero difference of sampled values, unless a root has already been found. After initialization it needs one new function evaluation per step. Unlike bisection, it need not retain a bracket. For a sufficiently smooth function near a simple root, its usual local order of convergence is

(1 + 5) /2

.^[ref]

Fixed-point iteration

Fixed-point iteration seeks a fixed point of a map

F

by repeatedly setting

x_{k + 1} = F (x_{k})

. Convergence requires additional hypotheses; for example, a contraction mapping on a complete metric space has a unique fixed point reached by this iteration.

Affine iteration ( $x_{k + 1} = M x_{k} + c$ )

Iterating an affine map gives

x_{k} = M^{k - 1} x_{1} + \sum_{j = 0}^{k - 2} M^{j} c

. A fixed point exists exactly when

c

is in the image of a linear map

I - M

. If the spectral radius of

M

is below one, its matrix powers tend to zero and every orbit converges to the unique fixed point

(I - M)^{- 1} c

Order of convergence ( $lim_{k \to \infty} ∥ e_{k + 1} ∥/∥ e_{k} ∥^{p} = C \in (0, \infty)$ )

An iteration converging to

x_{*}

has exact asymptotic order

p

if its nonzero errors

e_{k} = x_{k} - x_{*}

satisfy the displayed limit. A bound

∥ e_{k + 1} ∥ \leq C ∥ e_{k} ∥^{p}

instead asserts at least order

p

. Finite termination is treated separately because the error ratios become undefined.

Linear convergence ( $∣ e_{k + 1} ∣/∣ e_{k} ∣ \to ρ \in (0, 1)$ )

An iteration of a map has linear convergence to its limit when its nonzero errors

e_{k} = x_{k} - α

satisfy

∣ e_{k + 1} ∣/∣ e_{k} ∣ \to ρ

for some

0 < ρ < 1

. The error eventually shrinks by approximately the same factor at each step. This is exact order of convergence one, in contrast to quadratic convergence.

Q-linear convergence ( $lim sup_{k} ∥ e_{k + 1} ∥/∥ e_{k} ∥ < 1$ )

A sequence with errors

e_{k} = x_{k} - x^{*}

converges Q-linearly in a specified norm if there is a constant

q < 1

with

∥ e_{k + 1} ∥ \leq q ∥ e_{k} ∥

eventually. The limsup of consecutive nonzero error ratios is its asymptotic factor. This is stronger than merely having a bound

∥ e_{k} ∥ \leq C q^{k}

, which is called R-linear convergence. The norm and whether an iteration is one coordinate step or an entire sweep must be specified when quoting a factor.

Cubic convergence ( $∥ e_{k + 1} ∥ \leq C ∥ e_{k} ∥^{3}$ )

An iteration has at least cubic convergence when its errors eventually satisfy the displayed inequality with fixed

C > 0

. Exact order of convergence three requires a positive finite limit of

∥ e_{k + 1} ∥/∥ e_{k} ∥^{3}

. If

C^{1/2} ∥ e_{N} ∥ < 1

, setting

d_{k} = C^{1/2} ∥ e_{k} ∥

gives

d_{k + 1} \leq d_{k}^{3}

, hence

d_{N + j} \leq d_{N}^{3^{j}}

Fixed-point arithmetic

Rounding

Shooting method

Golden-section search

An optimization algorithm that brackets an extremum of a unimodal function on an interval. Two interior probes divide the bracket in the golden ratio; after discarding one side, one probe can be reused, so each subsequent contraction needs only one new function evaluation.

Significant figures

In numerical analysis, significant figures record how many meaningful digits are retained in a written numerical value. Leading zeros locate the decimal point rather than adding precision:

0.0079

and

79

each have two significant figures. Trailing zeros may require an explicit convention to distinguish measured precision from positional notation.

Numerical method

A numerical method is an algorithm for computing an exact or approximate solution to a mathematical problem through finitely many arithmetic operations and other prescribed steps. For example, Gaussian elimination solves a linear system exactly in exact arithmetic, whereas finite precision can introduce roundoff error.

Velocity Verlet method ( $x_{n + 1} = x_{n} + h v_{n} + \frac{1}{2} h^{2} a_{n}$ )

For an ordinary differential equation

\overset{x}{¨} = f (x, t)

with acceleration independent of velocity, the velocity Verlet method advances position and velocity by

a_{n} = f (x_{n}, t_{n}), x_{n + 1} = x_{n} + h v_{n} + \frac{h ^{2}}{2} a_{n}, v_{n + 1} = v_{n} + \frac{h}{2} (a_{n} + f (x_{n + 1}, t_{n} + h)) .

Here

h

is the time step and

t_{n + 1} = t_{n} + h

. For sufficiently smooth

f

, the global error on a fixed time interval is

O (h^{2})

, so halving

h

reduces the leading global error by a factor of four once the asymptotic regime is reached.

Numerical algorithm

A numerical algorithm implements a numerical method through specified finite computation steps, including arithmetic operations and any required function evaluations.

Floating-point arithmetic ( $\pm (d_{0} . d_{1} \dots d_{p - 1})_{β} β^{e}$ )

Floating-point arithmetic represents a real number by a sign, a fixed number of significant digits, and a scale. It trades exactness for a compact representation spanning many orders of magnitude. The picture shows the essential compromise: with the same precision, large numbers have wider gaps between representable neighbors. Thus accuracy is naturally relative to magnitude, rather than a fixed distance along the number line.

Figure 1.
Eight normalized four-bit binary values in each of three intervals, with gaps doubling from one eighth to one half
. **Details**.

In radix

β

, an integer at least two, a nonzero normalized value has the positional notation

y = \pm (d_{0} + \frac{d _{1}}{β} + \dots + \frac{d _{p - 1}}{β ^{p - 1}}) β^{e}, 1 \leq d_{0} < β, 0 \leq d_{j} < β (j \geq 1),

where the digits and exponent

e

are integers. The parenthesized factor is the significand;

p

counts all its significant digits, including the leading digit. Normalization means

d_{0} \neq = 0

, so the same value has one such representation. The allowed exponents satisfy

e_{m i n} \leq e \leq e_{m a x}

. Finite precision restricts the available digits, while the exponent bounds restrict the range. Zero needs a separate representation. This scales a digit string in positional notation by a power of its base.^[ref]

For positive normalized values in

[β^{e}, β^{e + 1})

, increasing the last digit by one gives the constant neighbor spacing

Δ_{e} = β^{e - p + 1} .

At the next scale, the spacing is multiplied by

β

. In the binary numeral system with

p = 4

, the values in

[1, 2)

are

1, 1.125, 1.25, \dots, 1.875

; those in

[2, 4)

are

2, 2.25, \dots, 3.75

; those in

[4, 8)

are

4, 4.5, \dots, 7.5

. The respective gaps are

1/8

1/4

, and

1/2

. Although absolute resolution deteriorates with magnitude,

Δ_{e} / y

stays between

β^{- p}

and

β^{1 - p}

. This scale invariance is why floating-point arithmetic suits both microscopic and astronomical calculations.

The IEEE 754 binary64 format uses the binary numeral system with

p = 53

, including an implicit leading one for normalized values, and normal exponents

- 1022 \leq e \leq 1023

. Its 64 stored binary digits comprise one sign digit, 11 exponent digits and 52 trailing significand digits. Subnormal numbers extend the range toward zero with a different leading-digit rule. In rounding to nearest, an operation returns the representable value closest to its exact result; a tie goes to the value whose last significant binary digit is even. Machine epsilon describes the gap just above one, whereas unit roundoff bounds a single relative rounding error in the normal range. These format and rounding rules are presented in MIT's floating-point lecture.^[ref]

Even an ordinary terminating decimal can require rounding in the binary numeral system. The rational number

0.1 = 1/10

has the repeating positional notation

0.0001100110011 \dots_{2}

. Every finite binary fraction has a reduced denominator that is a power of two, so the factor five in

10

prevents an exact finite binary representation. In binary64, the exact stored values of three decimal literals and their computed sum are

Python expression 0.1 0.2 0.3 0.1 + 0.2 exact stored rational number 3602879701896397/36028797018963968 3602879701896397/18014398509481984 5404319552844595/18014398509481984 1351079888211149/4503599627370496

The exact sum of the first two stored rational numbers is

10808639105689191/36028797018963968

, halfway between two binary64 values. The even-last-digit tie rule rounds it upward. The printed result is 0.30000000000000004, one step

2^{- 54}

above stored 0.3, so 0.1 + 0.2 == 0.3 is false. Both conversion of the inputs and rounding of the addition matter here.

To reproduce the table on a Python installation using binary64, inspect exact fractions rather than only their printed decimal representations:

import sys
assert (sys.float_info.radix, sys.float_info.mant_dig) == (2, 53)
for value in (0.1, 0.2, 0.3, 0.1 + 0.2):
    print(repr(value), value.as_integer_ratio())
print(0.1 + 0.2 == 0.3)  # False

Python's as_integer_ratio() returns a pair of integers whose quotient equals the stored finite value exactly.^[ref] Exact rational numbers are also useful for checking floating-point arithmetic: a short decimal printout can hide the approximation being made.

Historically, David Goldberg's March 1991 ACM Computing Surveys paper explained representation, rounding and catastrophic cancellation together with the IEEE standard's design choices. Oracle provides an authorized edited reprint. Its enduring lesson is practical: use the spacing to understand roundoff error, and inspect intermediate values when mathematically equivalent algorithms behave differently.^[ref]

Floating-point operation

A floating-point operation is one arithmetic operation on values represented in floating-point arithmetic. Algorithmic counts commonly count addition, subtraction, multiplication and division separately. A fused multiply-add can therefore count as two arithmetic operations even when executed by one hardware instruction. Such counts describe work, while elapsed time also depends on memory traffic, parallelism and the implementation.

Floating-point addition

To add two finite floating-point numbers, compare their exponents, shift the smaller significand to a common exponent, add or subtract the signed significands, normalize the result, and apply floating-point rounding. Extra low-order bits are needed to determine correct rounding when shifting or cancellation discards information. A carry beyond the leading digit increases the exponent; cancellation can decrease it. Unlike exact real addition, this operation need not be associative. Changing the child evaluation order in a fixed binary expression tree need not change its parenthesization, whereas reassociation changes which intermediate results undergo floating-point rounding.

Half-precision floating-point format (binary16)

 2  0

A sixteen-bit floating-point representation with one sign bit, five exponent bits and ten stored fraction bits. For sign

s

, stored exponent

E

and fraction integer

F

, normal finite values are

(- 1)^{s} (1 + F 2^{- 10}) 2^{E - 15}, 1 \leq E \leq 30, 0 \leq F < 1024.

The implicit leading one gives eleven significant binary digits. With

E = 0

, the value is

(- 1)^{s} F 2^{- 24}

, including signed zero at

F = 0

and subnormal numbers otherwise. The stored exponent

E = 31

is reserved for infinities at

F = 0

and NaN otherwise. The positive normal range is

2^{- 14}

through

65504

; the least positive subnormal number is

2^{- 24}

Double-precision floating-point format (binary64)

A floating-point representation occupying

64

bits, or

8

bytes, with one sign bit,

11

exponent bits and

52

explicitly stored fraction bits. The implicit leading bit of a normal number gives

53

bits of significand precision. Near one, adjacent numbers are separated by

2^{- 52}

; rounding to nearest has unit roundoff

2^{- 53} \approx 1.11 \times 1 0^{- 16}

in the normal range. These are relative precision scales, distinct from underflow limits or a universal absolute error bound.^[ref]

Floating-point overflow

Floating-point overflow occurs when the magnitude of a computed result exceeds the largest finite representable floating-point number. Large values of an exponential function are a common source. In a normalized weighted sum, replacing exponentials

e^{r_{i}}

e^{r_{i} - r_{m a x}}

before division prevents this particular overflow while leaving the normalized weights unchanged; this is an application of the log-sum-exp construction.

NaN

A special floating-point arithmetic value representing an undefined or unrepresentable numerical result. A displayed NaN is not a numerical zero and is not a valid p-value.

Floating-point rounding

Replacement of an exact real arithmetic result by a representable floating-point number under a specified rounding rule. Repeated products can depend on evaluation order, so a recursive rewrite intended to reproduce a floating-point loop should preserve its sequence of accumulator operations. An exact algebraic invariant does not by itself establish equality of rounded machine results.

Floating-point underflow

Floating-point underflow occurs when a nonzero numerical result is smaller in magnitude than the representable range or loses precision in the subnormal range. A long product of probabilities can become zero despite a mathematically positive result. Adding logarithms avoids forming the product; adding alternative path probabilities then requires log-sum-exp or controlled rescaling.

Floating-point representation

Significand ( $m$ )

The significant-digit factor

m

in a floating-point representation

(- 1)^{s} m β^{e}

, where

s

is the sign bit,

β

the radix and

e

the exponent. For a normal positive binary significand,

1 \leq m < 2

, and the leading one need not be stored. In a subnormal number, the leading digit is zero at the minimum normal exponent scale. The significand is a numerical factor, distinct from the stored trailing fraction bits that encode part of it.

Floating-point number

Compensated summation

Compensated summation retains a correction for low-order information lost to roundoff error in successive floating-point arithmetic additions. Kahan summation is a common implementation; it can substantially reduce accumulated error when many terms are added.

Kahan summation algorithm

A summation algorithm that keeps a correction for low-order information lost when adding floating-point numbers. Feeding that correction into subsequent additions can greatly reduce accumulated rounding error compared with straightforward summation.

Machine epsilon ( $ϵ_{mach} = β^{1 - p}$ )

Here machine epsilon means the gap from one to the next larger value in normalized floating-point arithmetic:

ϵ_{mach} = next (1) - 1 = β^{1 - p} .

It is a relative resolution scale near one, not the smallest positive representable value and not a fixed absolute error for all magnitudes. In binary64,

p = 53

gives

ϵ_{mach} = 2^{- 52} \approx 2.220446049250313 \times 1 0^{- 16}, u = 2^{- 53} \approx 1.1102230246251565 \times 1 0^{- 16},

where

u

is unit roundoff for rounding to nearest. Some authors instead call

u

their machine epsilon; MIT 18.335J uses this half-gap convention. Always check which definition an error bound uses.^[ref]

The gap below one is smaller: the preceding binary64 value is

1 - 2^{- 53}

, whereas the next value above one is

1 + 2^{- 52}

. Under nearest-even rounding, 1.0 + 2.0**-53 is exactly 1.0: the exact sum is a tie and one has the even final binary digit. By contrast, 1.0 + 2.0**-52 is the next value. Thus searching for the smallest increment that changes one by trial addition also depends on the rounding convention. At

2^{53}

the next larger binary64 value is

2^{53} + 2

, illustrating why a small absolute increment can disappear at large magnitude.

Scale-aware stopping criterion for Newton iteration ( $∣ x_{n + 1} - x_{n} ∣ \leq C ϵ max (1, ∣ x_{n} ∣)$ )

Near a well-conditioned simple root, a numerically stable Newton root-finding iteration cannot usually resolve steps far below the local spacing of floating-point numbers. That spacing has relative scale comparable to machine epsilon

ϵ

, motivating the displayed position-based criterion with a modest constant

C

. It combines an absolute scale near zero with a relative scale for large inputs. A residual-only test is changed by rescaling the equation, whereas the root and exact Newton step are unchanged. A tiny step alone is not a universal accuracy certificate: derivative errors, stagnation and cycles require separate consideration.

Unit roundoff ( $u = \frac{1}{2} β^{1 - p}$ )

For floating-point arithmetic with precision

p

in radix

β

and rounding to nearest, unit roundoff is the half-gap scale

u = β^{1 - p} /2

. If a nonzero exact real number

z

lies in the normal finite range, its nearest rounded value satisfies

fl (z) = z (1 + δ), ∣ δ ∣ \leq u .

Indeed, the absolute distance to the nearest value is at most

Δ_{e} /2

, and

∣ z ∣ \geq β^{e}

gives the bound. For a correctly rounded addition, subtraction, multiplication or division, apply this estimate to the exact result of that one operation on the stored operands. The bound does not include earlier roundoff error or automatically bound a whole algorithm. Results outside the finite range and tiny subnormal numbers require separate treatment. Rounding in a directed mode also changes the bound.^[ref]

Subnormal number ( $k β^{e_{m i n} - p + 1}, 1 \leq k < β^{p - 1}$ )

A subnormal number is a nonzero representable value smaller in magnitude than the least positive normalized value

β^{e_{m i n}}

. Keep the exponent scale fixed at

e_{m i n}

and allow a zero leading digit: the positive values are

k β^{e_{m i n} - p + 1}

for integers

1 \leq k < β^{p - 1}

. This fills the gap toward zero with constant absolute spacing, at the cost of fewer significant digits. Goldberg describes this gradual approach to zero and why it matters for differences of small numbers.^[ref]

For binary64 the least positive normalized value is

2^{- 1022}

, and subnormal numbers are spaced by

2^{- 1074}

. With nearest rounding and gradual underflow enabled, an exact result

z

with

∣ z ∣ < 2^{- 1022}

has absolute rounding error at most

2^{- 1075}

. A uniform relative bound by unit roundoff fails here:

z = 2^{- 1075}

rounds to zero under the even-last-digit tie rule, giving relative error one. Machine epsilon measures resolution near one; it does not measure the distance from zero to the first positive value.

Catastrophic cancellation ( $1 + x - 1 = \frac{x}{1 + x + 1}$ )

Catastrophic cancellation occurs when subtraction of nearly equal approximations removes their matching leading digits and exposes earlier roundoff error. The subtraction can itself be exact on the stored operands: the information may already have been lost while obtaining them. Goldberg distinguishes this from cancellation of accurately known quantities, which can be harmless.^[ref]

For

x > 0

, consider the single real-valued function

f (x) = 1 + x - 1 = \frac{x}{1 + x + 1} .

The rationalization of a square-root difference proves equality in exact arithmetic. The two evaluations use the same input

x

, but the direct version first forms a value near one and then subtracts one. The rearranged version keeps the small input in its numerator and divides by a denominator near two.

Figure 1.
Direct and rationalized evaluation of a small square-root difference against a high-precision reference for the same binary64 inputs
. **Details**.

For binary64 x = 1e-16, the rounded sum 1.0 + x is one, so sqrt(1.0 + x) - 1.0 is zero. The rationalized evaluation returns approximately

5 \times 1 0^{- 17}

. A 100-digit Decimal calculation using the exact stored input gives relative errors of one and about

2.50 \times 1 0^{- 17}

respectively. At x = 1e-12, the direct result is about

5.0004445029 \times 1 0^{- 13}

, with relative error

8.89 \times 1 0^{- 5}

; the rearranged result is about

4.99999999999875 \times 1 0^{- 13}

, with relative error

1.21 \times 1 0^{- 17}

The plot uses Decimal.from_float(x) to preserve the binary input exactly, computes

x / (1 + x + 1)

at 100 decimal digits, and verifies the reference at 160 digits. Exact conversion of a stored binary value is documented by Python.^[ref] Each error is

∣ f - f ∣/∣ f ∣

, with the computed result converted exactly to Decimal before comparison. Hence the curves measure evaluation error for one common input, rather than a mixture of input conversion error and evaluation error. The plotted interval

1 0^{- 16} \leq x \leq 1 0^{- 1}

is far above the subnormal number range.

The original problem is well conditioned with respect to relative changes in

x

. Using the derivative

f^{'} (x) = 1/ (2 1 + x)

and rationalization of a square-root difference, its relative condition number is

κ_{f} (x) = \frac{x f ^{'} ( x )}{f ( x )} = \frac{1 + x + 1}{2 1 + x}, \frac{1}{2} < κ_{f} (x) < 1 (x > 0) .

It tends to one as

x

tends to zero. The direct evaluation instead introduces the intermediate difference

s - 1

with

s = 1 + x

; an error of order unit roundoff in

s

is large relative to

f (x) \approx x /2

. This is a failure of the chosen algorithm's accuracy despite a small condition number for the same real-valued function. The rationalized evaluation avoids that intermediate catastrophic cancellation, retaining relative accuracy near the unit roundoff scale throughout the plotted interval. Recognizing this structure can improve a calculation far more than merely printing extra digits.

Newton root-finding iteration ( $t_{r + 1} = t_{r} - F (t_{r}) / F^{'} (t_{r})$ )

The Newton root-finding iteration seeks a zero of a function

F

by repeatedly solving its tangent line approximation. At a current guess

t_{r}

with derivative

F^{'} (t_{r}) \neq = 0

, set

0 = F (t_{r}) + F^{'} (t_{r}) (t_{r + 1} - t_{r}), t_{r + 1} = t_{r} - \frac{F ( t _{r} )}{F ^{'} ( t _{r} )} .

Geometrically, follow the tangent line to the horizontal axis, then return to the function graph for the next tangent line. Close to a suitable zero of a function, the error becomes dramatically smaller.^[ref]

For a concrete square root calculation, take the polynomial

F (t) = t^{2} - 2

and start at

t_{0} = 1

. The Newton root-finding iteration becomes

t_{r + 1} = (t_{r} + 2/ t_{r}) /2

, giving

r 01234 t_{r} 1 1.5 1.416666666667 1.414215686275 1.414213562375 ∣ t_{r} - 2 ∣ 4.14214 \times 1 0^{- 1} 8.57864 \times 1 0^{- 2} 2.45310 \times 1 0^{- 3} 2.12390 \times 1 0^{- 6} 1.59486 \times 1 0^{- 12}

The errors use the unrounded iterates. After the first step, each iterate overestimates

2

, and the exact identity

t_{r + 1} - 2 = \frac{( t _{r} - 2 ) ^{2}}{2 t _{r}}

exhibits quadratic convergence: once close, the number of correct decimal digits roughly doubles per step. This speed makes the Newton root-finding iteration useful when accurate derivatives are available.

Figure 1.
Newton tangent steps toward a square root and a nonconvergent two-cycle
. **Details**.

Global convergence requires additional care. For the polynomial

F (x) = x^{3} - 2 x + 2

, starting at

x_{0} = 0

gives

F (0) = 2

F^{'} (0) = - 2

, hence

x_{1} = 1

. Then

F (1) = 1

F^{'} (1) = 1

gives

x_{2} = 0

. The sequence repeats

0, 1, 0, 1, \dots

forever, although a real polynomial root lies in

(- 2, - 1)

by the intermediate value theorem. Neither point in the cycle is a zero of a function, and neither has zero derivative. A maintained sign-changing interval and fallback bisection method can protect a root solver from unsuitable Newton steps.

Broyden's method ( $B_{k + 1} = B_{k} + (y_{k} - B_{k} s_{k}) s_{k}^{T} / (s_{k}^{T} s_{k})$ )