Implicit Runge-Kutta method

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

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^{\prime \prime }{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.

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 $hL{\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.

The implicit midpoint rule advances an autonomous ordinary differential equation byIt is a second-order, A-stable, time-symmetric numerical method.

The implicit midpoint rule is the midpoint method with $y_{n+1}=y_n+hf(t_n+h/2,(y_n+y_{n+1})/2)$. It is a one-stage implicit Runge-Kutta method.

For a uniformly Lipschitz continuous state field $f$ with constant $L$, an implicit midpoint rule step is uniquely solvable when $hL/2<1$, by the contraction mapping theorem. Write its step map as ${\Psi }_h(u)=u+h{\Phi }_h(u)$, where ${\Phi }_h(u)=f(t+h/2,u+h{\Phi }_h(u)/2)$. Subtracting the increment equations givesTherefore $\mathrm{Lip}({\Phi }_h)\le L/(1-hL/2)$, and the triangle inequality gives the bound for ${\Psi }_h$. Combining it with the implicit midpoint consistency bound yieldsfor $L>0$. This follows by summing the geometric error recurrence with forcing $K{h}^3/(1-hL/2)$. For $hL\le 1$ and $nh\le T$, $R_h^n\le e^{2LT}$, proving second-order global error with exact starting data. If $L=0$, the direct bound is $\Vert e_n\Vert \le \Vert e_0\Vert +KT{h}^2$.

For an exact solution $y^{\prime }=f(t,y)$, define the normalized local truncation error of the implicit midpoint rule byIf $\Vert y^{\prime \prime }\Vert \le M_2$, $\Vert y^{\prime \prime \prime }\Vert \le 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 giveFor the first bound, integrate the second-order remainder of $y^{\prime }(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^{\prime }(t)=f(t,y(t))$ finish the proof by the triangle inequality. The unnormalized step residual is $h{\tau }_n=O(h^3)$.

The Crank--Nicolson method applies the trapezoidal rule to a semidiscrete evolution equation. For a linear system $u^{\prime }=Lu$ its update is

For symmetric negative semidefinite matrices $A,B$, each factor $C_A(h)=(I-hA/2{)}^{-1}(I+hA/2)$ has norm at most one, since its eigenvalues are $(1+h\lambda /2)/(1-h\lambda /2)$ with $\lambda \le 0$. The same holds for $B$, so their ordered product is a contraction for every $h\ge 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$.

Given distinct nodes $c_i$, a collocation Runge-Kutta method uses the Lagrange basis ${\ell }_j$ with $a_{ij}={\int }_0^{c_i}{\ell }_j(\tau )d\tau$ and $b_j={\int }_0^1{\ell }_j(\tau )d\tau$. 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.

For distinct nodes $c_1,c_2\in [0,1]$, the two-stage collocation Runge-Kutta method has stability functionThe 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\ge r$, equivalently $c_1+c_2\ge 1$, is necessary and sufficient for A-stability, using the maximum modulus principle and the bound at infinity.

For $s$ distinct collocation nodes, let ${\pi }_s(t)={\prod }_{i=1}^s(t-c_i)$. If ${\int }_0^1{t}^j{\pi }_s(t)dt=0$ for $0\le 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\le k\le 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 $2s$, and Radau IIA methods reach $2s-1$.

For $s$ collocation nodes, the residual of the collocation polynomial vanishes at each node and contains the factor $\omega (\tau )={\prod }_i(\tau -c_i)$. Quadrature exactness through degree $q-1$ implies ${\int }_0^1\omega (\tau )p(\tau )d\tau =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=2s$, while the first nonexact quadrature moment bounds the order from above.

For a collocation Runge-Kutta method, $a_{ij}={\int }_0^{c_i}{\ell }_j(t)dt$, where the Lagrange interpolation basis obeys ${\sum }_j{\ell }_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.

A Lobatto IIIA method collocates at the Lobatto quadrature nodes, including both endpoints of the time interval. Its three-stage version hasThe Butcher order conditions give order four. Its stability function iswhich is A-stable but not L-stable. It is not algebraically stable: the first diagonal entry of the defining matrix is $2b_1{a}_{11}-b_1^2=-1/36$. This distinguishes linear A-stability from nonlinear algebraic stability.

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 $2s$.

Choose $c_1,\dots ,c_s$ as the polynomial roots of $P_s(2c-1)$ in $(0,1)$, and form the Lagrange interpolation polynomials ${\ell }_j(t)={\prod }_{m\ne j}(t-c_m)/(c_j-c_m)$. The Gauss collocation method has $a_{ij}={\int }_0^{c_i}{\ell }_j(t)dt$ and $b_j={\int }_0^1{\ell }_j(t)dt$. This gives explicit coefficient formulas after finding polynomial roots: if ${\ell }_j(t)={\sum }_r{d}_{jr}{t}^r$, these integrals are ${\sum }_r{d}_{jr}{c}_i^{r+1}/(r+1)$ and ${\sum }_r{d}_{jr}/(r+1)$. Gaussian quadrature gives order $2s$.

The two-stage Gauss--Legendre Runge-Kutta method uses nodes $c_{1,2}=1/2\mp \sqrt{3}/6$, weights $b_1=b_2=1/2$ and coefficients $a_{11}=a_{22}=1/4$, $a_{12}=1/4-\sqrt{3}/6$, $a_{21}=1/4+\sqrt{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.

An $s$-stage Radau IIA method collocates at the right-endpoint Radau nodes. It has order $2s-1$ and is both algebraically stable and L-stable.

Suppose an $s$-stage Runge-Kutta method has order $2s-1$, an invertible matrix $A$, and $b^T{A}^{-1}1=1$. Its stability function is $R(z)=1+z{b}^T(I-zA{)}^{-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^{2s})$, identifying the $[s-1/s]$ Padé approximant. Uniqueness follows because the cross-multiplied difference of two candidates has degree at most $2s-1$ but a zero of order at least $2s$, 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.

A Runge-Kutta method is algebraically stable when $b_i\ge 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.

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 $\Vert y_{n+1}{\Vert }^2-\Vert y_n{\Vert }^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.

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.

For differences $D_i$ between corresponding stages and differences $F_i$ between their vector fields, a Runge--Kutta step satisfiesTo 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\ge 0$, and the second sum is nonnegative when the matrix $(m_{ij})$ is a positive semidefinite matrix. Thus algebraic stability implies B-stability.

As a quadrature rule, the trapezoidal rule replaces an integrand by the linear interpolation joining its two endpoint values: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)=\alpha x+\beta$ agrees with its linear interpolation, so the trapezoidal rule integrates it exactly. On $[0,2]$, $f(x)=2x+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$ by $1/6$. Its degree of precision of a quadrature rule is exactly one.

If $f\in C^2[a,b]$, its signed approximation error satisfiesA 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 identityobtained 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.

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)|\le \frac{1}{2}M(x-a)(b-x)$, where $M=\max |f^{\prime \prime }|$. Integrating this inequality and using the triangle inequality for an integral proves the bound because ${\int }_a^b(x-a)(b-x)dx=(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$.

On an equal-panel partition of an interval, with $h=(b-a)/n$, sum the trapezoidal rule over adjacent panels: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^{\prime \prime }(x)|$, summing panel quadrature errors givesfor some $\xi \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\le K(b-a)H^2/12$.

For the polynomial $f(x)=x^2$ on $[0,1]$, direct summation yieldsBoth 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 $\delta =0.002$ to the composite trapezoidal rule. That mathematical model converges to $1/3+\delta$, showing why smaller panels alone do not ensure accuracy for noisy samples. It illustrates noise sensitivity of a quadrature rule, rather than measured data.

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.

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.

If $f$ is twice continuously differentiable on $[a,b]$, take a uniform bound $|f^{\prime \prime }|\le M$. On one panel $[u,v]$, two integrations by parts give the signed quadrature errorThe 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.

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=\sqrt{2}$, $T=0$ and implicit differentiation gives $T^{\prime }=-2$, so nearby states expand rather than contract. This is a concrete nonlinear obstruction, rather than merely failure of the sufficient algebraic stability criterion.