Runge-Kutta methods

A Runge-Kutta method can be described by the following equations:

y n+1 =yn +hS i=1 q bi ki ,

ki =f(tn +ci h,yn +hS j=1 q a ij kj ), 1<= i<= q

where ci =S i=1 q a ij ,i=1,...,q. The values ci and bi are the abscissas, respectively the weights of the method. The Butcher matrix associated to the Runge-Kutta method is the following table
c1 a11 ... a1q
... ... ...
cq aq1 ... aqq
b1 ... bq
or, more simple,
c A
b

Another form to write the Runge-Kutta method is the following:

y n+1 =yn+hS i=1 qbif(tn+cih,Yi ), Yi=yn+hS j=1 qa ij f(tn+cjh,Yj ), i=1, ... ,q

If a ij =0, j>= i,i=1,...,q the Runge-Kutta method is explicit (ERK), otherwise is implicit (IRK). If a ij =0, j>= i, it is a diagonal implicit schema (DIRK).

The barrier order for a q-stage Runge-Kutta methods is 2q. The ( Butcher 's) simplified order condition are the following:

S j=1 q a ij cj k-1 = cik/ k , i=1, ... ,q, k=1, ... ,p, denoted by A(p)

S i=1 q bici k-1 = 1/ k , k=1, ... ,p, denoted by B(p)

Pq(2ci-1)=0, i=1,...,q, denoted by C(q)

S i=1 q bia ij ci k-1 =bj(1-cjk), j=1, ... ,q, k=1, ... ,p, denoted by D(p)

S i,j=1 pbici k-1 a ij cj l-1 = 1/ l(k+l) , k=1, ... ,m, l=1, ... ,p, denoted by E(m,p)

where Pq is the Legendre polynomial. If A(m)+D(x)+B(p), p<= x+m+1, p<= 2m+2, the method has at least p-order; A(q)+D(q)+B(q)+C(q) iff the method order is 2q. The methods of 2q order are the Gauss methods.

The stage order of a Runge-Kutta method is the value s=min(p1, ... ,p q+1 ),, where pi is the i-stage order, i.e. the smallest integer for which there is a value di such that ||ei||<= dih pi+1 , where

ei=y(tn-cih)-y(tn)-hS j=1 qa ij f(tn+cjh), i=1, ... ,q,

e q+1 =y(t n+1 )-y(tn)-hS j=1 qbjf(tn+cjh).

The stability function R is defined by R(z)=1+zbT(I-zA) -1 e, where I is the q x q identity matrix and e=(1, ... ,1)T. The Runge-Kutta method has the order p iff ez-R(z)=Cz p+1 + O (z p+2 ). The absolute stability region can be computed using the identity A={z in C| |R(z)|<= 1 }. The Runge-Kutta method is A-stable if the negative half complex plane is included in A. It is L-stable if it is A-stable and lim Re (z) to -infinity |R(z)|=0. The Daniel-Moore conjuncture for Runge-Kutta methods is similar with the Dahlquist ' second barrier: the order of a A-stable Runge-Kutta method is resctricted by the inequality p<= 2q.

The algebraic stability is a more strong restriction as the A-stability for Runge-Kutta methods. A Runge-Kutta methods is algebraic stable if bi>= 0, i=1, ... ,q and M=BA+ATB-bbT, where B=diag(bi), is a non-negative defined matrix. The algebraic-stable Runge-Kutta methods are superior to the A-stable Runge-Kutta methods in the case of a nonlinear stiff problem.

Let C=(cij/j) i,j , V=(ci j-1 ) i,j , N=(1/i) i,j , D=diag(bi), F=b1(e,0, ... ,0).

The highest order 2q of a Runge-Kutta method is attained by the Gauss methods, for which A=CV -1 . The 2q-1 order is attained by Radau methods -- there are two classes of methods, for which the abscissas are the roots of I: d q-1 / dx q-1 (xq(1-x) q-1 )=0, or II: d q-1 / dx q-1 (x q-1 (1-x)q)=0. Moreover, for Radau IA method, A=D -1 (VT) -1 (N-C)TD, respectively, for Radau IIA method, A=CV -1 . The abscissas of the 2q-2 order Lobatto methods are the roots of d q-1 / dx q-1 (xq(1-x)q)=0. Moreover, for Lobatto IIIA methods, A=CV -1 , for Lobatto IIIB methods, A=D -1 (VT) -1 (N-C)TD, and for Lobatto IIIC methods, A=(C-F)V -1 . The Gauss method and the type A methods are A-stable, and those of type IA, IIA and IIIC are L-stable. The Gauss , Radau IA, Radau IIA and Lobatto IIIC methods are algebraic-stable.

In a singly-implicit Runge-Kutta method (SIRK) the matrix A has one-point real spectrum. In multi-implicit Runge-Kutta methods (MIRK) the matrix A has all eigenvalues real and distinct so that A matrix is similar to a diagonal matrix with the eigenvalues on the diagonal. An A-stable q-stage MIRK method has a most q+1 order, and an algebraic-stable SIRK of the order q+1 cannot overpass the barrier q<= 4. Also, an algebraic-stable DIRK method cannot exceed order 4.