Efficient Implementation Of Fully Implicit Runge-Kutta Methods And Its Application To Time-discretization PDEs

S-stage Runge-Kutta methods for solving stiff initial value problems of ordinary differential equations a3 b~f(t~i-c~h.g3). (1.2b) 31 S where h > Ois the stepsize. a23. b3. c3 are real coefficient satisfying >3 ~ = j=1 1.2 s). y~, is the approximation of the true solution y(t) of problem (1.1) at the point t~ a -~- nh.. g~ y(t~ c~h)(i ?1.2. ...s). (1.2~ can be written in a more compact form { G= ey~ h(A ?Im)F(G). (1 .3a) Yn-~- 1 Yn ?h(bT ?Im)F(G). (1.3b) where A Lag,: ~ RS~s.b = ~b1.b9. 5iT e ~1.i... .11 R5 , G = T T~晻 . , F(G~ = ~ . f(t~梒5h.g5)7l ~91 9~ I,~, denotes m x rn identity matrix, the symbol ~ denotes Kronecker product. We call the method fully implicit Runge-Kutta method only when the coefficient ma- trix A is full. In this paper.we are devoted to build up a uniform approach for efficiently implementing all different type of fully implicit Runge-Kutta methods. Formlv Hairer et al only proposed the discrete scheme for implementing Radau hA Runge-Kutta methods. Because there are many kinds of fully implicit Runge-Kutta methods. the first step of our work is to compute the coefficients of the method of each different type and different order. In the paper. we give the computational schemes in detail for Gauss, Radau hA. Radau IA and Lobatto IIIC methods. which are widely discussed in papers. 8 Secondly, we diagonalize the matrix A of s-stage Runge-Kutta method through analogical transformation: AT. (1.4) where A =diag()q.A2,...,As). Th1 = [x12~,x5~. (1.5) )~ ~ , )~ and ~i. x2,?, r5 are the eignvalues and eignvectors of matrix A re- spectively. The purpose of the diagonalization is the preparation for Butcher抯 transformation. It can be proved that for the abovementioned Runge-Kutta meth- ods. the matrix A is destined to have s different eignvalues, and at least have a complex conjugate eignvalue pair, so the diagonalization is workable. Usually, we apply implicit double QR decomposition and opposite product method to compute all the eignvalues and eignvectors of matrix A. Thirdly we should give the start values and all kinds of needed parameters for computation. For this reason, we let t~ = a. y~ = ~. n = 0. and appropriately set the start stepsize h = h0, local error tolerance tol, and decide the iteration tolerance toli accordingly. Moreover, we need to give some additional informa- tion in order to make our code can freely switch among the functions listed below according to the character of th? solved problem and the request of the users: 1. Solving ordinary differential equations. differential algebra equations or initial boundary value problems of Partial Differential equations through the method of lines: 2. Fix-stepsize/fix-order computation. variable-stepsize/ftx-order computa- tion or variable-stepsize/varia...