Application of explicit methods with extended stability regions for solving stiff problems Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 519.622

Application of Explicit Methods with Extended Stability Regions for Solving Stiff Problems

Eugeny A. Novikov*

Institute of Computational Modeling SB RAS Akademgorodok, 50/44, Krasnoyarsk, 660036

Russia

Mikhail V. RybkoV

Institute of Mathematics and Computer Science Siberian Federal University Svobodny, 79, Krasnoyarsk, 660041

Russia

Received 10.11.2015, received in revised form 15.01.2016, accepted 20.02.2016 An algorithm is developed to determine coefficients of the stability polynomials such that the explicit Runge-Kutta methods have a predetermined shape and size of the stability region. Inequalities for accuracy and stability control are obtained. The impact of the stability control on efficiency of explicit methods to solving stiff problems is shown. Numerical calculations confirm that the three-step method of the first order with extended stability region is more efficient than the traditional three-stage method of the third order.

Keywords: stiff problem, explicit methods, stability region, accuracy and stability control. DOI: 10.17516/1997-1397-2016-9-2-209-219.

Introduction

In some cases large-scale stiff problems need to be solved with the algorithms based on explicit methods as shown in [1, 2]. The L-stable methods suffer from the decomposition of the Jacobi matrix [3,4] and in case of high dimension of a problem the calculation of the inverse to the Jacobi matrix defines overall computational efforts. At the same time, the algorithms based on explicit formulas are more efficient if stiffness of the problem allows to get an approximation to a solution in a reasonable time [5].

At the present time, algorithms of variable step and order are developed [1,2,6]. They use inequalities for stability control as a criterion for choosing between methods of high and low order of accuracy at the integration step. The numerical formulas of low-order methods are based on the same stages as high-order methods but its stability regions are much wider then stability regions of high-order ones. In a settling region there is no point in using high-order methods because the integration step is restricted by the condition of stability. The efficiency can be increased by applying some low-order methods with extended stability regions there. Further raising of efficiency may be achieved by developing not only variable order and step algorithms but algorithms with variable number of stages. We need some pool of low-order methods with extended stability regions to obtain them [1, 7].

The stability polynomials of degree up to m = 13 are constructed in [1,8]. Here an algorithm to determine the coefficients of the stability polynomial coefficients is developed such that the

* [email protected] 1 [email protected] © Siberian Federal University. All rights reserved

corresponding explicit Runge-Kutta methods have a predetermined shape and size of the stability region. Formulas for the coefficients of Runge-Kutta methods with extended stability regions are obtained. The stability and accuracy control inequalities for the first order methods are constructed. The three-stage explicit Runge-Kutta method is considered as an example. The results of comparison of the three-step method of the first order with extended stability region and the traditional three-stage method of the third order applied to solving stiff problems are given.

1. Explicit Runge-Kutta methods

We consider the Cauchy problem for the stiff system of ordinary differential equations

y = f (t,y), y(to) = yo, to < t < tk, (1)

where y and f are real W-dimensional vector functions, t is an independent variable. In [9] the authors propose to solve (1) with explicit Runge-Kutta methods

m i — 1

yn+1 = yn + Piki, ki = hf(tn + aih,yn + fa j, (2)

i=1 j = 1

where ki, 1 < i < m are the stages of the method, h is an integration step, pi, ai, Pij, 1 < i < m, 1 ^ j ^ i — 1 are numerical coefficients that define stability and accuracy characteristics of the scheme (2). For the sake of simplicity further we consider the Cauchy problem for the autonomous system of ordinary differential equations

y' = f (y), y(to) = yo, to < t < tk■ (3)

For solving (3) we also may write formulas (2) in the following form:

m i—1

yn+1 = yn + Piki, ki = hf(yn + Pj j ■ (4)

i=1 j = 1

The results given below can be used for non-autonomous systems if in (2) we assume

i—1

a1 =0, ai ^^ fy, 2 ^ i ^ m■ j=1

2. Order conditions

Below we need matrix the Bm with elements bj of the form [1,3]

&1i = 1, ^ i ^ m, bki =0, 2 < k < m, 1 < i < k — 1,

i—1 (5)

bki /3ijbk—\j, 2 ^ k ^ m, k ^ i ^ m,

j=k—1

where Pij are numerical coefficients of the scheme (2) or (4). The stability of one-step methods is usually investigated by applying a Runge-Kutta method to a linear scalar equation known as the Dahlquist equation [10,11]

y' = Xy, y(0)= yo, t > 0

with complex A, K(A)< 0. The variable A is considered as a certain eigenvalue of the Jacobi matrix of the problem (1) or (3). Applying numerical scheme (4) to solve the Dahlquist equation

we get

m

Vn+1 = Qm (z)yn, Qm(z) = 1+^ CiZ%,

m i=1 (6)

ci ^^ bij pj, 1 ^ i ^ m, j=i

where z = hA. Denoting c = ( c1, ..,cm)T and p ={p1, ..,pm)T, the latter equality (6) we can rewrite in the form

BmP = c, (7)

where the elements of the matrix Bm are defined in (5).

Expanding exact and approximate solutions in Taylor series in degrees of h, we get

y(in+i)= y(tn)+hf + 2h2 f'f + 1 h3 [f ''f2 + f'2f]+O(h4),

yn+1 = yn + b1jPj) hfn + b2jPj) h2ffn +

j=1 j=2

m m

+ (E b3jPj) h3fnfn + 0.5 fc b2j Pj) h3fnfn + O(h4),

j=3 j=2

(8)

where elementary differentials are calculated on the exact y(tn) and the approximate yn solutions, respectively. If we compare equalities (8) assuming that yn = y(tn) we can see that the numerical scheme (4) has the first order of accuracy if

m

Yjb1jPj = 1

j=1

holds. Hence, we need to put c1 = 1 in the linear algebraic system (7) to construct m-stage methods of the first order. The remaining coefficients ci, 2 ^ i ^ m, can be chosen according to the required stability characteristics. The stability function of an m-stage explicit Runge-Kutta method is known to be a polynomial of degree m. So we need the coefficients of this polynomial in order to use (7). Conversely, the coefficients of the polynomial can be calculated using (7) if the coefficients of the method (4) are given.

3. Stability polynomials

Let two integer number k and m, k ^ m be given. Consider the polynomial

k m

Qm,k (s) = 1 + 53 ciXi + ^ ciXi, (9)

i=1 i=k+1

where the coefficients ci, 1 ^ i ^ k, are given and ci, k +1 ^ i ^ m, are free. The coefficients ci, 1 ^ i ^ k, are usually defined from the approximation requirements. Therefore, to be precise, we assume ci = 1/i!, 1 < i < k below.

Let |7m| denote the length of the stability interval of an m-stage explicit Runge-Kutta method, i.e. on the interval [ym, 0 the condition lQm,k (x)1 < 1 holds. It is easy to show that the more m, the faster coefficients of the stability polynomial considered on the interval

[im, 0 tend to zero [1]. Using the algorithm [1, 8] the coefficients ci, k + 1 < i < m, are obtained up to the degree m = 13. Therefore, we consider an algorithm of obtaining the polynomials with prescribed properties on the interval [—1,1]. In this case the coefficients ci grow not so fast and it is possible to construct polynomials of degree m > 13. We transform the interval [lm, 0 to the interval [—1,1] with the change of variables x =1 — 2z/^m and get the polynomial

m

Qm(z)=J2 diZi. (10)

i=0

Then, the coefficients di, 0 ^ i ^ m, of the polynomial (10) are related to the coefficients ci,

0 ^ i ^ m, of the polynomial (9) by

c = UVd, (11)

where d = (d0, ..,dm)T, c = (c0, ..,cm)T, U is a diagonal matrix with uii = (—2/Ym)% \

1 ^ i ^ m + 1, on the main diagonal and the elements vij of the matrix V are given by the formulas

vlj = 1, 1 < j < m + 1, vij = vij-1 + vi-1'j-1, 2 < i < j < m + 1, vij =0, i > j.

It is easy to see that V is the Pascal triangle and its elements can be easily calculated by a recurrent formula. Having constructed the polynomial (10) on the interval [—1,1, we can calculate the coefficients of the polynomial (9) on the interval \_Ym, 0 using (11).

Now, turn to the polynomial (10). Let z1, z2,.., zm-1 denote the points of extremum of (10), assuming that z1 > z2 > ... > zm-1. We can find the coefficients di, 0 < i < m, such that the polynomial (10) assumes some prescribed values at the points of extremum zi, 1 < i < m — 1, i.e.

Qm(zi) = Fi, 1 < i < m — 1,

where F(z) is a prescribed function Fi = F(z^j. For this purpose consider the algebraic system of equations in variables zi, 1 ^ i ^ m — 1, and di, 0 ^ i ^ m, of the from

Qm(zi) = Fi, Q'm(zi) =0, 1 < i < m — 1;

1

idiz

i=1

Qm(z) = E idizi-1, (12)

where the normalization requirements Qm( — 1) = ( — 1)m and Qm(1) = 1 hold.

We rewrite (12) in the form that is convenient for computations. Let y, w and g denote vectors with components

yj = zj, 1 ^ j ^ m — 1; wi = di-1, 1 ^ i ^ m +1;

gi = Fi, 1 < i < m — 1; gi = 1, i = m; gi = ( — 1)m, i = m + 1.

Let E1 and E2 denote matrices of dimensions (m +1) x (m + 1) and (m — 1) x (m + 1), respectively, with the following elements on the main diagonal

j = j — 1, 1 ^ j ^ m +1; = 1/yi, 1 ^ i ^ m — 1.

Let A denote an (m + 1) x (m + 1) matrix with the following elements

aij = yj 1, 1 ^ i ^ m — 1, 1 ^ j ^ m +1;

amj = 1, am+1'j = (-1)j+1, 1 < j < m + 1. Then, we can rewrite the problem (12) in the form

Aw - g = 0, E2AE1w = 0. (13)

We use the relaxation method for numerical solution of (13) [1,12]. After obtaining the coefficients of the polynomial (10) we can calculate the coefficients of the polynomial (9) using (11). The value of jm is found from the condition that the desired polynomial corresponds to a first order method, i.e. c1 = 1. If we take the second row from (13) and make necessary computations, we get

2 m+1

Ym = - — • E V2j dj .

j = 1

4. Construction of stability regions

Numerical results [1,8] show that the more m, the less the coefficient cm of the polynomial (9) and, in particular, if m =13 and k =1 the value of cm is about 10~26. Rounding errors make difficulties for computation of polynomial coefficients when m > 13. When solving the system (13), the coefficients di, 0 < i < m, of the polynomial (10) are established to grow in absolute value as m increases and if m = 13 the maximum of |di| is about 105. When m = 25 we have the maximum of |cii\ to be about 109, i.e. the coefficients d grow with lower speed. We can get the coefficients of the polynomial (9) from the coefficients of the polynomial (10) by formula (11) after solving the problem (13). It allows to calculate coefficients of stability polynomials up to m = 27.

The size and shape of stability region depend on the values of the function F. If we assume Fi = ( — 1)i, k ^ i ^ m +1, then the length of the stability interval is known to be |Ym| = 2m2 [1]. In this case we have the maximal length of stability interval along the real axis for given m. Such stability regions are almost multiply connected and rounding errors can provoke appearing of small imaginary parts of eigenvalues of the Jacobi matrix so that stability region shrinks.

In order for rounding errors to not reduce the stability region, we need to "stretch" it along the imaginary axis in the points of extremum of the stability polynomial. For this purpose we assume Fi = ( —1)i • ¡j, 1 < i < m — 1, 0 < j < 1. Calculations show that if we choose j = 0.9, then the length of the stability interval reduces by 5-8% in comparison with the maximal possible length that equals 2m2. In addition, such stability region stretches along the imaginary axis in the points of extremum. It provides better properties of stability of the method to rounding errors while stability region is reduced insignificantly. In case of j = 0.95 the length of the stability interval reduces only by 3-4%. If j = 0.9, the stability region of the five-stage Runge-Kutta method has the shape shown in Fig. 1. The length of its stability interval is equal to |Ym| = 46.79. For better visualization of the roots of polynomial (9) in complex plane all the figures include level lines |Qm,k(x)| = q, if q is equal to 1, 0.8, 0.6, 0.4, and 0.2.

The less j from one to zero, the closer to each other the polynomial roots are situated on the real axis. So it is natural that the length of the stability interval of corresponding method decreases. Ellipses that are strongly pronounced at j =1 approach each other. Besides, it does not provide sufficiently large stretching along the imaginary axis. Therefore, depending on the characteristics of the problem to be solved, it is appropriate to use values of j in the interval between 0.8 and 0.95.

In case of the first-order methods, i.e. if k = 1, this condition can be met by choosing the values of the function F. For instance, if m = 4, k = 1, F = {0.85,0.95,0.85} we get the polynomial that meets this condition. Since m is even and all the values of F are positive, the graph of the polynomial does not cross the real axis, and the polynomial has two couples of

complex conjugate roots. That is why the stability region stretches along the imaginary axis in the complex plane and catches a part of it. At the same time the length of the stability interval along the real axis is not so big and equals |Ym| = 2.18. The less the values of F, the more the length of the stability interval is, and, in particular, if F = {0.55,0.65,0.55} the stability region becomes practically rectangular with |Ym| = 5.30. Further decreasing of the values of function F provides extension of the length of stability interval |Ym|, but stability region catches a smaller part of the imaginary axis. Thus, when constructing the first-order methods for solving problems with oscillating solutions, it is reasonable to choose stability polynomials that have a couple of complex conjugate roots near the origin in the complex plane {ft-A}. Besides, the values of function F that correspond to these roots should be chosen close to 1 in order that the stability region catches the maximal part of the imaginary axis.

If m is odd then at least one root of the polynomial (9) is real. We can use this in different ways. For example, if we place a real root of the polynomial of degree 5 between two couples of complex conjugate roots by choosing respective values of function F, then we can construct a stability region that is similar to rectangle. If the stability polynomial assumes the values F = {0.2,0.5, —0.5, —0.2} at the points of extremum, then the stability region of the respective 5-stage method has the shape shown in Fig. 2. In the general case, we can construct the stability regions of different shape and size by choosing parameters m, k and the values of function F.

Fig. 2. The stability region for m = 5, k = 1, F = {0.2, 0.5, —0.5, —0.2}, Ym| = 17.21

5. Accuracy and stability control

Now we derive the inequality for accuracy and stability control of the first-order explicit Runge-Kutta methods. These schemes are supposed to be used in a settling region where an integration step is restricted by stability but not by accuracy of the method. Accuracy control is subsidiary there and to obtain the respective inequality we use the estimate of the local truncation error. Comparing Taylor expansions (8) of the exact solution y(tn+1) and the approximate solution yn+1 up to terms with h2, provided that yn = y{tn), we get that the local truncation error of the scheme (2) has the following form:

m

Sn,1 = (0.5 — J2 b2iPi) h2f 'f + O(h3) .

i=2

Using (5), it can be rewritten in the form

Sn,1 =0.5(1 — 2c2)h2f 'f + O(h3),

where c2 is the coefficient at x2 of the stability polynomial (9).

The value of Sn1 can be estimated using stages ki, 1 < i < m, that have already been calculated in many ways. We use the following notation:

A'n = gm1\\k2 — kl||, An = g'm1\\hf(yn+1) — kl|, where | | | |

= |1 — 2cm2\/2, g'm1 = g'm1/\321\.

Then, taking into account that

k2 — k1 = ¡321k2 ffn + O(h3),

hf (yn+1) — k1 = h2fn fn + O(h3),

we can use the following inequalities to control accuracy and choose the integration step, respectively

An < e,A<n < e. (14)

The vector k1 depends linearly on the step size. So using the first inequality of (14) the repeated computation of the solution is accompanied with just one extra calculation of the right hand side of the system of differential equations (1). If the integration step is successful, then the second inequality does not provide the growth of computational efforts because the vector f(yn+0 is used only at the next integration step. If the second inequality of (14) does not hold, then the return will be expensive in terms of computational efforts. Besides, the more m, the higher computational efforts. Nevertheless, preliminary control of An usually allows to avoid extra computations.

Now we derive the inequality for stability control by the method described in [6]. We write the stages k1 , k2 and k3 for the problem y = Ay, where A is a matrix with constant coefficients. As a result, we get ( )

k1 = X • yn, k2 = (X + 321X2) • yn,

k3 = [X + (331 + 332)X2 + 321332 X3] • yn, where X = hA. We find the coefficients s^ s2, and s3 assuming that the following equality

sk + S2k2 + S3k3 = X3 • yn

holds. It is easy to see that this requirement holds on the conditions that

031 + 032 — 021 031 + 032 1

«1 =

ß2lß32

«2 =

ß2lß32

S3 =

ß2lß32

It is also easy to see that the equality

1

— {k2 — ki) = X2 • yn.

021

holds. Then, according to [6] the maximal eigenvalue vn = hAn max of the Jacobi matrix of the system (1) can be estimated by the power iterations. I.e. we use the following inequality for stability control of the first-order methods (2):

vn < \lm,i\, Vn = ß2i • max

1 ' 1 1<i<N

«ikl + S2k2 + S3^i3

k2 — k1

(15)

where the positive constants Ymi depend on the size of the stability regions of numerical schemes. Inequality (15) can be used on each step for choosing efficient numerical scheme.

6. First-order method

The problem of construction of first-order explicit Runge-Kutta methods with given stability regions comes to solving the system of linear algebraic equations (7) with nondegenerate matrix Bm. In this case the values of the vector c = (1,c2, ...,cm) define the size and the shape of the stability region. As an example, we consider 3-stage explicit Runge-Kutta method in the form

yn+1 = yn + Pikl + P2k2 + P3k3,

(16)

ki = hf{yn), k2 = hfiyyn + 02k), k3 = hfiyyn + 03iki + 032k^ .

Assume that 021 = 0.5, 031 = —1 and [332 = 2. Then, on each step the increments k1, k2 and k3 are calculated at the points tn, tn + 0.5h and tn + h, respectively. The numerical results show that calculation of stages at these points increase the reliability of calculations when integrating stiff problems.

The coefficients p1, p2 and p3 are obtained by means of equality (7), where the vector of the right part of the system includes coefficients of the stability polynomial of the method. These are calculated by solving the problems (11) and (13). If we assume Fi = 0.95 • ( — 1)®, 1 < i < 2, then the length of the stability interval of the method decreases insignificantly in comparison with the maximal possible length 2m2 that equals 17.46. The values of the coefficients of the stability polynomial are c1 = 1, c2 = 0.15209292726978 and c3 = 0.00580524400854 for the given F®. And the stability region is simply connected in this case. The requirement of the first order leads to the equalities

Pi + P2 + P3 = 1, 02iP2 + (03i + ^32)P3 = c2, 021032P3 = c3. As a result, we have coefficients

c3 c2 — (031 + 032) P3

P3 = 0210322 , P2 = -021-, P1 = 1 — P — P3.

We use inequalities (14) for the accuracy control of the numerical schemes of the first order. The length of the stability interval of the first-order method equals 17.46. So we can use inequality v3 ^ 17.46 where v3 is defined in formula (15) for the stability control of this scheme.

7. Third-order method

For numerical solution of (3) we consider an explicit three-stage Runge-Kutta method in the form

yn+1 = yn + r1 k1 + T2k2 + T3k3,

k1 = hf(yn), k2 = hf(yn + 321k1), (17)

k3 = hf(yn + 33k + 332k2),

where the coefficients 321, 331 and 332 are defined in the description of the first-order method and the coefficients T1 , T2 and T3 need to be defined. The conditions of the third order of accuracy are the following

1 11 1 1

T1 + T2 + T3 = 1, 2 T2 + T3 = 2, 4 T2 + T3 = 3, 2T3 = 3.

As a result, we have the coefficients of the scheme (17) of the third order

1 12 1

321 = X , 331 = —1, 332 = 2, T1 = -, T2 = -, T3 = -. 2 6 3 6

We derive the accuracy control inequality for a third-order method using the idea of embedded methods. We consider an auxiliary scheme yn+\y1 = yn+k2 that has the second order of accuracy. Then, the accuracy control inequality is the following [6]

\\yn+1 — yn+1,1 \\ = \\k1 — 2k2 + k3\/6 < e,

where || • || is some norm in RN, e is the required accuracy of calculations.

8. Numerical results

The calculations were performed on AMD A6-3420M APU with double precision. The required accuracy was e = 10~2. In concrete computations the norm ||£n|| from the inequalities for the accuracy control was calculated by the formula

UU = max |e|/(|yn|+ t),

where i is the number of a vector component, t is a positive parameter. If the inequality \yln\ <t holds on the component with the number i, then the absolute error T • e is controlled, otherwise we control the relative error e. During calculations the parameter p was chosen in a way that the practical accuracy was not worse than the required one.

We compared the efficiency of two methods. The first one is the Runge-Kutta method of the first order with accuracy and stability control described in section 6. And the second one is the traditional three-stage Runge-Kutta method of the third order with accuracy control. Below numbers is, iw and if denote the total number of integration steps, the total number of repeated solution calculations (returns) because of required accuracy violation and the total number of calculations of the right hand side of the problem (1), respectively.

As a test example, we consider simplified model with periodic solution that corresponds to a Belousov-Zhabotinsky reaction (oregonator) in the form [14]

y1 = 77.27(y2 — yy + y1 — 8.375 • 10~6y\), y2 = (—y2 — yxy2 + y3)/77.27, y'3 = 0.161 y — y3), t G [0, 300], y1 (0) = y3(0) = 4, y2(0) = 1.1, h0 = 10-3.

The given problem is "too" stiff for explicit methods. Nevertheless, we present exactly this example in order to emphasize advantages of explicit methods with stability control.

Using the first-order method the computational efforts are is = 427 018, iw = 16 149 and if = 1 725 219. For the third-order method respective computational efforts are is = 2 903 722, iw = 769 200 are if = 10 249 566. The results show fivefold benefit of the first-order method with extended stability region in comparison with the traditional method of the third order. Both decreasing of number of returns because of stability control and extending the stability region lead to the growth of efficiency. At the end of the integration interval the practical accuracy is better than the required one for all the algorithms. This tendency holds when solving other examples [14].

Conclusion

The coefficients of stability polynomials of degree up to m = 27 and the respective methods of the first order are obtained using the algorithm described above. The maximal length of the stability interval of m-stage first-order Runge-Kutta method is known to be 2m2. So each computation of the right part of the system corresponds to 2m units of the stability interval. This means that if the integration step is restricted by stability of the scheme then the more m, the more efficient method becomes. Using the inequality for stability control does not lead to increase of computational efforts because the estimation of the maximal eigenvalue of the Jacobi matrix of the system (1) is calculated via stages computed in advance and it does not lead to the growth of number of calculations of function f. Thus, the first-order methods with extended stability regions allow to significantly increase the efficiency in a settling region where the step is restricted by stability.

The developed set of the first-order methods with stability and accuracy control can be applied in the algorithms with variable step, variable order and variable number of stages. In this case the efficiency increases not only because of switching from the method of high order into the method of low order in a settling region, but by means of "boosting" an integration step in this region using the first-order methods with extended stability regions.

This work is partially supported by Russian Foundation of Basic Research (grant 14-0100047).

References

[1] E.A Novikov, Explicit methods for stiff systems, Novosibirsk, Nauka, 1997 (in Russian).

[2] E.A.Novikov, Yu.V.Shornikov, Computer simulation of stiff hybrid systems, Novosibirsk, Izdat. NSTU, 2012 (in Russian).

[3] E.Hairer, S.P.Norsett, G.Wanner, Solving ordinary differential equations. I. Nonstiff problems, Berlin, Springer, 1987.

[4] E.Hairer, G.Wanner, Solving ordinary differential equations. II. Stiff and differential-algebraic problems, Berlin, Springer, 1996.

[5] V.A.Novikov, E.A.Novikov, Raising the efficiency of algorithms for the integration of ordinary differential equations at the expense of loss of stability, USSR Computational Mathematics and Mathematical Physics, 25(1985), no. 4, 39-43.

[6] E.A.Novikov, Variable order and step algorithm based on a stages of Runge-Kutta method of third order of accuracy, Izv. Saratov Univ. (N.S.), Ser. Mat. Meh. Inform., 11(2011), no. 3(1), 46-53 (in Russian).

[7] E.A.Novikov, M.V Rybkov, The numerical algorithm of constructing of stability regions for explicit methods, Sistemy upravleniya i informatsionye tehnologii, 55(2014), no. 1.1, 173177 (in Russian).

[8] E.A.Novikov, The construction of the stability regions for the explicit Runge-Kutta methods, Vychislitelnye metody i programirovanie, 10(2009), 248-257 (in Russian).

[9] V.A.Novikov, E.A.Novikov, Control of the stability of explicit one-step methods of integrating ordinary differential equations, Soviet Math. Dokl., 30(1984), №1, 211-215.

[10] G.Dahlquist, A special stability problem for linear multistep methods BIT, 3(1963), 23-43.

[11] K.Dekker, J.G.Verwer, Stability of Runge-Kutta methods for stiff nonlinear differential equations, Amsterdam, North Holland, 1984.

[12] G.Korn, T.Korn, Mathematical handbook, Moscow, Nauka, 1970 (in Russian).

[13] A.E.Novikov, E.A.Novikov, L-stable (2,1)-method for stiff nonautonomius problem solving, Vychislitelnye texnologii, Tom. 13, Vestnik KazNU, 58(2008), no. 3, 477-482 (in Russian).

[14] W.H.Enright, T.E.Hull, Comparing numerical methods for the solutions of stiff systems of ODE's, BIT, 15(1975), 10-48.

Применение явных методов с расширенными областями устойчивости для решения жестких задач

Евгений А. Новиков

Институт вычислительного моделирования СО РАН Академгородок, 50/44, Красноярск, 660036

Россия

Михаил В. Рыбков

Институт математики и фундаментальной информатики Сибирский федеральный университет Свободный, 79, Красноярск, 660041

Россия

Разработан алгоритм определения коэффициентов полиномов устойчивости, при которых явные методы типа Рунге-Кутта имеют заданные форму и размер области устойчивости. Получены неравенства для контроля точности и устойчивости численных формул первого порядка. Исследовано влияние контроля устойчивости на эффективность явных методов применительно к решению жестких задач. Приведены результаты расчетов, показывающие повышение эффективности трехстадийного метода первого порядка с расширенной областью устойчивости в сравнении с традиционным трехстадийным методом третьего порядка.

Ключевые слова: жесткая задача, явные методы, область устойчивости, контроль точности и устойчивости.