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

DOI: 10.17516/1997-1397-2020-13-2-242-252 УДК 519.622

First-Order Methods With Extended Stability Regions for Solving Electric Circuit Problems

Mikhail V. Rybkov* Lyudmila V. Knaub Danil V. Khorov*

Siberian Federal University Krasnoyarsk, Russian Federation

Received 16.01.2020, received in revised form 06.02.2020, accepted 25.03.2020 Abstract. Stability control of Runge-Kutta numerical schemes is studied to increase efficiency of integrating stiff problems. The implementation of the algorithm to determine coefficients of stability polynomials with the use of the GMP library is presented. Shape and size of the stability region of a method can be preassigned using proposed algorithm. Sets of first-order methods with extended stability domains are built. The results of electrical circuits simulation show the increase of the efficiency of the constructed first-order methods in comparison with methods of higher order.

Keywords: stiff problem, explicit methods, stability region, accuracy and stability control.

Citation: M.V.Rybkov, L.V.Knaub, D.V.Khorov, First-Order Methods With Extended Stability Regions for Solving Electric Circuit Problems, J. Sib. Fed. Univ. Math. Phys., 2020, 13(2), 242-252. DOI: 10.17516/1997-1397-2020-13-2-242-252.

Introduction

Systems of ordinary differential equations (ODEs) describe various dynamic processes in chemistry, physics, etc. One of the areas where ODEs may be effectively applied is the electric circuit theory. Any changes in electric circuit lead to transient processes where some voltage swells, electromagnetic oscillations, extra currents may occur. They can damage electrical devices. At the same time transient processes occur in electrical generators and other electric circuits. Many electric circuits problems are described by stiff systems of ODEs.

In some cases explicit methods are required to solve initial value problems of stiff ODEs because L-stable methods involve inversion of the Jacobi matrix of a system. This defines overall computational costs [1-2]. At the same time explicit methods do not require the Jacobi matrix computation and they are more preferable to use for problems with not so high stiffness ratio.

At present time various explicit and implicit methods were developed [3]. The former are used on transition regions where the integration step is restricted by the accuracy criterion and there is no requirements for large stability interval. The latter are for the regions where large stability interval gives an opportunity to pass the integration interval in "several steps". Nevertheless these

* mixailrybkov@yandex.ru https://orcid.org/0000-0002-6560-1435

t lvknaub@yandex.ru https://orcid.org/0000-0003-4857-2078

idanilkhorov@gmail.com https://orcid.org/0000-0001-8967-8341 © Siberian Federal University. All rights reserved

algorithms are not so efficient to solve high dimension systems of ODEs because of mentioned above reasons.

Variable order algorithms based on explicit schemes were developed [4]. They are applied to the regions where there is no need to use high-order methods. Low computational cost can be achieved by using there low-order methods with extended stability intervals which in fact play part of implicit methods from the point of view of the length of stability interval.

Low-order methods with large stability interval are needed to develop such algorithms. In addition the greater number of stages of a method (and therefore the higher the degree m of stability polynomial ), the large the stability interval is. The stability polynomials of degree up to m = 13 were constructed [2]. The algorithm to determine the stability polynomial coefficients was developed such that the corresponding explicit Runge-Kutta methods have a predetermined shape and size of the stability region [7].

Here implementation of the algorithm to obtain the stability polynomial coefficients with the use of the library for arbitrary precision arithmetic GMP is presented. Set of the firstorder methods with extended stability intervals is developed. Numerical simulation of Van der Pol oscillator shows that proposed algorithms are more efficient in comparison with the Merson method of fourth order of accuracy.

1. Explicit Runge-Kutta methods

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

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

where y h f are real W-dimantional vector functions, t is independent variable. To solve (1) the following explicit Runge-Kutta methods were proposed [2]

m i — 1

yn+1 = yn + Pmih, ki = hf(tn + aih,yn + 3ij j , (2)

i=1 j=1

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

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

To solve (3) we can also write formulas (2) in the following form:

i m

yni = yn + Pi+ij kj, 1 ^ i ^ m — 1, yn+1 = yn + Pmiki, (4)

j=1 i=1

where ki = hf(yrhi—1^J, 1 < i < m, yn,0 = yn. The results given below can be used for non-autonomous systems if we assume in (2) that

i—1

ai =0, ai ^^ ¡3ij, 2 ^ i ^ m. (5)

j=1

Below we need matrix Bm with elements bij in the form [2]

< i < k - 1,

(6)

b1i = 1, 1 < i < m, bki = 0, 2 < k < m, 1 < i < k - 1,

i-1

bki = Pijbk-i,j, 2 < k < m, k < i < m, j=k-i

where Pij are numerical coefficients of schemes (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 Dahlquist's test equation

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

with complex X, Re(X)< 0. Variable X is considered as a certain eigenvalue of the Jacobi matrix of problems (1) or (3). Applying numerical scheme (4) to Dahlquist's equation we get

m m

yn+1 = Qm(z)yn, Z = hX, Qm(z) = 1 + ^ CmiZi, Cmi = ^ bijPmj, 1 < i < m. (8)

i=1 j=1

Denoting Cm = (cm1,..., cmm) and Pm = {pm1,...,pmm) , we can rewrite the latter equality

(8) in the form

BmPm Cm: (9)

where the elements of matrix Bm are defined in (6). For internal numerical schemes (4) we have ( ) ( ) k k

yn,k = Qk( z)yn, Qk( z) = 1 + ^3 CkiZi, Cki =^2 bij I3k+1,j, 1 < k < m - 1. (10)

i=1 j=1

Using designations pk = ( fik+1,1,..., fik+1ik)T and ck = (ck1, ...,ckk)T we obtain that coefficients ¡3ij of internal schemes (4) and the coefficients of corresponding stability polynomials are related by the equation

Bk Pk = ck, 1 < k < m - 1. (11)

It follows from of (6) and (10) that bki = ci-1,k-1, i.e., the elements of (k+1)-th column of matrix Bm are equal to coefficients of stability polynomial Qk(z). Hence, if the coefficients of the stability polynomials of the basic and intermediate numerical schemes are defined then the coefficients of methods (4) are uniquely determined from linear systems (9) and (11) with upper triangular matrices Bi, 1 ^ i ^ m.

Expansions of the exact and approximate solutions in the Taylor series in powers of h have the form

y(tn+1)= y(tn)+hf + y f 'f + O(h3),

m m (12)

yn+1 = yn+^2 bjPmo^j hf b2 jPmj^J h^f'jn + O{h?) ,

where the elementary differentials are taken with respect to exact y(tn) and approximate yn solutions, respectively. Comparing relations (12) under assumption that y(tn) = yn, one can

m

show that numerical scheme (4) has the first order of accuracy if bijpmj = 1. Hence, to

j=i

design m-stage methods of the first order of accuracy it is necessary to set cmi = 1.

2. Stability polynomials

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

k m

Qm,k (x) = 1 + ^ Oixi + ^ Ci%i, (13)

i=1 i=k+1

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

Denote extremum points of (13) as xi,..., x^fy—i, where x1 > x2 > • • • > xm—1. Unknown coefficients c^,k +1 ^ i ^ m can be obtained from the condition that polynomial (13) takes on predefined values at extremum points xi,k < i < m — 1, i.e.,

Qm,k{xi) = Fi, k < i < m — 1, (14)

where F(x) is a preassigned function, Fi = F{x^j. For this purpose let us consider the algebraic system of equations with respect to variables xi,k < i < m — 1, and cj,k +1 < j < m,

m

Qm,k(xi) = Fi, Q'm,k(xi) =0, k < i < m — 1, Qm,k = J2 i^xi—1 ■ (15)

1

We rewrite (15) in the form that is convenient for computations. Let us introduce vectors y, z, g and r with components

yi = xk+i-i, zi = ck+i, gi = Fk+i-i - 1 - Y cjyj, ri = - ^ jcjyj 1,

j=i j=i

1 ^ i ^ m - k,

We also introduce diagonal matrices E1,... ,E5 with elements on the main diagonal

(16)

4 = k + i, e2 = 1/yi, e3 = £ jcj yj-1 +$> + j )zj vk

%—k

--3- ¿^ jj yj—1 + Y1(k+j )zj yk+j—1 j=1 j=1

k m—k

e4 = £ j (j — 1)cj yj—2 + ^ (k + j)(k + j — 1)zj yk+j—2, j=2 j=1

(17)

e5

( — 1)k+i—1, 1 < i < m — k.

Consider matrix A with elements aij = y,k+j, 1 ^ i, j ^ m — k. The elements of vectors (16), matrices (17) and A depend on numbers m and k, where

g = g(y), r = r(y), E2 = E2(y), E3 = E3(y,z), EA = EA(y,z), A = A(y).

Then, we can rewrite problem (15) in the form

Az - g = 0, E2AElz - r = 0. (18)

System (18) is ill-conditioned that leads to some difficulties in applying the fixed-point iterations for its solution. For convergence of the Newton method it is necessary to obtain good initial values but in this case it is difficult problem in its own right.

If we assume in (15) that Fi = ( — 1)i, k ^ i ^ m — 1, we can find the polynomial with the maximal length of stability interval. In this case the problem of finding initial value y0 is solved by using values of the Chebyshev polynomial at extremum points over interval [—2m2,0], where m is the degree of polynomial (13). These values are

yi = m2[cos(i^/m) — 1], 1 ^ i ^ m — 1. (19)

Substituting (19) in system (17), we obtain coefficients of the Chebyshev polynomial for which |Qmi(x) | < 1 on x € [—2m2,0]. Then for any k values given in (19) can be taken as the initial values, and according to numerical calculations there is good convergence rate in this case. If Fi = ( —1)i, k < i < m — 1, then the choice of initial values is quite a difficult problem.

Let us describe a way to solve (18) that does not require good initial values. One can apply relaxation to solve system (18). The main idea of relaxations is to solve unsteady-state problem which solution converges to the steady-state solution of the initial problem. Let us consider the Cauchy problem

y' = Es( E2AEiA-ig — r), y(0)= yo. (20)

Apparently, upon finding the stationary point of (20), the coefficients of a stability polynomial can be determined from system (18). Let us notice that because matrix Es is used all eigenvalues of the Jacobi matrix of (20) have negative real parts, i.e., problem (7) is stable. It follows from numerical results that (20) is a stiff problem. Methods for solving such problems use calculation of the Jacobi matrix which cause difficulties in solving (20). Therefore, we apply the method of the second order of accuracy using numerical calculation and freezing the Jacobi matrix to solve (20) [5-6].

It can be shown that values of the polynomial coefficients tend to zero as m increases. Coefficients ci, k + 1 < i < m, were calculated for the polynomial degree up to m = 13 using algorithm [2]. Moreover, the algorithm of obtaining polynomial coefficients on the interval [—1,1] was described in [7]. In this case coefficients ci grow with slower rate and it is possible to construct polynomials of degrees m > 13.

3. Calculation of coefficients of stability polynomials using the GMP library

It is not difficult to see that coefficient cm of stability polynomial (13) tends to zero as m increases and in particular if m =13 and k = 1 the value of cm is about 10-26. Solution of problem (20) where m > 13 with double precision is very hard to realize because of roundoff errors. In order to compute coefficients of polynomial of higher degrees m algorithm was implemented using qd library [8, 9].

The qd library allows one to perform computations with higher accuracy. Standard data type double which allows one to perform computations with double precision is restricted by 53 bites of binary mantissa and provides accuracy of 16 decimal digits. Whereas qd data type dd_real has 106-bit mantissa that provides accuracy of 32 decimal digits. In fact, the number of data type dd_real is a programmed concatenation of two double numbers, where mantissa becomes doubled but the range of values that can be represented using new data type stays the same (from 10 -308 to 10308). Despite this restriction accuracy of number representation increases.

With the use of this library the coefficients of polynomials up to degree m = 35 were computed [8]. Nevertheless, the qd library has some disadvantages. Firstly, accuracy of number

representation is restricted because of program implementation of data types. Secondly, it can be used only in Unix systems. Moreover, the qd library is written in the C + + programming language. That is why numerical codes that use this library could be slower than codes written in low-level programming languages (for example, C).

Here we show numerical results of calculations of polynomial coefficients with the help of the library for arbitrary precision arithmetic GMP. This library provide accuracy of computations that is restricted only by the size of random access memory. It is cross-platform library, and it supports operations on integer, rational and real numbers. Besides, the GMP library is written in the C programming language which potentially increase the speed of computations.

Using the GMP library we obtain coefficients of polynomials up to degree m = 40. At higher degrees there are some difficulties that may be related to the choice of initial conditions for problem (20).

4. Construction of stability regions

Let us now describe the effect of function F on the size and shape of the stability region. If we assume Fi = ( — 1)®, k < i < m — 1, than the stability interval length is \jm| = 2m2. In this case, we have the maximum length of stability interval along the real axis for given m. The stability region of such methods is almost multiply connected which leads to the reduction of stability interval length because rounding errors may cause small imaginary parts of Jacobi matrix eigenvalues to appear. Fig. 1 shows the stability region of 5-stage method, where the stability interval length is |yto| = 50.

In order to avoid the stability region reduction because of rounding errors it should be "stretched" along the imaginary axis at the extremum points of the stability polynomial. For that we can assume F® = ( — 1)®^, 1 < i < m — 1, 0 < n < 1. For example, if we choose H = 0.95 then the stability interval length is reduced by only 3-4% in comparison with the maximal possible length that is equal to 2m2. Then it becomes equal to |yto| = 48.39 (Fig. 2). The stability region of the 5-stage method at n = 0.8 is shown in Fig. 3. In this situation, the stability interval length is reduced to |yto| = 43.55 with conjoined "stretching" along; the imaginary axis. For better visualization of the roots of polynomial (13) level lines |Qm,k (x)| = 1, ^mk^x^ = 0.8, Qmk(x) = 0.6, ^mk(x)! = 0.4, |Qm,k(x)| = 0.2. in the complex plane Are shown in all figures.

lm<z)

5

Fig. 2. Stability region at m = 5, k = 1, F = {—0.95,0.95,0.95, 0.95}, |7m| = 48.39

5. First-order method

For numerical solution of Cauchy problem (1) we consider the explicit five-stage Runge-Kutta method

yn+i = yn + piki + P2k2 + P3k3 + P4k4 + psks, ki = hf(yn), k2 = kf(yn + $2iki),

k3 = hf (yn + Aiiki + p32k2), (21)

k4 = hf(yn + ^4iki + ^42 k2 + ,

ks = hf (yn + Psiki + Ps2k2 + @S3k3 + ^S4k^,

where y and f are real N-dimensional vector functions, t is independent variable, h is the integration step, ki,k2,k3,k4 and ks are stages of the method, pi,p2,p3,p4,ps, P2i, P3i, P32, P4i, P42, P43, Psi, Ps2, Ps3, Ps4 are numerical coefficients that define accuracy and stability of (21). Applying the algorithm, we obtain coefficients of the stability polynomial:

csi = 0.1e1, cs 2 = 0.164341322127140896342e0, cs,3 = 0.948975952580473808808e — 2, cs,4 = 0.223956930863224544258e — 3, cs,s = 0.18509727522235334153e — 5

In this case the stability interval length is \Ym\ = 48.39. Upon solving (9) and (11), we obtain the coefficients of the first-order method:

321 = 0.0413243016210550, 331 = 0.0805823881610573, 332 = 0.0805823881610573, 341 = 0.11916681511228434, 342 = 0.1597820013984078, 343 = 0.0819394878966193, 351 = 0.1570787892802991, 352 = 0.2379583021959820, 353 = 0.1631711307360486, 354 = 0.0822916178203657,

p1 = 0.1945277188657676, p2 = 0.3151822878089125, p3 = 0.2437005934695969,

p4 = 0.1641555613805598, p5 = 0.0824338384751631.

Accuracy control of numerical scheme is based on local error estimation [10]. The magnitude of

An = [(0.5 - cm2)/«2p2 - ki) (22)

is used as preliminary estimation of local error. To we estimate the final accuracy the magnitude of

A'n = (0.5 - cm2)(hf (vn+1) —ki) (23)

is used. Thus, the following inequalities

An < e, An < e. (24)

are used for the accuracy control and for the choice of of integration step. As k1 linearly depends on integration step then omission of inequality (24) leads to just one additional computation of the right hand side of the problem. If the step of integration is successful the second inequality

(24) does not lead to the increase of computational cost because f(yn+^ is not used at the next step. At the same time if the second inequality (24) is used for accuracy control the repeat computations in the case of violating accuracy criterion are quite expensive. Moreover, the greater m the higher computational cost is. Nevertheless, in most cases preliminary estimation of An allows one to avoid repeat computations. The following inequality

Vn < Ym,i (25)

is used for stability control of method (2), where

«2332

1

max

[«2k3 + «3k2 — («2 + a3)k1]j/[k2 — k1]j , (26)

and positive constants Ym,1 define the size of stability regions [10].

6. Merson method

One of the most effective explicit fourth-order Runge-Kutta type methods is the Merson method [8]

121

Vn+1 = Vn + t k1 + - k4 + - k5,

636

k1 = hf (vn), k2 = hf (vn + 1 k1), k3 = hf (Vn + 1 k1 + 1 k2), (27)

( 13 ) ( 13 )

k4 = hf(vn + ^ k1 + ^ ky, k5 = hf(vn + ^ k1 — ^ k3 + 2k^.

Vn =

The fifth computation of function f does not result in the fifth order of accuracy but allows one to extend the stability interval to 3.5 and estimate truncation error 6n,4 using stages ki i.e.,

5nA = (2ki - 9k3 + 8k4 - 2k5)/30.

We use inequality ||£n,4|| ^ 5e5/4 for accuracy control. Despite the fact that inequality for accuracy control is obtained with the help of a linear equation, it shows high reliability in solving non-linear problems.

Now let us construct the inequality for stability control. Applying to k3 — k2 the first order Taylor's formula with the remainder term written in the Lagrangian form, we have

k3 — k2 = h[df (Mn)/dy](k2 — ki)/6,

where vector ¡in is calculated in some vicinity of solution y{tn). Taking into account that

k2 — ki = hfn fn/3 + O(h3),

the inequality

vn 4 = 6 • max

p _ p

k3 k2

kj — kj

< 3.5 (28)

can be used for stability control of (27), where 3.5 is the approximate length of stability interval. Let en,4 = Sn,4/5. Then inequalities en,4 < 5e5/4 and vn,4 < 3.5. can be used for accuracy and stability control of scheme (27), respectively.

7. Numerical results

The computations were performed on Intel(R) Core(TM) i7-8550U CPU. However, coefficients of stability polynomial were computed with the help of the GMP library whereas solution of differential problem was determined with double precision. The norm ||£n|| in the inequalities for the accuracy control was calculated by the formula

||£n|| = max |e\/(\ynI + r),

where i is a number of vector component, r is a positive parameter. If inequality \yln \ < r is satisfied for the component with number i then the absolute error r ■ e is controlled. Otherwise we control the relative error e, where e is the required accuracy.

We chose the Van der Pol oscillator (29) as a test example. This problem has the stiffness ratio approximately equal to 106:

y'l = y2, y2 =((1 — y2)y2 — yi)/10-6,

V J (29)

0 < t < 1, h0 = 10^3, yi(0)=2, y2(0)=0, e = 10~2.

The efficiency of two algorithms are compared. The first algorithm is the first order 5-stage Runge-Kutta method described in Section 5. The second algorithm is the traditional 5-stage Merson method of the forth order of accuracy (27). Both algorithms were applied in two modes: with stability control and without it. We counted total numbers of steps, repeat computations of a solution (due to omission of the defined accuracy), and the number computations of the

right hand side the problem. The accuracy e = 10-2 was supported with the Merson method, whereas for the first-order method we needed to use e = 10-5 in order to provide 10-2 in fact. Nevertheless, under these conditions the constructed algorithm shows better efficiency (Fig. 4).

The comparison of two algorithms shows that stability control increases efficiency because extra repeat computations of solution originating from instability of numerical scheme are eliminated. In addition, the constructed first-order algorithm has less computational cost estimated by the number of computations of the right hand side. Simulations of other test examples show similar pattern.

Criterion First-order method without stability control First-order method with stability control Merson method without stability control Merson method with stability control

Number of integration steps 69 433 51 414 549 241 556 114

Number of repeated solution calculations 20 001 1 052 187 120 6 464

Number of right part computations 452 683 309 948 3 494 685 2 806 426

Fig. 4. Numerical results for the Van der Pol oscillator problem

Conclusion

Implementation of the algorithm to obtain coefficients of stability polynomial with the use of the GMP library allowed one to build stability polynomial up to degree m = 40. It provides a possibility to develop methods with extended stability regions with respective number of stages. The greater number of stages the larger stability interval is, and therefore the higher efficiency of numerical scheme is achieved in the case of stiff problems.

Comparing two five-stage methods (the proposed first-order method and the Merson method), one can see that at the same number of stages extending the stability interval decreases the overall computational cost.

It is important to say that the first-order methods with extended stability regions allow one to significantly increase the efficiency in the region where the step is restricted by stability. So methods described here can be used in adaptive algorithms where number of stages may vary from one integration step to another. It provides large stability interval where it is needed and decreases computational cost when numerical scheme is unconditionally stable.

The study was funded by Russian Foundation for Basic Research (project no. 18-31-00375).

References

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

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

[3] E.A.Novikov, A.E.Novikov, Explicit-Implicit Variable Structure Algorithm for Solving Stiff Systems, International Journal of Mathematical Models and Methods in Applied Sciences, 9(2015), no. 1, 62-70.

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

[5] A.E.Novikov, E.A.Novikov, L-stable (2,1)-method for stiff nonautonomius problem solving, Computing technologies, 13(2008), 477-482 (in Russian).

[6] E.A.Novikov, Yu.A.Shitov, Integration algorithm for stiff systems based on a second-order accuracy (m, k)-method with numerical calculation of the Jacobi matrix, Krasnoyarsk: Preprint of the Exhibition Center of the Siberian Branch of the USSR Academy of Sciences no. 20, 1988 (in Russian).

[7] E.A.Novikov, M.V.Rybkov, The numerical algorithm of constructing stability polynomials of first order methods, Bulletin of the Buryat State University, (2014), no. 9-2, 80-85 (in Russian).

[8] E.A.Novikov, M.V.Rybkov, The numerical algorithm of constructing of stability regions for explicit methods, Control systems and information technologies, 55(2014), no. 1.1, 173-177 (in Russian).

[9] Yozo Hida, Xiaoye S Li, David H Bailey, Quad-double arithmetic: algorithms, implementation, and application, Technical Report LBNL-46996, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, 2000.

[10] L.V.Knaub, P.S.Litvinov, A.E.Novikov, M.V.Rybkov, Solving Problems of Moderate Stiffness Using Methods of the First Order with Conformed Stability Domains, University Scientific Journal, (2016), no. 22, 49-58.

[11] R.H.Merson, An operational methods for integration processes, Proc. of Symp. on Data Processing, Weapons Research Establishment, Salisbury, Australia, 1957.

Методы первого порядка с расширенными областями устойчивости для расчета задач электрических цепей

Михаил В. Рыбков Людмила В. Кнауб Данил В. Хоров

Сибирский федеральный университет Красноярск, Российская Федерация

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

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