Introduction
Implicit Runge-Kutta methods have a number of desirable properties.
The Gauss-Legendre methods for example are self-adjoint meaning that they provide the same solution when integrating forwards or backwards in time.
This loads packages defining some example problems and utility functions.
Coefficients
A generic framework for implicit Runge Kutta methods has been implemented.
The focus so far is on methods with interesting geometric properties and currently covers the following schemes:
ImplicitRungeKuttaGaussCoefficients
ImplicitRungeKuttaLobattoIIIACoefficients
ImplicitRungeKuttaLobattoIIIBCoefficients
ImplicitRungeKuttaLobattoIIICCoefficients
ImplicitRungeKuttaRadauIACoefficients
ImplicitRungeKuttaRadauIIACoefficients
The derivation of the method coefficients can be carried out to arbitrary order and arbitrary precision.
Coefficient generation
Start with the definition of the polynomial defining the abscissae of the 
stage coefficients. For example the abscissae for Gauss-Legendre methods are defined as 
.
Univariate polynomial factorization gives the underlying irreducible polynomials defining the roots of the polynomials.
Root objects are constructed to represent the solutions (using unique root isolation and Jenkins-Traub for the numerical approximation).
Root objects are then approximated numerically for precision coefficients.
Condition estimates for Vandermonde systems governing the coefficients yield the precision to take in approximating the roots numerically.
Specialized solvers for non-confluent Vandermonde systems are then used to solve equations for the coefficients (see [GVL96]).
One step of iterative refinement is used to polish the approximate solutions and to check that the coefficients are obtained to the requested precision.
This generates the coefficients for the two stage fourth order Gauss-Legendre method to 50 decimal digits of precision.
The coefficients have the form {a,bT
,cT
}.
This generates the coefficients for the two stage fourth order Gauss-Legendre method exactly. For high order methods generating the coefficients exactly can often take a very long time.
This generates the coefficients for the six stage tenth order RaduaIA implicit Runge-Kutta method to 20 decimal digits of precision.
Examples
Load an example problem.
This problem has two invariants which should remain constant. A numerical method may not be able to conserve these invariants.
This solves the system using an implicit Runge-Kutta Gauss method. The order of the scheme is selected automatically depending on the values of the AccuracyGoal and PrecisionGoal options of NDSolve.
A plot of the error in the invariants shows an increase as the integration proceeds.
The method ImplicitRungeKutta has options AccuracyGoal and PrecisionGoal that specify the absolute and relative error to aim for in solving the nonlinear system of equations.
These options have the same default values as the corresponding options in NDSolve, since often there is little point in solving the nonlinear system to much higher accuracy than the local error of the method.
The first invariant is the Hamiltonian of the system and the error is now bounded, as it should be since the Gauss implicit Runge-Kutta method is a symplectic integrator.
The second invariant is conserved exactly (up to roundoff) since the Gauss implicit Runge-Kutta method conserves quadratic invariants.
This defines the implicit midpoint method as the one stage implicit Runge-Kutta method of order two.
For this problem it can be more efficient to use a fixed point iteration instead of a Newton iteration to solve the nonlinear system.
At present, the implicit Runge-Kutta method framework only uses fixed step sizes and does not use banded Newton techniques for uncoupling the nonlinear system.
