RUNGE—KUTTAULASWITHTHEABILITYOFERRORESTIMATION.doc

RUNGE-KUTTA FORMULAS WITH THE ABILITY OF ERROR ESTIMATION ESTIMATIONBy Masatsugu TanakaAssistant ProfessorFaculty of Engineering, University of Yamanashi1. PreferenceSome trials to give the error-estimating ability to Runge-Kutta methods have been made by several researchers:a. As to the fourth-order methods.(1) The research of R. E. Scraton 1.(2) The research of R. H. Merson 2, 3, 4, 5. This is so-called Kutta-Merson Process. In this paper we take into consideration an ordinary differential equation dy/dx=f(x,y) with the initial consideration y(x0)=y0. Merson dealt with the case only where it is linear. But it is known that if it is linear, the truncation error of classical Runge-Kutta formula is obtained without sacrificing the computations of the function.b. As to the third-order methods(1) The research of R. H. Merson above stated. In spite of the misunderstanding of the deviser and introducers of this process, this method is in itself one of the third-order methods in the case of non-linear functions. (see section 4 and 6)(2) The research of F. Ceschino 6c. As to the second-order methods.(1) The research of F. Ceschino 6. In both cases Ceschino adds two more computations to the ordinary methods.In this paper, the author studies second, third and fourth-order formulas which have the ability of error estimation. To develop this ability the author, in view of truncation error, optimizes the paramaters which characterize runge-kutta formulas. And here the accuracy of truncations in each formula is evaluated by the measures independent of the function f.In the case of second and third-order methods the author computes the function only while the above-mentioned researchers do twice for the error evaluation. In the case of fourth-order methods, the author makes approximately a fourth-order method for the general function without making kutta-merson process so complicated as R. E. Scraton did, but only by optimizing the coefficients. This is the way which merson tried to devise but failed.Further ,author shows in this paper the fact that it is possible to obtain the accuracy of fifth-order method when five function values are used, though F. Ceschino has already proved it impossilble to make fifth-order methods in a strict sence 7.2. Second-order FormulasThe set of formulask1=hf(x0, y0) (2.1)k2=hf(x0+h, y0+k1) (2.2)k3=hf(x0+1h, y0+1k1+1k2) (2.3)yx0+h=yx0+i=12iki (2.4)T= i=13iki (2.5)Is the general form of the second-order method, where i, i, i, i and i are scalars. T denotes the estimation of the truncation error of (2.4).In this section, we determine the coefficients so as to make (2.4) represent a second-order formula and (2.5) its truncation error. For this purpose we expand k2 and k3into series:k2=hf+hD1f+h22!D12f+h33!D13f+h44!D14f+ (2.6) k3=hf+hD2f+h22!D22f+h33!D23f+h44!D24f+h21fyD1f+h2!fyD12f+ hD1fD2fy+h23!fyD13f+h22!D12fD2fy+h22!1fy(D1f)2+h22!D1fD22fy+ (2.7)WhereD1=x+f0y (2.8)D2=1x+(1+1)f0y (2.9)And the zero subscript means that the functions are evaluated at the point (x0,y0). On the order hand, if we setD=x+f0y (2.10)Then yx0+h=yx0+hf+h22!D2f+h33!(D2f+fyDf+h44!(D3f+fyD2f+fy2Df+ 3DfDfy)+h55!(D4f+6DfyD2f+4D2fyDf+D2ffy2+Dffy3+3(Df)2fy+ D3ffy+7fyDfDfy)+ (2.11)In the order that (2.4) can be a second-order formula, the right-hand side of (2.4) must agree completely with the right-hand side of (2.11) through the terms up to order h2. Thus the following equations are obtained:1+2=1 (2.12)2D1f=12!Df (2.13)The truncations error of (2.4) ish312!2D1f-13!(D2f+fyDf)+ (2.14)In order to let (2.5) represent the truncation error, the expansion of (2.5) must agree with (2.14). Accordingly1+2+3=0 (2.15)2D1f+2D2f=0 (2.16) 12!2D12f+12!3D22f=12!2D1f-13!D2f (2.17)31fyD1f=-13!fyDf (2.18)In order to let the equation above be satisfied independent of f, if we set =, 1=1+1 (2.19)Or D1=D, D2=1D (2.20)Then1+2=1 (2.21) 2=12 (2.22) 1+2+3=0 (2.23)2+31=0 (2.24) 1222+12312=1222-16 (2.25)31=-16(2.26)As these simultaneous equations are the system of six equations among eight unknown parameters, they have two degree of freedom. If we let 0,10, 1 , 23and solve these equations using and 1 as parameters, we obtain 1=2-12 2=12 1=3-261 2=3-26(-1) 1=1(-1)(3-2) 3=3-261(1-) (2.27)Any pair of and 1, chosen arbitrarily, correspond with a formula. Several examples are given below.(i) in case of =12,1=1 1=0 2=1 1=-16 2=13 3=-16 1=2So we get k1=hf(x0, y0) k2=hf(x0+12h, y0+12k1) k3=hf(x0+h, y0-k1+2k2) yx0+h=yx0+k2 T=-16k1+13k2-16k3 (2.28) (ii) in case of=1 1=12 1=2=12 We have 1=13 2=13 3=-23 1=14So we obtain k1=hf(x0, y0) k2=hf(x0+h, y0+k1) k3=hf(x0+12h, y0+14k1+14k2) yx0+h=yx0+12(k1+k2) T=13(k1+k2-2k3) (2.29)(iii) in case of =23 1=14 2=34 The second-order formula in this case was found by A. Ralston and has the minimum error bound 8. The substitution of the given value into the simultaneous equations of (2.23) to (2.26) yields the following equations;1+2+3=0 232+31=0 292+12312=0 2331=-16 (2.30)The solution of these equations leads to the following two formulas.k1=hf(x0, y0) k2=hf(x0+23h, y0+23k1) k3=hf(x0+23h, y0+1112k1-14k2) yx0+h=yx0+14k1+34k2 T=-k2+k3 (2.31)k1=hf(x0, y0) k2=hf(x0+23h, y0+23k1) k3=hf(x0, y0-14k1+14k2) yx0+h=yx0+14k1+34k2 T=k1-k3 (2.32)These formulas are able to evaluate the error at each step and reliable even when the function f makes sudden changes, but, are as a matter of fact, not so practical.3. Third-order FormulasThe general form of the formula is k1=hf(x0, y0) (3.1)k2=hf(x0+h, y0+k1) (3.2)k3=hf(x0+1h, y0+1k1+1k2) (3.3)k4=hf(x0+2h, y0+2k1+2k2+2k3) (3.4)y1=yx0+i=13iki (3.5)y2=yx0+i=14iki (3.6)T=y1-y2 (3.7)Where i, i, i, i, i and i are all scalars. (3.5) is a third-order formula while (3.6) is a fourth-order formula, and T denotes the estimated truncation error of (3.5). if we expand (3.5) into scalars:y1=yx0+i=13iki=yx0+1hf+2hf+hD1f+h22!D12f+h33!D13f+h44!D14f+3hf+hD2f+ h22!D22f+h33!D23f+h44!D24f+h21(fyD1f+h2!fyD12f+hD1fD2fy+h23!fyD13f+ h22D12fD2fy+h22!1fyD1f2+h22!D1fD22fy+=y0+1+2+3hf+h22D1f+3D2f+h32!2D12f+3D22f+h331fyD1f+h43!2D13f+3D23f+h42!31fyD12f+h331D2fyD1f+(3.8)On the conditions of (3.5) being a third-order formula, and also with (2.11) and (3.8), we get1+2+3=12D1f+3D2f+4D3f=12!Df 2D12f+3D22f+4D32f=2!3!D2f 31D1f=13!Df (3.9)On the other hand, for (2.11) and (3.8), we get the truncation error of (3.5) as follows:E1=h43!2D13f+3D23f-14D3f+h4fy2!31D12f-112D2f-h44!fy2Df+h431D2fyD1f-18DfDfy+(3.10)If we let D1=D, D2=1D and D3=2D, (3.9) is 1+2+3=1 (3.11) 2+31=12 (3.12) 22+312=13 (3.13)31=16(3.14)And (3.10) is expressed asE1=h43!D3f33+313-14+h4fy2!D2f321-112-h44!fy2Df+h4DfDfy311-18+ (3.15)Likewise, expanding (3.6) into series,y2=yx0+i=14iki=yx0+1+2+3+4hf+2D1f+3D2f+4D3fh2+h32!2D12f+3D22f+4D32f+h43!2D13f+3D23f+4D33f+h54!2D14f+3D24f+4D34f+h3fy31D1f+42D1f+2D2f+h4fy2!31D12f+42D12f+2D22f+h4fy2412D1f+h431D1fD2fy+42D1f+2D2fD2fy+(3.16)On condition of (3.6) being a fourth-order formula and with (2.11) and (3.16),1+2+3+4=1 (3.17)2D1f+3D2f+4D3f=12!Df (3.18)2D12f+3D22f+4D32f=12!D2f (3.19)31D1f+42D1f+2D2f=13!Df (3.20)2D13f+3D23f+4D33f=14D3f (3.21)31D12f+42D12f+2D22f=112D2f (3.22)412D1f=14!Df (3.23)31D1fD2fy+42D1f+2D2fD2fy=34!DfDfy (3.24)If we let D1=D, D2=1D and D3=2D as before, the question from (3.17) to (3.24) are1+2+3+4=1 (3.25)2+31+42=12 (3.26)22+312+422=13 (3.27)31+42+12=16 (3.28)23+313+423=13 (3.29)321+422+122=112 (3.30)412=124 (3.31)311+42+122=18 (3.32)When we obtain the solutions of equations from (3.11) to (3.1) and from (3.25) to (3.32), the formulas which are made using the coefficients of these solutions are the ones we desire. These are twelve simultaneous equations among thirteen parameters, so they have one degree of freedom. But as we need three degrees of freedom in order to determine freely the evaluating point of a third-order formula, it is quite doubtful whether the equations have solutions. So we test it as follows:When0, 1 , 2, 10, 23, 20, 12 3-41+610with (3.26), (3.27) and (3.29), we obtain2=3-41+2+621121-(2-) (3.33)3=-3-42+62121-2-1(3.34)4=-3-41+611221-2(2-)(3.35)With (3.25)1=1-(2+3+4) (3.36)While (3.11) to (3.14) 2=1(121-1/3)11-3=2-3611-1=611-312+21-2+32611-1=1(1-)2-3(3.37)And with (3.31), (3.35) and (3.37)2=2(2-1)(2-)(2-3 )211-3-41+61(3.38),With (3.38)2=121216-31-12(3.39)Using the equations from (3.33) to (3.38), the 2 of (3.39) is expressed by , 1 and 2. After the substation of these results into (3.30) and (3.32) and the collection of terms, both (3.30) and (3.32) are 3-46-9=12(3.40)Accordingly (3.30) and (3.32) are satisfied when 0. Hence the following conclusion is obtained. First we determine, 1 and 2 freely so as to let them satisfy a few conditions above stated. The system of equations from (3.11) to (3.14), from (3.25) to (3.29) and (3.31) is satisfied by the parameters which were computed using the equations from (3.33) to (3.39). so if is sufficiently small, the right-hand sides of (3.30) and (3.32) can be equated to their left-hand side with sufficient accuracy.Now, if we =160 1=12 and 2=1 , we obtain 1=10 2=-30029 3=3929 1=29039 1=16 2=0 3=23 4=16 2=-342251131 2=11758Then we get k1=hf(x0, y0) k2=hf(x0+160h, y0+160k1) k3=hf(x0+12h, y0-541178k1+29039k2) k4=hf(x0+h, y0+191832165598k1-342251131k2+11758k3) y1=y0+10k1-30029k2+3929k3 y2=y0+16k1+23k2+16k3 E1=y1-y2 (3.41)Where E is the estimated truncation error of y.Next, in order to increase the accuracy of the error estimation, we try to optimize the coefficients from the standpoint of truncation error. With this purpose, we take up that measures of the accuracy for truncation error which are independent if the given differential equations. When we set 1 and 2 arbitrarily and compute each parameter taking the above procedure, we obtain the truncation error of (3.6) by the use of (2.11) and (3.16) as follows:E2=h4a1fyD2f2+a2DfDfy+h5b1D4f+b2DfD2fy+b3fD3fy+b4DfD2fy+b5fD2fy2+b6Df2fyy+b7f(fyD)2+b8Dffy3+(3.42)Where ai, i=1, 2, bi, j=1, 2,-,8 are the functions of 1 and 2 and their functional parameters. So they are expressed as below:a1=321+422+122-112a2=311+42+122-18 b1=14!(24+314+424-15)b2=12!3121+42+1222-110 b3=13!331+432+13222-120 b4=12!3211+422+1222-115 b5=12!(4212-160) b6=12!3212+42+122-120 b7=41+212-7120 b8=-1120 (3.43)In order to search for the combination of coefficients which minimize the terms of the truncation error of order h4 and h5 as well as possible, the following are used the measures of the accuracy for truncation error above stated；A=i=12|ai|(3.44)B=i=18|bi| (3.45)C=i=18bi2 (3.46)D=16b1+4b1+b2+3b3+2b2+3b3+b2+b3+8b4+b3+b5+2b5+b7+b5+b6+b7+b6+2b6+b7+b7+2|b8|(3.47)Among these, D was used as the measures of truncation accuracy by A. Ralston to optimize Runge-Kutta formula, and so were A, B, and C by T. E. Hull and R. L. Johnson 9.When we set the coefficients in the equations from (3.1) to (3.7) as follows:=0.001 =0.001 1=0.7 1=-244.3175262 1=225.0175262 2= 0.8 1= 136.1510201 1=-136.0025668 2=9.6515466956 1=-23.52380952 2=23.84358067 3=0.6802234484 1=-53.31547619 2=53.71521268 3= 0.3392601675 4=0.2610033375 We can obtain one of the formulas we desire. In this case, A=6.2510-5, B=2.7010-2, C=1.6910-4 and D=7.5210-2. A. Ralston studied the values of D concerning to the representative fourth-order method 8. From these values, it is clear that (3.6) has the accuracy of fourth order.Table 1y1=y1-yActual ErrorT=y1-y2Error EstimateT/1.6093414971-0.0011685030-0.00104194900.892The table 1 shows the first step in the numerical solution of the differential equationdy/dx = 5y/(1+x)with initial condition y0=1 when integrated at pitch 0.1, and also shows its true error, its evaluated error and the ratio of the evaluates error to the true one. The computation was made at 10 dicimal digits using a digital computer, FACOM 231.4. Kutta-Merson ProcessThe general expression of the formula of this type is shown as yn+1=yn+i=15iki (4.2)k1=hf(xn, yn) (4.2)k2=hf(xn+h, yn+k1) (4.3)k3=hf(xn+1h, yn+1k1+1k2) (4.4)k4=hf(xn+2h, yn+2k1+2k2+2k3) (4.5)k5=hf(xn+3h, yn+3k1+3k2+3k3+3k4) (4.6)The e stimation formula of the truncation error n is n=i=15iki(4.1)Where i, i, i, i, i, i and i are scalars. In Kutta-Merson process each coefficient in the above formula is set as follows: 1=16 2= 3=0 4=23 5=16 =13 1=13 2=12 3 =1 = 13 1=16 2=18 3=12 1 =16 2=0 3=0 2=38 3=-32 3=2 1= 115 2=0 3=-310 4=415 5=-130 This process, which enables us to estimate the error by adding only one more computation of function to the ordinary fourth-order methods, may look rather attractive when compares with the complicatedness of the former method of estimation in which we had to integrate twice changing the pitch each time. But, in spite of the misunderstanding of its deviser and the introducers who took this process for a fourth-order method with the error-estimating ability, it is no more than a third-order formula in general cases where the assumption of linearity does not hold, though it had this ability. It was because Merson approximated the arbitrary function by a linear expression in (x, y) in order to derive the formula easily. This fact will be demonstrated in section 6 by the result of the experiment and by the measures of the accuracy for truncation error.Further, F. Ceshino also made a researching giving the ability of error evaluation to a third-order method by means of adding two more computations to the ordinary formula 6.5. Runge-Kutta Formulas Which uses Five Function ValuesIn this section, in preparation for making fourth-order methods with the ability of error estimation, we study the optimization of the Runge-Kutta formula that uses five function values. It is impossible to make fifth-order method in a strict sense if we use five function values. But, in fact, it is possible to obtain the accuracy of fifth-order method by means of optimizing the coefficients-that is ,by means of making each term of order h5 of the truncation error smaller than that of orderh6 . The general expression of the formula isk1=hf(x0, y0) k2=hf(x0+h, y0+k1) k3=hf(x0+1h, y0+1k1+1k2) k4=hf(x0+2h, y0+2k1+2k2+2k3)k5=hf(x0+3h, y0+3k1+3k2+3k3+3k4) yx0+h-yx0=i=15iki (5.1)The question is how to determine i, i, i, i, i, i and i. So we set them to let Taylors expansion of (5.1) agree with (2.11) up to terms of order h5 as well as possible. In order to obtain Taylors expansion at the point (x0, y0) of (5.1), we define the operaters D1, D2, D3 and D4 as follow:For the expansion of k2 D1=x+f0yFor the expansion of k3 D2=1x+(1+1)f0yFor the expansion of k4 D3=2x+(2+2+2)f0yFor the expansion of k5 D4=3x+3+3+3+3f0y (5.2)Each of ki in (5.1) is expanded in turn. Firstly we let k1=hf0 (5.3)Next, in order to obtain the expansion of k2, we let h1=hx+hf0y=hx+k1yD11So we get k2=hf+D11f+D112f2!+D113f3!+D114f4!+=hf+hD1f+h22!D12f+h33!D13f+h44!D14f+ (5.4)Thirdly, to obtain the expansion of k3, we let1hx+1k1+1k2y=h(1x+1f0y)+1k2y=hD2+1(k2-hf0)yAnd using (5.4), we get =hD2+h21D1f+h2!D12f+h23!D13f+yD21Therefore,k3=hf+D21f+D212f2!+D213f3!+D214f4!+=hf+hD2f+h22!D22f+h33!D23f+h44!D24f+h21fyD1f+h2!fyD12f+hD1fD2fy+h23!fyD13f+h22D12fD2fy+h22!1fyyD1f2+h22!D1fD22fy+ (5.5)Fourthly, to obtain the expansion of k4, we let 1hx+2k1+2k2+2k3y=hD3+2k2-hf0+2k3-hf0yD31If we use (5.4) and (5.5), hD3+h22D1f+h2!D12f+h23!D13f+2D2f+h2!D22f+21D1ffy+h23!D23f+h22!D23f+h22!1fyD12f+h21D1fD2fy+yD31Thereforek4=hf+D31f+D312f2!+D313f3!+D314f4!+=hf+hD3f+h2fy2D1f+h2!D12f+h23!D13f+2D2f+h1D1ffy+h2!D22f+h22!1fyD12f+h21D1fD2fy+h23!D23f+12!(h2D32f+2h3D3fy2D