INTEVEP, S.A. Centro de Investigaci´ on y Apoyo Tecnol´ ogico. Filial de PDVSA. Los Teques, Edo. Miranda. VENEZUELA.
METODO RUNGE-KUTTA IMPLICITO DE SEXTO ORDEN TIPO LOBATTO PARA RESOLVER SISTEMAS DE ECUACIONES DIFERENCIALES ORDINARIAS CON CONTROL DEL PASO ANDRES L. GRANADOS M. (EPPR/322).
Departamento de Producci´ on (EPPR). Secci´ on Manejo de Crudos (EPPR/3).
PROYECTO No. 5707. REPORTE TECNICO No. INT-EPPR/3-NT-92-003.
MARZO, 1992.
PANAMERICAN WORKSHOP FOR APPLIED AND COMPUTATIONAL MATHEMATICS. Sim´on Bol´ıvar University. January 10-15. 1993. Caracas, Venezuela.
LOBATTO IMPLICIT SIXTH ORDER RUNGE-KUTTA METHOD FOR SOLVING ORDINARY DIFFERENTIAL EQUATIONS WITH STEPSIZE CONTROL ANDRES L. GRANADOS M. INTEVEP, S.A. Production Department (EPPR/32). Venezuela Government Research Facility for Petroleum Industry. Los Teques, Estado Miranda. Venezuela.
ABSTRACT.-
A method for solving ordinary differential equations has been developed using implicit Runge-Kutta methods. The implicit Runge-Kutta methods used are based in two quadratures of Lobatto type. The first quadrature produces the principal Runge-kutta method which is of sixth order, while the second quadrature produces a Runge-Kutta method of third order which is embedded in the former. Both implicit RungeKutta methods constitute the Lobatto embedding form of third and sixth orders with four stages. The most important advantage of this method is that it is implicit only in the second and third stages, which reduces considerably the costs of computer calculations. The Butcher notation is used here for the analysis of the studied methods. In order to solve, for each step, the system of non-linear equations in the implicit auxiliary variables k2 and k3 , an explicit Runge-Kutta method of four stages and fourth order is defined for the same intermediate points, such as the implicit sixth order Runge-Kutta method. This explicit method estimates the inicial values for the aforementioned auxiliary variables, and then, an iterative method of the type “fixed point” is used to solve the system of non-linear equations for each step. With the third order Runge-Kutta method, an estimation of the local truncation error may be calculated using a comparison with the sixth order method. This aspect is used to control the step size when tolerances for the relative and absolute global errors are specified. An algorithm is presented to do this step control automatically. The implicit method, as is exposed here, is really useful and has demonstrated to be efficient to solve huge and stiff systems of ordinary differential equations. Finally, convergence criteria and stability analysis are studied for the implicit Runge-Kutta methods presented here.
JORNADAS PANAMERICANAS DE MATEMATICAS APLICADAS Y COMPUTACIONALES. Universidad Sim´on Bol´ıvar. 10-15 de Enero. 1993. Caracas, Venezuela.
METODO RUNGE-KUTTA IMPLICITO DE SEXTO ORDEN TIPO LOBATTO PARA RESOLVER SISTEMAS DE ECUACIONES DIFERENCIALES ORDINARIAS CON CONTROL DE PASO ANDRES L. GRANADOS M. INTEVEP, S.A. Departamento de Producci´ on (EPPR/32). Centro de Investigaci´ on y Apoyo Tecnol´ ogico. Filial de PDVSA. Los Teques, Estado Miranda. Venezuela.
RESUMEN.-
En este trabajo se ha desarrollado un m´etodo para resolver ecuaciones diferenciales ordinarias usando m´etodos Runge-Kutta impl´ıcitos. Los m´etodos Runge-Kutta aqu´ı empleados est´ an basados en dos cuadraturas del tipo Lobatto. La primera cuadratura conforma el m´etodo Runge-Kutta principal que es de sexto orden; la segunda cuadratura conforma un m´etodo de Runge-Kutta de tercer orden y est´a acoplado con el primero. Ambos m´etodos de Runge-Kutta impl´ıcitos constituyen la forma acoplada de Lobatto de tercer y sexto orden con cuatro etapas. La ventaja m´ as importante de este m´etodo es que es impl´ıcito s´olo en la segunda y tercera etapa, lo que reduce considerablemente los costos de c´alculo en la computadora, dentro de sus pocas etapas. La notaci´on de Butcher es usada para el an´ alisis de los m´etodos tratados aqu´ı. Con la finalidad de resolver en cada paso el sistema de ecuaciones no lineales originado en las variables k2 y k3 , un m´etodo de Runge-Kutta expl´ıcito de cuatro etapas y cuarto orden es definido para los mismo puntos intermedios, como en el m´etodo de Runge-Kutta impl´ıcito de sexto orden. Este m´etodo expl´ıcito estima los valores iniciales de las mencionadas variables auxiliares, y, luego, un m´etodo iterativo del tipo “punto fijo” es usado para resolver el sistema de ecuaciones no lineales en cada paso. Con el m´etodo de Runge-Kutta de tercer orden puede ser calculado una estimaci´on del error local de truncamiento, compar´ andolo con el m´etodo de sexto orden. Esto es usado para controlar el tama˜ no del paso cuando las tolerancias para los errores relativo y global absoluto son especificados. Se presenta un algoritmo para realizar este control de paso autom´ aticamente. El m´etodo impl´ıcito, tal como se expone aqu´ı, es realmente u ´ til y ha demostrado ser eficiente para resolver sistemas de ecuaciones diferenciales ordinarias r´ıgidas y de gran tama˜ no. Finalmente, son elaborados los criterios de convergencia y los an´ alisis de estabilidad para los m´etodos Runge-Kutta impl´ıcitos aqu´ı presentados.
INTRODUCTION
The principal aim of this work is the development of an algorithm for solving ordinary differential equations based on known implicit Runge-Kutta methods but where the selection of the methods and the application of iterative procedures have been accurately studied in order to produce the least number of numerical calculations and the highest order of exactitude. As it is well known, the implicit Runge-Kutta methods are more stable in the solution of system of ordinary differential equations than the explicit ones. However, the solution of the system of non-linear equations with the auxiliary variables “k”, originated by the implicit dependence, generates an iterative process that can be extremely slow and that can produce a restriction in the size of the step in the method. In order to improve the speed of the iterative process, a good initial value for each auxiliary variable k is estimated by the algorithm using a fourth order, four stages explicit Runge-Kutta method. This method is new and its best characteristic is that the auxiliary variables are evaluated in the same intermediate points as in the principal implicit Runge-Kutta method of sixth order and four stages. The restriction of the size of the step can be overcome by the selection of a special implicit RungeKutta method based on the Lobatto quadrature and that is implicit only in the second and the third stages. This characteristic makes the iterative process less restrictive respective to the size of the step, but, at the same time, reduces considerably the number of numerical calculations. This method conforms the principal implicit Runge-Kutta method. Additionally to the aforementioned two aspects, a third order, three stages implicit Runge-Kutta has been selected to control the step size in each jump of the method. This method is also based on the Lobatto quadrature and its best characteristic is that it is embedded in the principal implicit Runge-Kutta method. This makes the algorithm more efficient for the stepsize control purpose, and, as a consecuence, calculations necessary for solving the ordinary differential equations are the same to those necessary to control the size of the step. Finally, convergence criteria for the iterative process and stability analysis are made for all the Runge-Kutta methods presented here. Additional important information applied for more general implicit Runge-Kutta methods can be found in the appendixes.
3
IMPLICIT RUNGE-KUTTA METHODS
DIFFERENTIAL EQUATIONS
As it is well known, every system of ordinary differential equations of any order may be transformed, with a convenient change of variables, into a system of first order ordinary differential equations [1][2]. This is the reason why only this last type of differential equations will be treated. Let the following system of M first order ordinary differential equations be expressed as
dy i = f i (x, y) dx
i = 1, 2, 3, . . . , M
(1)
being y a M -dimensional function with each component depending on x. This is
dy = f (x, y) dx
(2)
y = y(x) = (y 1 (x), y 2 (x), y 3 (x), . . . , y M (x))
(3)
where
When each function f i (x, y) depends only on each y i the system is said to be uncoupled, otherwise it is said to be coupled. If the system of ordinary differential equations is uncoupled then every differential equation can be solved separately. When the conditions of the solution y(x) are known at a unique specific point, for example y i (xo ) = yoi
x = xo
(4)
both expressions (1) and (4) are said to be an “Initial Value Problem”, otherwise it is a “Boundary Value Problem”. In deed, the system of differential equations (1) is a particular case of a general one stated in the next form [2][3] dy = f (y) ≡ dx
dy i /dx = 1 dy i /dx = f i (y)
but with the additional condition yo1 = xo in (4).
4
if i = 1 if i = 2, 3, . . . , M + 1
(5)
RUNGE-KUTTA METHODS
Trying to make a general formulation, a Runge-Kutta method of order P and equiped with N stages is defined [3] with the expression i yn+1 = yni + cr kri
(6.a)
where the auxiliar M -dimensional variables kr are calculated by kri = h f i (xn + br h, yn + ars ks )
(6.b)
for i = 1, 2, 3, . . . , M
r, s = 1, 2, 3, . . . , N
(6.c)
Notice that index notation has been used. Thus every time an index appears twice or more in a term, this should be summed up to complete its range (In this context, it is not important the number of factors with the same index within each term). A Runge-Kutta method (6) has order P , if for a sufficiently smooth problem (2) and (4) y(xn + h) − yn+1 ≤ Φ(ζ) hP +1 = O(hP +1 )
ζ ∈ [xn , xn + h],
(7)
i.e., if the Taylor series for the exact solution y(xn + h) and for yn+1 coincide up to (and include) the term with hP [6]. The Runge-Kutta method thus defined can be applied for solving initial value problems, and it is used recurrently. Given a point (xn , yn ), it can be obtained the next point (xn+1 , yn+1 ) using the expressions (6), being xn+1 = xn + h
(8)
and h the step of the method. Every time that this is made, the method goes forward (or backward if h is negative) an integration step h in x, offering the solution in consecutive points, one for each jump.
In this way, if the method begins with the conditions (x0 , y0 ) stated by (4), it can calculate
(x1 , y1 ), (x2 , y2 ), (x3 , y3 ), . . . , (xn , yn ), and continue this way, up to the desired boundary in x. In each integration or jump the method reinitiates with the information from the adjacent point inmediately preceding, that is why the Runge-Kutta is within the group of one step methods. It should be notice, however, that the auxiliary variables kri are calculated for all r up to N in each step. These calculations are no more than evaluations of the functions f i (x, y) for intermediate points x + br h in the interval [xn , xn+1 ] (0 ≤ br ≤ 1), multiplied by h. The evaluation of each M -dimensional auxiliary variable kr , represents one stage of the method. Let it now be introduced a condensed representation of the generalized Runge-Kutta method, formerly developed by Butcher [4] and that is presented systematically in the book of Lapidus and Seinfeld [5] and in the book of Hairer, Norsett and Wanner [6]. This last book has a numerous collection of Runge-Kutta 5
methods using the Butcher’s notation and an extensive bibliography. After the paper of Butcher [4] it became customary to symbolize the general Runge-Kutta method (6) by a tableau. In order to illustrate this representation, consider the expresions (6) applied to a method of four stages (N = 4). Accomodating the coefficients ars , br and cr in a adequate form, it is obtained b1 b2 b3 b4
a11 a21 a31 a41
a12 a22 a32 a42
a13 a23 a33 a43
a14 a24 a34 a44
c1
c2
c3
c4
(9)
The aforementioned representation allows for the basic distinction of the following types of RungeKutta methods, according to the characteristics of the matrix ars : If ars = 0 for s ≥ r, then the matrix ars is lower triangular, excluding the principal diagonal, and the method is said to be completely explicit. If, ars = 0 for s > r, then the matrix ars is lower triangular, including the principal diagonal, and the method is said to be semi-implicit. If, conversely, none of the previous conditions are satisfied, the method is said to be implicit. In the cases of implicit methods, it can be noticed that an auxiliary variable kr may be dependant on itself and on any other variables not calculated before in the same stage. That is why the methods are named implicit in these cases. Additionally, the condense representation above described, permits to verify very easily certain properties that the coefficients ars , br , and cr must fulfill. They are 0 ≤ br ≤ 1
(10.a)
ars δs = br
(10.b)
cr δ r = 1
(10.c)
The vector δ is unitary in every component, i.e., δr = 1 ∀r = 1, 2, 3, . . . , N . Those mentioned properties may be interpreted in the following manner: The property (10.a) expresses that the Runge-Kutta is a one step method, and the functions f i (x, y(x)) in (6.b) should be evaluated for x ∈ [xn , xn+1 ]. The property (10.b) results in applying the Runge-Kutta method (6) to the system of differential equations (5), where ks1 = 1 ∀s = 1, 2, 3, . . . , N , and thus the sum of ars in each line r offers the value of br . The property (10.c) means i is obtained from the value of yni , projecting with h an average that in the expression (6.a) the value of yn+1
of the derivatives dy i /dx = f i (x, y), calculated using weighted coefficients, and obviously the sum of all cr should be equal to unit. The coefficients ars , br and cr are determined applying the properties (10) and using some relations that are deduced in the following form: 6
Let the system of ordinary differential equations be expressed according to (5) as an initial value problem of the type dy = f (y) dx x = x0
(5 ) (4 )
y(x0 ) = y0
The Runge-Kutta method applied to this problem is formulated as yn+1 = yn + cr kr
(6.a )
kr = h f (yn + ars ks )
(6.b )
where the auxiliary variable kr is defined as
If now the expansion in Taylor series is made to the component kri of (6.b ), around the point (xn , yn ), being yn = y(xn ), it results that
kri =h f i [δr ] + h fji [ars ksj ] +
h i f [ars ksj ] [art ktk ] 2 jk
+
h i f [ars ksj ] [art ktk ] [aru kul ] 6 jkl
+
h i f [ars ksj ] [art ktk ] [aru kul ] [arv kvm ] + O(h6 ) 24 jklm
(11.a)
where the index rules and the following notation has been used
f i = f i (xn )
fji =
∂f i ∂y j yn
i fjk =
∂ 2 f i ∂y j ∂y k yn
···
(11.b)
Here the fuctions are supposed to be C ∞ (Analytical Functions), and therefore the subindexes in (11.b) are permutable. The variable ksj in the second term of the right side of (11.a) can again be expanded in Taylor series as ksj =h f j [δs ] + h fkj [asα kαk ] +
h j f [asα kαk ] [asβ kβl ] 2 kl (11.c)
h j + fklm [asα kαk ] [asβ kβl ] [asγ kγm ] + O(h5 ) 6 7
In the same way kαk can be expanded as
kαk = h f k [δα ] + h flk [aαδ kδl ] +
h k f [aαδ kδl ] [aα km ] + O(h4 ) 2 lm
(11.d)
and thus sucessively l kδl = h f l [δδ ] + h fm [aδϕ kϕm ] + O(h3 )
(11.e)
kϕm = h f m [δϕ ] + O(h2 )
(11.f )
up to
If a recurrent backward substitution is finally made, as Appendix A decribes in details, it is obtained that
1 i j k 2 kri =h f i [δr ] + h2 [fji f j br ] + h3 [fji fkj f k ars bs + fjk f f br ] 2
1 j k l + h4 [fji fkj flk f l ars ast bt + fji fkl f f ars b2s 2 1 i j k l 3 i + fjk flj f k f l br ars bs + fjkl f f f br ] 6
1 l m k + h5 [fji fkj flk fm f ars ast atu bu + fji fkj flm f l f m ars ast b2t 2 1 j k l m j + fji fkl fm f f ars bs ast bt + fji fklm f k f l f m ars b3s 6 1 i j k l m i l m + fjk flj f k fm f br ars ast bt + fjk flm f f f br ars b2s 2 1 i j k l m 1 i j k l m 2 + fjk fl fm f f ars bs art bt + fjkl fm f f f br ars bs 2 2 +
1 i f f j f k f l f m b4r ] + O(h6 ) 24 jklm 8
(11.g)
Inserting this last expression of the components of kr into the equation (6.a ), and comparing with the next expansion in Taylor series of yn+1 , around yn
i yn+1 =yni + h f i +
h2 i j h3 i (fj f ) + (fji fkj f k + fjk f j f k) 2 6
+
h4 i j k l j k l i i (f f f f + fji fkl f f + 3fjk flj f l f k + fjkl f j f k f l) 24 j k l
+
h5 i j k l m j k l m j k i l m (f f f f f + fji fkj flm f l f m + 3fji fkl fm f f + fji fklm f k f l f m + 4fjk flj f k fm f 120 j k l m
(11.h)
j i i k l m i j k l m i + 4fjk flm f k f l f m + 3fjk flj fm f f + 6fjkl fm f f f + fjklm f j f k f l f m ) + O(h6 )
it results the following relations that should be fulfilled by the coefficients ars , br and cr h
h2
cr δ r = 1
(12.a)
cr br = 1/2
(12.b)
cr ars bs = 1/6 h
3
cr b2r = 1/3
(12.c)
cr ars ast bt = 1/24 cr ars b2s = 1/12 cr br ars bs = 1/8 h4
(12.d)
cr b3r = 1/4
cr ars ast atu bu = 1/120 cr ars ast b2t = 1/60 cr ars bs ast bt = 1/40 cr ars b3s = 1/20 cr br ars ast bt = 1/30 cr br ars b2s = 1/15 cr ars bs art bt = 1/20 cr b2r ars bs = 1/10 h5
cr b4r = 1/5 9
(12.e)
In these relations, br has been defined according to the property (10.b). Notice also that in the development of the aforementioned relations expansion in Taylor series were used only up to the term of fifth order (with h5 ). For higher order the deduction is very tedious. This is why the terms with h6 were not included in the series, though this term indicates the order of the method that this work is dealing with. Therefore, the relations (12) are valid for Runge-Kutta methods, whether implicit or explicit, from first order method (e.g. Euler method) to fifth order method (e.g. Fehlberg method [7]). In all the case the indexes r, s, t and u vary from one up to N , which is the number of stages. Gear [3], has presented a deduction similar to (11), but only for explicit methods. In Hairer et al. [6], relations similar to (12) appear, but only for explicit methods and only up to the term of order h4 . Also in Hairer et al. [6], a theorem that deduces the equivalence of the implicit Runge-Kutta methods and the orthogonal collocation method is presented. Ralston (1965) made a similar analysis to (11), to obtain the relations of the coefficients, but for an explicit Runge-Kutta method of fourth order and four stages, and found the following family of solutions of the relations (12), in term of b2 and b3 b1 = 0
b4 = 1
(13.a)
ars = 0
s≥r
(13.b)
a21 = b2
(13.c)
a31 = b3 − a32
(13.d)
b3 (b3 − b2 ) 2 b2 (1 − 2 b2 )
(13.e)
a41 = 1 − a42 − a43
(13.f )
a32 =
a42 =
a43 =
(1 − b2 )[b2 + b3 − 1 − (2 b3 − 1)2 ] 2 b2 (b3 − b2 )[6 b2 b3 − 4 (b2 + b3 ) + 3]
(13.g)
(1 − 2 b2 )(1 − b2 )(1 − b3 ) b3 (b3 − b2 )[6 b2 b3 − 4 (b2 + b3 ) + 3]
(13.h)
1 1 − 2(b2 + b3 ) + 2 12 b2b3
(13.i)
c2 =
2 b3 − 1 12 b2(b3 − b2 )(1 − b2 )
(13.j)
c3 =
1 − 2 b2 12 b3(b3 − b2 )(1 − b3 )
(13.k)
c1 =
10
c4 =
2 (b2 + b3 ) − 3 1 + 2 12 (1 − b2 )(1 − b3 )
(13.l)
Notice that if b2 = 1/2 and b3 = 1/2, then the well known classical Runge-Kutta method of fourth order is obtained 0 1/2
0 1/2
0 0
0 0
0 0
1/2 1
0 0
1/2 0
0 1
0 0
1/6 1/3 1/3 1/6
11
(14)
LOBATTO QUADRATURES
The explicit Runge-Kutta methods are of direct application, while the implicit Runge-Kutta methods require the resolution of a system of simultaneous equations with the variables kr in each step of integration of the differential equations, as it is suggested by the expression (6.b). This system of equations is generally not linear, unless the function f (x, y) be linear, and can be solved applying the iterative method of fixed point, explained farther on. The implicit Runge-Kutta method to be used here is a sixth order method (P = 6) with four stages (N = 4), developed on the bases of Lobatto quadrature [8] (for more details see [5] or [6]). The coefficients of this method expressed in the Butcher notation are 0 √ (5 − 5)/10 √ (5 + 5)/10 1
0 0 0 √ √ (5 + 5)/60 1/6 (15 − 7 5)/60 √ √ (5 − 5)/60 (15 + 7 5)/60 1/6 √ √ 1/6 (5 − 5)/12 (5 + 5)/12 1/12
5/12
5/12
0 0 0
(15.a)
0 1/12
This method will be named the principal implicit Runge-Kutta method. Within the coefficients of the principal implicit Runge-Kutta method it can be detected that a part of them forms another implicit Runge-Kutta method embedded in the first one. This method is of third order (P = 3), has three stages (N = 3) and in the Butcher notation can be expressed as
0 √ (5 − 5)/10 √ (5 + 5)/10
0 0 0 √ √ (5 + 5)/60 1/6 (15 − 7 5)/60 √ √ (5 − 5)/60 (15 + 7 5)/60 1/6 √ √ 1/6 (5 − 5)/12 (5 + 5)/12
(15.b)
Both methods, the principal and the embedded, form what is named the Lobatto embedding form of third and sixth orders with four stages. Notice that these methods are implicit only in the auxiliary variables k2 and k3 , and therefore the system of non-linear equations should be solved only in the mentioned variables. The other variables are of direct solution. In order to apply an iterative process to solve the system of non-linear equations, initial estimations of the values of the auxiliary variables are required. The best way to carry out the latter is to obtain these values from an explicit Runge-Kutta method, where the auxiliary variables kr are evaluated in the same intermediate point in each step, i.e. an explicit method having the same coefficient br of the implicit method. Observing the method (15.a), it is clear that the aforementioned explicit method is rapidly obtained from the 12
√ √ relations (13) assuming b1 = 0, b2 = (5 − 5)/10, b3 = (5 + 5)/10 and b4 = 1. Notice that these coefficients are consistent with the characteristics of an explicit method. This last aspect makes the implicit method (15.a) ideal for the desired purpose. Thus, if b2 and b3 are substituted in the relations (13), it is obtained that
a21 =
a31 = −
a32 =
√ (5 − 5) 10
(16.a)
√ (5 + 3 5) 20
(16.b)
√ (3 + 5) 4
(16.c)
These are the coefficients of the new explicit Runge-Kutta method that in the Butcher notation can be expressed as 0 √ (5 − 5)/10 √ (5 + 5)/10 1
0 0 √ (5 − 5)/10 0 √ √ −(5 + 3 5)/20 (3 + 5)/4 √ √ (−1 + 5 5)/4 −(5 + 3 5)/4 1/12
0
0
0
0
0 √ (5 − 5)/2
0 0
5/12
1/12
5/12
(17)
The obtained explicit Runge-Kutta method is not reported in the speciality literature and does not correspond to any known quadrature, but pertains to the family of solutions (13) of the fourth order, four stages explicit Runge-Kutta methods. This method will be used to obtain the initial estimations k2 and k3 for the iterative process in the following form k2,(0) = h f (xn + b2 h , yn + a21 k1 )
(18.a)
k3,(0) = h f (xn + b3 h , yn + a31 k1 + a32 k2 )
(18.b)
Once these initial estimations are obtained, a rapid convergence to the solution of the system of non-linear equations (6.b), with the coefficients (15.a), can be expected.
13
ITERATIVE PROCESS
As it was mentioned before, the system of non-linear equations (6.b), that is originated by any implicit Runge-Kutta method, can be solved in the auxiliary variables kr applying the iterative method of fixed point (see Gear[3]) i kr,(m+1) = h f i (xn + br h , yn + ars ks,(m) )
(19)
which is the easiest method to be used due to the form of the mentioned system of equations. Here m = 0, 1, 2, 3, . . . is the number of iteration in the iterative process. The global error in the iterative process is defined as i εir,(m) = kr,(m) − kri
(20.a)
where kri is the exact solution of the system of non-linear equations. The local error in the iterative process is defined as i i ir,(m) = kr,(m+1) − kr,(m)
(20.b)
cr r,(m) < max
(20.c)
and it is stopped when
i of the differential equations where max is the tolerance imposed to the local error to find the solution yn+1
in one step, and the norm of the local error r,(m) is supposed to be euclidean. If now expression (6.b) is substracted from expression (19) i kr,(m+1) − kri = h [f i (xn + br h , yn + ars ks,(m) ) − f i (xn + br h , yn + ars ks )]
(21)
Then, if the Lipschitz condition (with h > 0 for convenience) is applied, it results j i − kri | ≤ h lji |ars | |ks,(m) − ksj | |kr,(m+1)
(22)
|εir,(m+1) | ≤ h lji |ars | |εjs,(m) |
(23)
where lji is the maximum of the absolute value of each element in the jacobean matrix of f . This is |fji | ≤ lji Thus, if it is satisfied that ε(m) = max
max |εjs,(m) |
1≤j≤M 1≤s≤N
(24)
(25)
then |εir,(m+1) | ≤ h lji |ars | |εjs,(m) | ≤ h lji δj |ars | δs ε(m) 14
(26)
max
max |εir,(m+1) | ≤ max max h lji δj |ars | δs ε(m)
1≤i≤M 1≤r≤N
1≤i≤M 1≤r≤N
ε(m+1) ≤ h L A ε(m) where
(28)
L = max lji δj
(29.a)
A = max |ars | δs
(29.b)
1≤i≤M
and
(27)
1≤r≤N
The expression (28) means that for a high number of iterations, m → ∞, the global error ε(m) → 0 when
h≤
1 LA
(30)
and the iterative process is convergent locally (also globally) in the form cr r,(m+1) < cr r,(m)
(31)
The expression (30) is the limit of the stepsize, for the iterative method of fixed point to be convergent, when the system of non-linear equation is being solved in the implicit Runge-Kutta method. This is the only restriction of the implicit Runge-Kutta methods, compared to the explicit methods which are much less stable. In Appendix B a general formulation of an iterative process may be found, where the NewtonRaphson method, applied for any implicit Runge-Kutta method, is used instead of the fixed point procedure of expression (19). Thus the iterative procedure may converge more rapidly, but the additional numerical calculations must be evaluated, otherwise, the method may become inefficient. Sometimes the system of differential equations does not appear in the form (1), but in an implicit form of the type dy i dy = f i x, y, dx dx
i = 1, 2, 3, . . . , M
(32)
In this cases the iterative procedure is the same that was described before, but for the system of differential equations (32) it is necessary to define initial estimations for the auxiliary variables kr , including for the explicit Runge-Kutta method. These estimations should be close to the final solution, otherwise the number of iterations may be high. A form to obtain the above mentioned estimations is by making extrapolations with Lagrange polynomials. This is, to estimate k2 , it is made an extrapolation with a polynomial of degree 0 starting on k1 . To estimate k3 , it is made an extrapolation with a polynomial of degree 1 defined by k1 15
and k2 . To estimate k4 , an extrapolation is made with a polynomial of degree 2 defined by k1 , k2 and k3 . The procedure described yield the following results k2 = k1 b2 − b3 b3 k1 + k2 b2 b2
(33.b)
(b4 − b2 )(b4 − b3 ) b4 (b4 − b3 ) b4 (b4 − b2 ) k2 + k3 k1 + b2 b3 b2 (b2 − b3 ) b3 (b3 − b2 )
(33.c)
k3 =
k4 =
(33.a)
Notice that the coefficients br have been used as the points where the polynomials are defined. If expressions (33) are applied to the particular case of Runge-Kutta method of Lobatto type, it is finally obtained k2 = k1
k3
√ √ 1+ 5 3+ 5 k1 + k2 =− 2 2 k4 = k1 −
√
5(k2 − k3 )
(34.a)
(34.b)
(34.c)
The reader must bear in mind that the variable k4 , obtained one step before, may be used as an estimation of the initial point, necessary for the calculation of k1 .
16
STEPSIZE CONTROL
ERROR ANALYSIS
The implicit Runge-Kutta method of sixth order with four stages that is defined by the coefficients (15.a), in fact represents two different embedded methods, one of third order within the other of sixth order, i.e. the coefficients (15.b) are included within the coefficients (15.a). This aspect is relevant to control the stepsize, because solving the system of non-linear equations (6.b) for the same coefficients, it can be obtained two solutions of different order in the local truncation error, reducing to a minimum the number of numerical calculations to be made. Fehlberg [7] reported this aspect to control the stepsize in his explicit Runge-Kutta methods of fourth and fifth orders in a embedded form. ˜ n+1 be the solutions of the system of differential equations (1) or (2), offered by the Let yn and y Runge-Kutta methods type Lobatto of sixth and third orders, respectively, embedded in only one formulation as it was described before. This is 1 i (k + 5k2i + 5k3i + k4i ) 12 1
(35)
√ √ 1 [2k1i + (5 − 5)k2i + (5 + 5)k3i ] 12
(36)
i yn+1 = yni +
i y˜n+1 = yni +
The auxiliary variables k1 , k2 , k3 and k4 are the same for both expressions and are obtained using the equation (6.b) with the coefficients (15.a) and (15.b). i It will be denoted as En+1 the difference of the equation (35) minus the equation (36). This is
i i i En+1 = yn+1 − y˜n+1 =
√ 1 [−k1i + 5(k2i − k3i ) + k4i ] 12
(37)
If y(xn ) is the exact solution of the differential equation (2) in the value x = xn , the local truncation errors of the numerical solutions (35) and (36) are defined respectively as ein = yni − y i (xn ) = O(h7n−1 )
(38)
e˜in = y˜ni − y i (xn ) = O(h4n−1 )
(39)
i i i En+1 = yn+1 − y˜n+1 = ein+1 − e˜in+1 = O(h4n )
(40)
and then
Remember that the local truncation error is of order P + 1 if the Runge-Kutta method is of order P . 17
If the expression (40) is organized in the following form
i = En+1
yi
− y i (xn+1 ) i i yn+1 − y i (xn+1 )] y (xn+1 ) − [˜ y i (xn+1 )
n+1
(41)
it is obtained i = ei(r),n+1 y i (xn+1 ) − e˜in+1 En+1
(42)
where ei(r),n+1 =
yi
− y i (xn+1 )
i y (xn+1 )
n+1
(43)
is the relative local truncation error. i If now it is assumed that y i (xn+1 ) is approximated by yn+1 in the denominator of (43), it can be
applied the Cauchy-Schwartz and the triangular desequalities to the expression (42), and this results i i | ≤ |ei(r),n+1 | |y i (xn+1 )| + |˜ ein+1 | ≤ e(r),max |yn+1 | + e˜max |En+1
(44)
where e(r),max and e˜max are respectively the tolerances for the relative and absolute local truncation errors for the implicit Runge-Kutta methods of sixth and third orders. The expression (44) also means that, for the solution of the differential equations in one step be accepted, it should be verified that
Qin =
i | |En+1 ≤1 i e(r),max |yn+1 | + e˜max
(45)
being the tolerances for the relative and absolute local truncation errors proposed by the user of the algorithm.
CONTROL ALGORITHM Let hn+1 be the stepsize in the next step that tends to make Qin ∼ = 1. Taken into account the order i defined by (40), the parameter Qn may be redefined as of the difference En+1
Qn =
where
h 4 n hn+1
Qn = max Qin 1≤i≤M
18
(46)
(47)
and thus, solving for hn+1 , it is obtained
hn+1 = hn
1 1/4 = h n Sn Qn
with Sn =
1 1/4 Qn
(48)
(49)
Here it is convenient to mention that Shampine et al. [9] uses expressions similar to (48) and (49) to control the size of the step of integration in the Runge-Kutta method of fourth and fifth orders with five stages developed by Fehlberg [7], but with some modifications, in order to guarantee that Sn always be bound in the interval [Smin , Smax ], and that hn+1 always be greater than a limit value hmin . Additionally, the mentioned authors multiply Sn by a coefficient Cq less than unit, to make hn+1 tends to hn , and thus to make Qn ∼ = 1, but a little lower. All the aforementioned modifications are resumed in continuation as
Sn = Cq
1 1/4 Qn
Cq = 0.9 ∼ 0.99
(50)
Sn = max(min(Sn , Smax ), Smin )
(51)
hn+1 = hn Sn
(52)
hn+1 = max(hn+1 , hmin )
(53)
While in [9] the exponent is 1/5 in the expression (50), here the exponent is 1/4, and they also recommend for the coefficients and limits the values Cq = 0.9, Smin = 0.1 and Smax = 5. The value of the minimum stepsize, hmin , is determined by the precision of the computer to be used. In this work they are recommended the same values for the coefficients in the expressions (50) to (53). The procedure to calculate the optimal value of the size of the integration step, that permits to satisfy the tolerances e(r),max and e˜max , is described in continuation: • Estimated an initial stepsize hn , the implicit Runge-Kutta method type Lobatto is used to calculate the auxiliary variables k1i , k2i , k3i and k4i with the expression (6.b), using the coefficients (15.a) and with the iterative process (19), using the initial values (18). i i • The expressions (35) and (36) permit to find the solutions yn+1 and y˜n+1 of the methods of sixth
and third orders, respectively. i • The definition (37) permits to calculate the difference En+1 between both methods.
19
• With the equation (45) it can be calculated the parameters Qin , and with the equation (47) it can be obtained the maximum of them. • The relations (50) to (53) determine the value of the size of the next step hn+1 . • If Qn ≤ 1, the integration with the step hn (or the application of the Runge-Kutta method from xn to xn+1 ) is accepted and the step hn+1 is considered the step for the next integration (or the next application of the Runge-Kutta method from xn+1 to xn+2 ). • If Qn > 1, the integration with the stepsize hn is rejected and it is repeated all the algorithm but with hn = hn+1 obtained from (53). This procedure, sometimes increases the stepsize, and other times decreases the stepsize, in a optimal form, in order to guarantee that the relative error ei(r),n+1 of the sixth order Runge-Kutta method be less than the tolerance e(r),max , and the error e˜in+1 of the third order Runge-Kutta method be less than the i , i.e. the solution with the tolerance e˜max . In any case, the solution of the Runge-Kutta method will be yn+1
implicit Runge-Kutta method of sixth order.
20
METHOD’S ANALYSIS
PRECISION
The order of precision of any Runge-Kutta method comes from the comparison between this and the expansion in Taylor series of y(xn+1 ) around y(xn ). From this comparison was obtained the relations (12) that should be fulfilled by the coefficients of the Runge-Kutta methods. For a Runge-Kutta method of order P , it should be satisfied the relations (12) up to the term of the Taylor series that contains hP . The remainder terms in the Taylor series are they which determine the order of the local truncation error. Thus, a method of order P has a local truncation error of order P + 1, i.e. that depends on hP +1 . The global truncation error always is one order less than the local truncation error [3], i.e. that depends on hP . All these criterions can be applied to the implicit Runge-Kutta type Lobatto of third and sixth order. In this form it is obtained that, the method of third order satisfy the relations (12) up to term with h3 and the local and global truncation errors depend, respectively, on h4 and h3 . The method of sixth order satisfy the relations (12) up to term with h6 (the relations for this last term does not appear for the reasons explained there) and the local and global truncation depend, respectively, on h7 and h6 . The dependence of an error respect to hP is indicated as O(hP ) and is said that the order of precision is P , according to the definition (7). Also it is satisfied that O(hP ) + O(hQ ) = O(hP ) when Q ≥ P .
CONVERGENCY
In the section of Iterative Process it was indicated that the implicit Runge-Kutta methods generated a system of equations of the type (6.b), in general of non-linear characteristics, where the unknowns were the auxiliary variables kr . As it was said before, this system of non-linear equations may be solved using an iterative method of fixed point, which converges for the treated problem in particular, if it is satisfied that
h≤ where
1 LA
(30)
L = max lji δj
(29.a)
A = max |ars | δs
(29.b)
1≤i≤M
1≤r≤N
21
and where the local convergence is established according to cr r,(m+1) < cr r,(m)
(31)
For the specific case when it is being solved the problem with only one differential equation of the form dy = f (y) dx
f (y) = λy
(54)
it is obtained that L = |λ|. For the implicit Runge-Kutta methods type Lobatto of third and sixth orders, it √ is obtained that A = (5 + 5)/10 and A = 1, respectively. Thus, to assure the convergence of the iterative process stated in these cases, particularly to the linear problem (54), it should be satisfied the following two conditions (assuming h positive)
h≤
√ 5 − 5 ∼ 1.38 = 2 |λ| |λ|
h≤
1 |λ|
(Implicit 3rd order method)
(Implicit 6th order method)
(55)
(56)
respectively. Notice that the condition (56) is the most restrictive of the two.
STABILITY
The stability of the Runge-Kutta methods is generally studied on the bases of their performance for solving the specific problem (54) with only one linear ordinary differential equation. In Appendix C it can be found the analysis of the stability for the more general case of a system of linear ordinary differential equations. Thus, if the expression (6.b) is applied, taken into account that the problem being solved is (54), it is obtained that kr = h f (yn + ars ks ) = h λ (yn + ars ks ) = h λ yn + h λ ars ks
(57)
Moreover, if kr is substituted as δrs ks and the therms are regrouped, it results δrs ks = h λ yn + h λ ars ks
(58.a)
δrs ks − h λ ars ks = h λ yn
(58.b)
and thus, if ks is factorized, then [δrs − h λ ars ] ks = h λ yn δr 22
r, s = 1, 2, 3, . . . , N
(59)
This last expression represent a system of N linear equations with N unknowns ks . If this system of linear equations is solved and if the solutions kr are substituted in the equation (6.a), then, it is found a relation for yn+1 depending only on yn and on the coefficients of the used Runge-Kutta method. The mentioned relation is of the form yn+1 = µ1 (hλ) yn
(60)
where µ1 (hλ) is found applying, for example, to the implicit sixth order Runge-Kutta method of Lobatto type, the procedure above explained. In this case, it is obtained
1 3 z + 1 + 23 z + 15 z 2 + 30 µ1 (z) = 1 1 2 1 − 3 z + 30 z
1 4 360 z
(61)
where the function µ1 (z) was deduced from the coefficients (15.a). For the case of the implicit third order Runge-Kutta method type Lobatto resumed in the coefficients (15.b), it is obtained y˜n+1 = µ ˜1 (hλ) yn
(62)
where
1 3 1 + 23 z + 15 z 2 + 30 z µ ˜1 (z) = 1 1 2 1 − 3 z + 30 z
(63)
The functions µ ˜1 (z) and µ1 (z) are denominated characteristic roots of the Runge-Kutta methods of third and sixth orders, respectively. The Runge-Kutta methods, to which these roots pertain, are considered stables if their absolute values are less than unit, for a determined value of z = hλ. Notice that, if |˜ µ1 (hλ)| yn+1 | or |yn+1 | is less than |yn | and the stability is or |µ1 (hλ)| is less than the unit, then it is satisfied that |˜ guaranteed. The figure 1 represent graphically the functions µ ˜1 (hλ) and µ1 (hλ) vs. hλ, and the region of stability for each one. There it can be observed that the methods are stable, if it is satisfied the following two conditions (assuming that h is positive and λ is negative).
h≤
6.8232 |λ|
(Implicit 3rd order method)
(64)
h≤
9.6485 |λ|
(Implicit 6th order method)
(65)
Comparing the conditions (55) and (56) with the conditions (64) and (65), respectively, the condition (56) continues being the most restrictive of all the mentioned conditions. 23
The function µ1 (z) constitutes an aproximation of Pad´e [2] for the function y = ex (see Lapidus and Seinfeld [5]), and, additionally, is always positive and less than unit in the interval [−9.648495252 , 0.0]. Notice in the upper graphic of the figure 1 that the function µ1 (z) approximates well to the function y = ez for the range z > −4. The function µ ˜1 (z), however, is not an approximation of Pad´e, has only one root in the point z = −2.706010973 . . ., and is, in absolute value, less or equal to the unit within the interval [−6.823183583 , 0.0]. This can be observed in the lower graphic of the figure 1. The conditions (64) and (65) reveal that the implicit Runge-Kutta methods type Lobatto of third and sixth orders are more stable than the explicit Runge-Kutta methods type Fehlberg of fourth and fifth orders, having these last methods the following stability conditions
h≤
2.785 |λ|
(Explicit 4th order method)
(66)
h≤
3.15 |λ|
(Explicit 5th order method)
(67)
The only limitation of the implicit Runge-Kutta methods comes from the convergency conditions (55) and (56), which are based on the way used to solve the system of equations (6.b), and on the assumption that the differential equation has a linear form (54). These convergency conditions permit an increase of the stepsize h much less than the stability conditions (64) and (65). This aspect brakes the advance of the implicit Runge-Kutta methods in a notable manner, but the difficulty may be compensated partially in two forms. First, estimating conveniently the initial values of the auxiliary variables kr for the iterative process. As it was treated before, this can be made using an explicit Runge-Kutta method with the expressions (18) and/or making polynomial extrapolations with the expressions (34). Second, it may be used the fixed point method in a efficient manner, similar to the Gauss-Seidel method, i.e. when an unknown is approximately calculated, it is substituted inmediately in the next equation, and thus on. These two modifications of the implicit Runge-Kutta methods here used may improve the performance of the method, in the sense of that the stepsize h is liberated from the convergency conditions (55) and (56) and thus it may be permitted to increase much more. A more sofisticated solution to the problem mentioned before can be developed if it is used the Newton-Rapson method to solve the system of non-linear equations (6.b) as it is proposed in Appendix B. This may increase the number of numerical calculations to be performed by the algorithm, but the convergence will be more rapid.
24
Figure 1. Stability regions for implicit methods. Curves of µ1 vs. h λ
25
REFERENCES
[1] Gerald, C. F. Applied Numerical Analysis. 2nd Edition. Addison-Wesley, 1970. [2] Burden R. L.; Faires, J. D. Numerical Analysis. 3rd Edition. PWS. Boston, 1985. [3] Gear, C. W. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, 1971. [4] Butcher, J. C. “On the Runge-Kutta Processes of High Order”. J. Austral. Math. Soc., Vol. IV, Part 2, 1964, pp. 179-194. [•] Butcher, J. C. The Numerical Analysis of Ordinary Differential Equations, Runge-Kutta and General Linear Methods. John Wiley, New York, 1987. [5] Lapidus, L.; Seinfeld, J. H. Numerical Solution of Ordinary Differential Equations. Academic Press, 1971. [6] Hairer, E.; Norsett, S. P.; Wanner, G. Solving Ordinary Differential Equations I. Nonstiff Problems. Springer-Verlag, 1987. [•] Hairer, E.; Wanner, G. Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems. Springer-Verlag, 1991. [7] Fehlberg, E. “Low-Order Classical Runge-Kutta Formulas with Stepsize Control”. NASA Report No. TR R-315, 1971. [8] Lobatto, R. Lessen over Differentiaal- en Integraal-Rekening. 2 Vol. La Haye, 1851-52. [9] Shampine, L. F.; Watts, H. A.; Davenport, S. M. “Solving Non-Stiff Ordinary Differential Equations - The State of the Art”. SANDIA Laboratories, Report No. SAND75-0182, 1975. SIAM Review, Vol. 18, No. 3, 1976, pp. 376-411.
26
APPENDIXES
Appendix A. Deduction of the relations fulfilled by the coefficients in the implicit Runge-Kutta method. The coefficients ars , br and cr are determined applying their properties and using some relations that are deduced in the following form: Let the system of ordinary differential equations be expressed as an initial value problem of the type
dy = f (y) dx x = x0
y(x0 ) = y0
(A.1) (A.2)
The Runge-Kutta method, applied to this problem, is formulated as yn+1 = yn + cr kr
(A.3)
where the auxiliary variables kr are defined as kr = h f (yn + ars ks )
(A.4)
and where the coefficients have the next properties 0 ≤ br ≤ 1
(A.5)
ars δs = br
(A.6)
cr δ r = 1
(A.7)
If now the expansion in Taylor series is made to the component kri of (A.4), around the point (xn , yn ), being yn = y(xn ), it results that
kri =h f i [δr ] + h fji [ars ksj ] +
h i f [ars ksj ] [art ktk ] 2 jk
+
h i f [ars ksj ] [art ktk ] [aru kul ] 6 jkl
+
h i f [ars ksj ] [art ktk ] [aru kul ] [arv kvm ] + O(h6 ) 24 jklm 27
(A.8)
where the index rules and the following notation has been used
f i = f i (xn )
fji =
∂f i ∂y j yn
i fjk =
∂ 2 f i ∂y j ∂y k yn
···
(A.9)
Here the functions are supposed to be C ∞ (Analytical Functions) and therefore the subindexes in (A.9) are permutable. The variable ksj in the second term of the right side of (A.8) can again be expanded in Taylor series as ksj =h f j [δs ] + h fkj [asα kαk ] +
h j f [asα kαk ] [asβ kβl ] 2 kl (A.10)
h j + fklm [asα kαk ] [asβ kβl ] [asγ kγm ] + O(h5 ) 6 but up to the term that is one order less. In the same way kαk can be expanded as
kαk = h f k [δα ] + h flk [aαδ kδl ] +
h k f [aαδ kδl ] [aα km ] + O(h4 ) 2 lm
(A.11)
and thus sucessively l kδl = h f l [δδ ] + h fm [aδϕ kϕm ] + O(h3 )
(A.12)
kϕm = h f m [δϕ ] + O(h2 )
(A.13)
up to
If now a recurrent backward substitution is made, then, it is obtained that l kδl = h f l [δδ ] + h fm aδϕ [h f m δϕ + O(h2 )] l m = h f l [δδ ] + h2 fm f bδ + O(h3 )
l m kαk = h f k [δα ] + h flk aαδ [h f l δδ + h2 fm f bδ ] +
h flm [aαδ hlδ ] [aα km ] + O(h4 ) 2
l m f aαδ bδ + = h f k [δα ] + h2 flk f l bα + h3 flk fm
h flm aαδ [h f l δδ ] aα [h f m δ ] + O(h4 ) 2
l m = h f k [δα ] + h2 flk f l bα + h3 flk fm f aαδ bδ +
h3 k l m 2 f f f bα + O(h4 ) 2 lm
28
(A.14)
(A.15)
l m ksj =h f j [δs ] + h fkj asα [h f k δα + h2 flk f l bα + h3 flk fm f aαδ bδ +
+
h j l m f asα [h f k δδ + h2 fuk f u bα ] asβ [h f l δβ + h2 fm f bβ ] 2 kl
+
h j f asα [h f k δα ] asβ [h f l δβ ] asγ [h f m δγ ] + O(h5 ) 6 klm
h3 k l m 2 f f f bα ] 2 lm
l m =h f j [δ] + h2 fkj f k bs + h3 fkj flk f l asα bα + h4 fkj flk fm f asα aαδ bδ
+
h4 j k l m h3 j k l 2 h4 j k l m fk flm f f asα b2α + fkl f f bs + fkl f fm f bs asβ bβ 2 2 2
+
h4 j k m l h4 j f k f l f m b3s + O(h5 ) fkl fm f f bs asα bα + fklm 2 6
(A.16)
and [ars ksj ] =h f j br + h2 fkj f k ars bs + h3 fkj flk f l ars asα bα +
l m + h4 fkj flk fm f ars asα aαβ bβ +
j k l m + h4 fkl fm f f ars bs asα bα +
h3 j k l f f f ars b2s 2 kl
h4 j k l m f f f f ars asα b2α 2 k lm
h4 j f f k f l f m ars b3s + O(h5 ) 6 klm
(A.17)
where j k l m j k l m fkl f fm f = fkl fm f f
(A.18)
because the function f (y) is analytical. Thus the term h fji [ars ksj ] results as
=h2 fji f j br + h3 fji fkj f k ars bs + h4 fji fkj flk f l ars asα bα
l m + h5 fji fkj flk fm f ars asα aαβ bβ +
h5 i j k l m f f f f f ars asα b2α 2 j k lm
j k l m + h5 fji fkl fm f f ars bs asα ba lpha +
29
h4 i j k l f f f f ars b2s 2 j kl
h5 i j f f f k f l f m ars b3s + O(h6 ) 6 j klm
(A.19)
The term
h 2
i fjk [ars ksj ] [art ktk ] results as
h i h3 j v w = fjk [h f j br + h2 fvj f v ars bs + h3 fvj fwv f w ars asα bα + fvw f f ars b2s + O(h4 )] 2 2
l m . [h f k br + h2 flk f l art bt + h3 flk fm f art atβ bβ +
h3 k l m f f f art b2t + O(h4 )] 2 lm
h i 2 j k 2 h4 l m k = fjk [h f f br + h3 f j flk f l br art bt + h4 f j flk fm f br art atβ bβ + f j flm f l f m br art b2t 2 2
+ h3 fvj f v f k br ars bs + h4 fvj f v flk f l ars bs art bt + h4 fvj fwv f w f k br ars asα bα
+
h4 j v w k f f f f br ars b2s + O(h5 )] 2 vw
h i 2 j k 2 l m [h f f br + 2 h3 flj f k f l br ars bs + 2 h4 flj f k fm f br ars asα bα = fjk 2
j k l m + h4 flm f k f l f m br ars b2s + h4 flj fm f f ars bs art bt + O(h5 )]
=
h3 i j k 2 i i l m f f f br + h4 fjk flj f k f l br ars bs + h5 fjk flj f k fm f br ars asα bα 2 jk
+
h5 i j k l m h5 i j k l m fjk flm f f f br ars b2s + fjk fl fm f f ars bs art bt + O(h6 ) 2 2
(A.20)
where i i fjk f j flk f l = fjk flj f k f l i l m i l m fjk f j flk fm f = fjk flj f k fm f j i k i fjk f j flm f l f m = fjk flm f k f lf m
30
(A.21)
The term
h 6
i fjkl [ars ksj ] [art ktk ] [aru kul ] results as
=
h i j k m f [h2 f j f k b2r + 2 h3 fm f f br ars bs + O(h4 )] 6 jkl
. [h f l br + h2 fwl f w aru bu + O(h3 )]
=
h4 i j k l 3 h5 i j k l w 2 f f f f br + fjkl f f fw f br aru bu 6 jkl 6
+
=
2 h5 i j k m l 2 f f f f f br ars bs + O(h6 ) 6 jkl m
h4 i j k l 3 h5 i j k l m 2 f f f f br + fjkl fm f f f br ars bs + O(h6 ) 6 jkl 2
(A.22)
where i i l m i k m fjkl f j f k fwl f w = fjkl f j f k fm f = fjlk f j f l fm f i k l m i k l m = fjlk f j fm f f = fjkl f j fm f f
(A.23.a)
i j k m l i j k l m i k j l m fm f f f = fjkl fm f f f = fkjl fm f f f fjkl i k l m i k l m = fkjl f j fm f f = fjkl f j fm f f
(A.23.b)
and where b2r aru bu = b2r ars bs . The term
h 24
i fjklm [ars ksj ] [art ktk ] [aru kul ] [arv kvm ] results as
=
h i f [h3 f j f k f l b3r + O(h4 )] [h f m br + O(h2 )] 24 jklm
=
h5 i f f j f k f l f m b4r + O(h6 ) 24 jklm 31
(A.24)
Finally an expression for kri is obtained as 1 i j k 2 f f br ] kri =h f i [δr ] + h2 [fji f j br ] + h3 [fji fkj f k ars bs + fjk 2 1 j k l + h4 [fji fkj flk f l ars ast bt + fji fkl f f ars b2s 2 1 i j k l 3 i + fjk flj f k f l br ars bs + fjkl f f f br ] 6 1 l m k f ars ast atu bu + fji fkj flm f l f m ars ast b2t + h5 [fji fkj flk fm 2 (A.25) 1 j k l m j + fji fkl fm f f ars bs ast bt + fji fklm f k f l f m ars b3s 6 1 i j k l m i l m + fjk flj f k fm f br ars ast bt + fjk flm f f f br ars b2s 2 1 i j k l m 1 i j k l m 2 + fjk fl fm f f ars bs art bt + fjkl fm f f f br ars bs 2 2 +
1 i f f j f k f l f m b4r ] + O(h6 ) 24 jklm
If this last expression of the components of kr is inserted into the equation (A.3), and this is compared to the expansion in Taylor series of yn+1 around yn
i yn+1 =yni + h f i +
h2 i j h3 i (fj f ) + (fji fkj f k + fjk f j f k) 2 6
+
h4 i j k l j k l i i (f f f f + fji fkl f f + 3fjk flj f l f k + fjkl f j f k f l) 24 j k l
+
h5 i j k l m j k l m j k i l m (f f f f f + fji fkj flm f l f m + 3fji fkl fm f f + fji fklm f k f l f m + 4fjk flj f k fm f 120 j k l m
(A.26)
j i i k l m i j k l m i + 4fjk flm f k f l f m + 3fjk flj fm f f + 6fjkl fm f f f + fjklm f j f k f l f m ) + O(h6 )
it results the following relations that should be fulfilled by the coefficients ars , br and cr h
h2
cr δ r = 1
(A.27)
cr br = 1/2
(A.28)
32
cr ars bs = 1/6 h3
cr b2r = 1/3
(A.29)
cr ars ast bt = 1/24 cr ars b2s = 1/12 cr br ars bs = 1/8 h4
(A.30)
cr b3r = 1/4
cr ars ast atu bu = 1/120 cr ars ast b2t = 1/60 cr ars bs ast bt = 1/40 cr ars b3s = 1/20 cr br ars ast bt = 1/30 cr br ars b2s = 1/15
(A.31)
cr ars bs art bt = 1/20 cr b2r ars bs = 1/10 h5
cr b4r = 1/5
In these relations br has been defined according with the properties (A.5), (A.6) and (A.7). Notice also that in the developing of the mentioned relations were used expansions in Taylor series only up to the term with h5 . For higher order the deduction is very tedious, but it can be made similarly as it was made in this appendix. As an example, there exist the full implicit Runge-Kutta methods based on the Gauss-Legendre quadrature. In these cases, additionally to the expressions (A.27) to (A.31), the coefficients have to fulfill the following relations ars bγ−1 = s
bγr γ
(A.32)
cs bγ−1 = s
1 γ
(A.32)
γ = 1, 2, 3, . . . , N and the coefficients br are found from the roots of the Legendre polynomial of order N , the number of stages, i.e. PN (2 br − 1) = 0
(A.33)
where the order of the Runge-Kutta method originated is twice the number of stages of the method (P = 2N ). 33
In the Butcher notation, the three first full implicit Runge-Kutta methods, based on the GaussLegendre quadrature, with 1, 2 and 3 stages, are respectively
1/2
1/2
(A.34)
1 √ 3)/6 √ (3 + 3)/6
(3 −
√ 1/4 (3 − 2 3)/12 √ (3 + 2 3)/12 1/4 1/2
(5 −
√ 15)/10
1/2 √ (5 + 15)/10
(A.35)
1/2
√ √ 5/36 (10 − 3 15)/45 (25 − 6 15)/180 √ √ (10 + 3 15)/72 2/9 (10 − 3 15)/72 √ √ (25 + 6 15)/180 (10 + 3 15)/45 5/36 5/18
4/9
34
5/18
(A.36)
Appendix B. Newton-Raphson method applied to solve the non-linear system of equations with auxiliary variables kri in the implicit Runge-Kutta methods. Let a system of ordinary differential equations be expressed as
dy i = f i (y) dx
(B.1)
with the initial conditions y i (xo ) = yoi
x = xo
(B.2)
For solving this initial value problem, the implicit Runge-Kutta method i yn+1 = yni + cr kri
(B.3)
kri = h f i (yn + ars ks )
(B.4)
can be used. The expression (B.4) should be interpreted as a system of non-linear equations with the auxiliary variables kri as unknowns. For this reason it is convenient to define a function gri (k) = kri − h f i (yn + ars ks )
(B.5)
that should be zero in each component when the solution for each kri has been found. In order to solve the system of non-linear equations (B.5), it is more efficient to use the NewtonRaphson method rather than using the easier fixed point method which can be applied with the expression (B.4) in the form i kr,(m+1) = h f i (yn + ars ks,(m) )
(B.4 )
However, for using the Newton-Raphson method, the jacobean tensor of the function gri (k), with respect to the variables ktj , has to be defined as
∂gri ∂ktj
= δ ij δrt − h
∂f i ars δ kj δst ∂y k yn +ars ks
= δ ij δrt − h fji (yn + ars ks ) art
(B.6)
Thus the Newton-Raphson method can be applied in the following manner i i i kr,(m+1) = kr,(m) + ω ∆kr,(m)
35
(B.7)
j where the variables ∆kt,(m) are found from the next system of linear equation
∂g i
r ∂ktj (m)
j ∆kt,(m) = −gri (k(m) )
(B.8)
and ω is a relaxation factor. The iterative process (B.7) and (B.8) is applied in a succesive form until cr r,(m) < max cr ∆kr,(m) <
max ω
(B.9) (B.9 )
where i i ir,(m) = kr,(m+1) − kr,(m)
(B.10)
is the local error in the auxiliary variables kri and where max is the tolerance for the local truncation error i in the numerical solution yn+1 .
For very complicated functions, it is convenient to express the partial derivatives of the jacobean tensor (B.6) in a numerical approximated form, as it is indicated in the following expression ∂g i
r ∂kti
f i (y + ∆y ) − f i (y(j) )
n n n ∼ art = δ ij δrt − h ∆ynj
(B.11)
where ∆ynj = ars ksj
(B.12)
f i (y(j) ) = f i (yn1 + ∆yn1 , yn2 + ∆yn2 , . . . , ynj , . . . , ynM + ∆ynM )
(B.13)
and
36
Appendix C. Stability of the Runge-Kutta methods for the general linear case. Let a system of ordinary differential equations be expressed as dy i = f i (y) dx
(C.1)
f i (y) = λij y j
(C.2)
where
is of linear dependency, i.e. the parameters λij are constant. Let now the Runge-Kutta method be applied to this particular case. The first step is for obtaining the auxiliar variables kri in the form kri = h f i (yn + ars ks ) = h λij [ynj δr + ars ksj ] = h λij ynj δr + h λij ars ksj
(C.3)
where expressing kri = δ ij hjs δsr it is obtained that δ ij ksj δsr = h λij ynj δr + h λij ars ksj
(C.4.a)
δ ij ksj δsr − h λij ars ksj = h λij ynj δr
(C.4.b)
δ ij δrs ksj − h λij ars ksj = h λij ynj δr
(C.4.c)
[δ ij δrs − h λij ars ] ksj = h λij ynj δr
(C.4.d)
ij , such that If each member in the last equation is multiplied by a tensor Rrs ij jk Rrs [δ δst − h λjk ast ] = δ ik δrt
(C.5)
[δ jk δst − h λjk ast ] ktk = h λjk ynk δs
(C.6)
then it results that
ij jk ij Rrs [δ δst − h λjk ast ] ktk = Rrs h λjk ynk δs
(C.7.a)
ij δ ik δrt ktk = Rrs h λjk ynk δs
(C.7.b)
ij kri = Rrs h λjk ynk δs
(C.7.c)
Notice that the equation (C.6) is the same expression as equation (C.4.d), but with the indexes changed. This was made in order to obtain the appropiated indexes in the result (C.7.c). If then, the obtained expression for kri is substituted in the second step of the Runge-Kutta method, it results i = yni + cr kri yn+1
37
(C.8)
i ij yn+1 = yni + cr Rrs h λjk ynk δs
(C.9.a)
ij = δ ik ynk + cr Rrs h λjk ynk δs
(C.9.b)
ij = [δ ik + cr Rrs h λjk δs ] ynk
(C.9.c)
Thus i yn+1 = µik ynk
(C.10)
ij µik = [δ ik + cr Rrs h λjk δs ]
(C.11)
where
are the characteristic roots of the system of ordinary differential equations (C.1) proposed. If the norm of the expression (C.10) is finally obtained as max yn+1 ≤ max |µik | δ k ynmax i
(C.12)
where ynmax = max |ynk |
max k yn+1 = max |yn+1 |
k
k
(C.13)
then,it is deduced that, the norm of µik must be less than unit, i.e. µ = max |µik |δ k < 1, i
for the method to be stable.
38
(C.14)
Appendix D. Listing of the subroutine “RKI36” to solve systems of ordinary differential equations. The code has been deleted intentionally. If you want to seek it, write to the author to the e-mail:
[email protected]
39
Appendix E. An illustrating example: The Arenstorf orbits. Let the methods described in this work be illustrated with an example from Astronomy, the restricted three body problem. Consider two bodies of masses η and µ in circular rotation in a plane and a third body of negligible mass moving around in the same plane. The differential equations in relative variables for the case of the moon and the earth are dx =u dt
x+µ du x−η = x + 2v − η −µ dt A B (E.1) dy =v dt
y dv y = y − 2u − η −µ dt A B
where [(x + µ)2 + y 2 ]3 B = [(x − η)2 + y 2 ]3 A=
µ = 0.012277471 η =1−µ
(E.2)
The initial values have been determined carefully
x(0) = 0.994
u(o) = 0
y(0) = 0
v(0) = −2.00158510637908252240537862224
(E.3)
so that the solution becomes periodic with period T = 17.0652165601579625588917206249
(E.4)
Such periodic solutions have fascinated astronomers and mathematicians for many decades (Poincar´e) and are now often called “Arenstorf orbits” in honour to Arenstorf (1963), who also did many numerical computations on high speed electronic computers. The problem is C ∞ with the exception of the two singular points x = −µ and x = η, y = 0, therefore the Euler polygons are known to converge to the solution. But are they really numerically useful here? It has been chosen ns = 24000 steps of step length h = T /ns and plotted the result in figure E.1. The result is not very striking. The performance of the four order (and four stages) explicit Runge-Kutta method is already much better and converges faster to the solution. It has been used ns = 6000, so that the numerical work becomes 40
equivalent. Clearly, most accuracy is lost in those parts of the orbit which are close to a singularity. Therefore, codes with automatic step size selection perform much better. For example, the code DOPRI5 computes the orbit with a precision of 10−3 in 75 steps. The step size becomes very large in some regions and the grafical representation as poligons connecting the solution points becomes unsatisfactory. The solid line is the interpolatory solution, which is also precise for all intermediate values and useful for many other questions such as delay differential equations, dense output or implicit output. The figure E.2 shows the results computed by Lobatto implicit sixth order Runge-Kutta method (RKI36) and by Fehlberg fifth order Runge-Kutta method (RKF45). In both cases it has been chosen ns = 2000 steps by orbit to printout the results, but using automatic stepsize control (with a tolerance of 10−12 and 10−10 for absolute and relative errors, respectively). Thus a total of five orbits were obtained. The upper left graphic shows the solution of the differential equations in three or less orbits with both methods while the upper right graphic shows that RKF45 is losing accuracy in the orbit four. The lower graphics shows the results in the orbit five for RKI36 and RKF45, respectively. It is clear that RKF45 continues the orbit four with a catastrophic orbit which is really far from the exact solution, like Euler method in the figure E.1. Contrarily to this RKI36 continues the orbit four with an orbit close to the exact solution and thus follows in a stable manner. Speaking in absolute terms, the method RKI36 with automatic stepsize control perform as good as the code DOPRI5, computing 75 steps in an orbit with a tolerance in the absolute local error of 10−3 . Finally for a time t = T it was obtained almost the initial location with errors of ∆x = 8∗ 10−5 and ∆y = 3∗ 10−3 . However, the method RKF45 with automatic stepsize control computed the results in 44 steps with errors of ∆x = −3∗ 10−1 and ∆y = −3∗ 10−1 for the final location, using the same tolerance. This mean that RKF45 is more rapid than RKI36, aproximately twice, but much less stable. It is interesting to mention that for the method RKI36 with equidistant step size, i.e. without the automatic stepsize control, to reach the final location at t = T with errors of ∆x = −2∗ 10−2 and ∆y = 8∗ 10−4 , it was necessary approximately 6000 steps. This indicates what efficient is the method when automatic stepzise control is used. The conclusions for the analysis of the aforementioned results are: The implicitness of the method makes the performance more stable. The algorithm of automatic stepsize control makes the method more efficient.
41
Figure E.1. An Arenstorf orbit computed by equidistant Euler, equidistant Runge-Kutta and variable step size Dormand and Prince (DOPRI5).
42
Figure E.2. Arenstorf orbits computed by implicit (RKI36) and explicit (RKF45) Runge-Kutta methods with ns = 2000 steps by orbit (for printout only) with automatic stepsize control.
43
