ISSN 20700482, Mathematical Models and Computer Simulations, 2012, Vol. 4, No. 5, pp. 455–463. © Pleiades Publishing, Ltd., 2012. Original Russian Text © N.G. Bandurin, N.A. Gureeva, 2012, published in Matematicheskoe Modelirovanie, 2012, Vol. 24, No. 2, pp. 3–16.

A Method and a Software Package for Numerical Solution of the Systems of Nonlinear Ordinary IntegroDifferentialAlgebraic Equations N. G. Bandurin and N. A. Gureeva Volgograd State University of Architecture and Civil Engineering, Volgograd, Russia email: [email protected]; natalya[email protected] Received December 14, 2010

Abstract—A method, an algorithm and a software package for automatically solving the ordinary non linear integrodifferentialalgebraic equations (IDAEs) of a sufficiently general form are described. The author understands an automatic solution as obtaining a result without carrying out the stages of selecting a method, programming, and program checking. Both initial and boundary value problems for such equations are solved. It is assumed that the complete set of boundary and initial conditions at the beginning of the integration interval are given. By performing differentiation, the system of IDAEs can be modified, in general, into a system of ordinary nonlinear differential equations (IDEs). The problem of finding the solution of the abovementioned system on the uniform grid on the integration interval is posed in two forms: as solving the system of IDAEs and as solving the appropriate system of IDEs, where the developed program is to be used. In order to reduce the system of IDAEs and the sys tem of IDEs to the systems of ordinary nonlinear algebraic equations, at every stage of the algorithm the integration and differentiation formulas obtained earlier by N.G. Bandurin are used. Systems sim ilar to those test systems of both nonlinear IDAEs and IDEs considered in this investigation are solved by using the computer programs. It is evident that the coincidence of the results for one and the same system of equations in its different forms can serve as good evidence of the correctness of the obtained results. Keywords: integraldifferentialalgebraic equations, numerical methods, nonlinear problems, com puter programs DOI: 10.1134/S2070048212050031

1. INTRODUCTION The systems of nonlinear ordinary integrodifferentialalgebraic equations (IDAEs), which constitute the most general form of representation of equations depending on one argument, have good prospects of being used in the area of mathematical simulation of various physical processes in science and technology. To date, methods for solving the equations have been developed and the properties of the solutions of the systems of differential equations (DEs) of the form

y (n) = f[t, y(t ), y '(t ),..., y (n −1)(t )] = 0, each of which is solved with respect to the highest derivative [1–9], and also of the systems of integral equations (IEs) [10, 11], have been studied. Software packages for solving such systems of equations (Maple, Macsyma, Mathematica, Reduce, Matlab, Macsyma) are widely used. The solutions of many types of differential equations with singularities that do not allow obtaining acceptable results for practical use by standard means have been investigated [12–21]. In comparison with the system of differential equations, a more general form of the equations depending on one argument is the integrodifferential equations (IDEs); methods for solving them have been proposed in many works in Russia and abroad, for example, in [22–26]. A large number of publications are devoted to solving rigid differential equations and the problems with boundary layer [27–30]. The highest generality is obviously possessed by the systems of IDAEs; a method of their numerical solution is given, for example, in [31]. 455

456

BANDURIN, GUREEVA

2. STATEMENT OF THE PROBLEM The problem of obtaining the solution on the uniform grid t1 = 0, t 2,..., t M = b of the segment [0, b] of the following system of N ordinary nonlinear, in the general case, IDAEs is formulated:

F[t, y(t ), y '(t ),..., y (n)(t ), I(t )] = 0 . Here, F is a column vector of size N;

(1)

(y(t ))T = (y1(t ), y 2(t ),..., y N (t )); (I(t ))T = (I 1(t ), I 2(t ),...); x

∫

I i (t ) = ψ i[t, u, y(t ), y '(t ),..., y(u), y '(u)...]du. 0

The results are obtained in the form of tables of values of functions and their derivatives at the grid nodes. Particular kinds of system (1) can be the systems containing, in various combinations, integrals, deriv atives, and algebraic expressions in terms of the unknown functions. The experience of using these pro grams has shown that in the case of correctly written equations, which adequately simulate real physical pro cesses, it is possible to get a sufficiently precise solution. 3. METHOD OF SOLUTION In this work, the method of steps, according to which the entire interval of integration [0, b] should be presented in the form of Q large steps of length τ, that is, we take b = Qτ , is used. For numerical imple mentation of the method, each large step is subdivided into m nodes, and the entire computation process consists of the consecutive performance of Q stages. As a result of carrying out each stage, the values of the sought functions at m nodes of the corresponding large step are obtained. While carrying out the next stage, the values of the functions calculated at the end of the previous one are taken as the initial condi tions. To reduce the system of IDAEs to a system of usual nonlinear algebraic equations, at each stage of solution, the integration and differentiation formulas which were first published in [32], are utilized in the algorithm. Below, we present, in a concise fashion, the derivation of these formulas as applied to the for mulated problem. To obtain an integration formula, we assume that on segment [c, d ] of the real axis t a uniform grid is given; at its P nodes the values of the k th order derivative of the function y = f (t ), defined on [c, d ], and the values of the lower order derivatives at the initial node t1 = c, are known. The problem of the calcula tion at the same nodes of the values of the function, as well as its derivatives and integrals, expressed by means of the nodal values of the derivative of the kth order and the derivatives at t1 = c, is formulated. We interpolate the derivative y (k )(t ) on the segment [t1, t n ] (n ≤ P ) by a polynomial of degree (n – 1) [32−34]:

y (k )(t ) = t T c,

(2)

where t T = (1, t,..., t n−1); c T = (c1, c2,..., cn). Setting t = t1, t = t 2,..., t = t n, in (2), we can obtain the matrix relation

y (k ) = Z c,

(3)

where y (k ) is a vector, whose components are the nodal values of the kth order derivative of function f(t), whereas the elements of matrix Z are equal to z ij = t i j −1. The determinant of matrix Z is the Vandermonde determinant; therefore, we can write

c = Z −1y (k ),

(4)

y (k )(t ) = t T Z −1y (k ).

(5)

and Eq. (2) assumes the form Integrating the last equality successively over the segments [t1, t 2 ] ,..., [t1, t n ] , we can write

y (k −1) = I y1(k −1) + BZ −1y (k ). MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

(6) Vol. 4

No. 5

2012

A METHOD AND A SOFTWARE PACKAGE

457

Here, y (k −1) is a vector whose components are the nodal values of the derivative of the (k – 1)th order; I is a vector of size n, whose elements are equal to one; and the elements of the matrix B are equal to bij = (t i j − t1j )/ j. From (6), we obtain a formula of integration over the segment [t1, t n ]

y (k −1) = I y1(k −1) + s y (k ),

(7)

where s = BZ −1 is a square matrix of integration of order n. If P > n, then, applying successively formula (7) for the segments [t j(n −1)+1, t j(n −1)+ n] ( j = 0,1,2,...), we can obtain the components of vector Y (k −1) of the derivative at all the nodes of the segment [c, d ]:

Y (k −1) = I y1(k −1) + S Y (k ), where S is a square matrix of integration of order P. For example, for c – d = 3, n = 3, P = 4, this matrix has the form

(8)

0 0 0 ⎤ ⎡ 0 ⎢0.37500 0.79167 − 0.20833 0.41667⎥ ⎥. S=⎢ 0 ⎥ ⎢0.33333 1.33333 0.33333 ⎢⎣0.37500 1.12500 1.12500 0.37500⎥⎦ Having performed the integration r times according to (8), we obtain the formulas to calculate the deriv atives at the interpolation nodes, expressed through the values of the derivatives of order k at the same nodes and the derivatives of lower orders at the initial node x = x1: Y

(k − r )

r −1

=

∑S

m

I y1(k

− r + m)

+ S r Y (k )

for

r ≤ k.

(9)

m =0

If r > k, the components of the vector Y (k − r ) in (9) will be the values of integrals, whereas the formula will assume the form

Y

(k − r )

r −1

=

∑

Sr

−k + m

I y1(m) + S r Y (k ).

m=r −k

It is taken into account in the last equality that y1(k − r ) ≡ 0 for r > k, because, according to the accepted notations, these quantities are the values of the integrals at the node t = t1. To obtain a differentiation formula, we represent the integration matrix in the form

⎡s s ⎤ (10) S = ⎢ 11 12 ⎥ , ⎣s 21 s 22 ⎦ where s11 = 0, the elements of row matrix s12 are equal to zero, whereas the square matrix s 22 of order (P − 1) is nonsingular. From this, it is possible to obtain the expression of the vector of the nodal values of the derivative of function y(t ) in terms of the initial values of the derivative and the nodal values of the function itself: Y ' = ( E − DS ) I y11' + D Y,

(11)

0⎤ ⎡ 0 where D = ⎢ −1 is the matrix of differentiation and E is the unit matrix. −1 ⎥ ⎣−s 22 I 2 s 22 ⎦ For c – d = 3 < n = 3, P = 4 this matrix has the form

0 0 0 ⎤ ⎡ 0 ⎢− 0.9444 0.5000 − 0.5000 −0.0555⎥ ⎥. D=⎢ ⎢ 0.7778 −2.0000 1.0000 0.2222 ⎥ ⎢⎣− 2.1667 4.5000 − 4.5000 2.1667 ⎥⎦ MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012

458

BANDURIN, GUREEVA

Having performed the differentiation r times according to (11), we obtain the formula for the calculation of the derivatives of the r th order at the interpolation nodes, expressed in terms of the derivatives at the initial point t = t1 and the values of the function at these nodes: r −1

Y

(r )

=

∑D

m

( − m)

(E − DS )I y1r

+ D r Y.

(12)

m =0

For any integer p ≥ 0, the following equalities, which express the results of integration and differenti ation of the function f(x) = 1 hold: P

∑ j =1

P

p Sij

= (t i − t1) p !, p

∑D

p +1 ij

= 0(i = 1,2,..., P ).

j =1

In contrast to other known numerical methods and algorithms, the use of formulas (9) and (12) in numerically solving the boundary value and the initialboundary value problems for the systems of IDAEs gives the possibility (1) to approximate with high precision IDAE systems, having sufficiently smooth solutions, by their discrete analogues, since the formulas are exact on the segment [c, d ] for all the polyno mials of degrees up to n; (2) to solve the equations with the given values of derivatives not only on the boundaries but also inside the integration domain, since in the derivation of the formulas, we can take any node t i (i = 1,2,..., P ) as the initial point on the segment [c, d]; and (3) to exclude the necessity of carrying out nonstandard algebraic transformations on the boundary of the domain, while solving the IDAEs of a high degree using interpolation by polynomials, as recommended, for example, in [35]; this is because the obtained formulas, (9) and (12), do not depend on the indices of the nodes that are situated outside of the integration domain, and explicitly contain the derivatives necessary for setting the boundary conditions. Utilizing Taylor’s formula, we can obtain the estimate O(h n +1) for the approximation error of the function and its derivatives, where h is the grid step. It should be indicated that the numerical algorithm can be based both on the integration formula (9) and the differentiation formula (12). The comparative calculations show that, in solving nonrigid correct problems, formula (9) ensures greater precision than formula (12) does. At the same time, while solving rigid problems, the calculation process based on formula (12) is more stable in comparison to (9). 4. THE SOFTWARE PACKAGE The software package considered below consists of two separate programs: the CASHY1 program, designed for numerical solving the initial value problems for systems of nonlinear IDAEs (possibly, with a delayed argument), and the BOUNDS1 program for solving the boundary value problems for IDAEs. The delay value τ in the CASHY1 program is assumed to be the same for all functions. The programs work on the basis of the minimum of the initial data, since, to start the iteration process, it suffices to input into the computer memory the texts of Eqs. (1), expressions under the integral sign, the initial (boundary) condi tions, and the necessary numerical parameters (selection of the solution method and the programming and debugging of each specific program are eliminated). In the CASHY1 problem, the initial continuous functions for all the sought functions on the delay interval of length τ should be defined. Besides the usual derivatives and integrals, the equations may also contain fractional derivatives and integrals [36, 37]. On various sections of the integration axis the system of equations may have different structures under the condition that the order of the highest derivatives of the functions remains the same for all sections. In the programs an interpolation procedure on the basis of a complete polynomial up to the ninth degree (n = 9) is implemented; thus, the precision is O(h10 ). No further increase of n is allowed in the program, since, in this case, the process of inversion of matrix Z in (3) becomes unstable, which is connected with the limited precision of the arithmetic unit of a PC. Because of the same reason, the excessive refinement of the grid becomes impossible. One of the specificities of the method and the programs is that there is a possibility to calculate, to print, and save the graphs of the dependences y = f (t ), y' = f (t ), y i' = f (y i ), y i = f (y j ), while the known numerical methods and programs for solving the nonlinear Cauchy problem are based on the calculation of the nodal values of the sought functions [5–8]. The main restrictions are the following: the system must be correct; the maximal number of the equa tions being solved is 9; and the maximal order of the derivative functions is 9. MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012

A METHOD AND A SOFTWARE PACKAGE

459

5. THE TEST EXAMPLES 5.1. The boundary value problem for the system of nonlinear IDAEs depends on two functions u(x) and v(x) and two unknown numerical parameters p1 and p2: 2

t

∫

2

∫

t

∫

∫

u "+ [t + u "(t ) v(z)dz]dt + p1x [t + 0,5 zv(z)dz]dt − x 3 6 − x 2 2 − 13 4 x + p1 − 19 6 = 0, −1

−1

−1

2

t

∫

∫

2

v − [t + 0,5u(t ) v(z)dz]dt = 0, −1

−1

p1 − p23 + 6 = 0,

∫

p2 + [tu'(t ) + p1 − 2]dt − 2 = 0.

−1

−1 2

∫

The boundary conditions are u(–1) + [tu'(t ) + p1 − 2]dt − p1 − 5 = 0, u(2) − 2 p1 = 0. −1

The solution is to be found on the segment [–1, 2]. The exact solution is u(x) = x 2, v(x) = x, p1 = 2, and p2 = 3. The initial approximation is u0(x) = 0, v 0(x) = 0, p10 = 1, and p20 = 1. Since the exact expressions for the sought functions are polynomials of not more than second degree, the usage of interpolation polynomials of a higher degree provides a zero error of the applied method; and, as a result of solving the problem with h ∈ [0.03;0.75], resulted in a deviation of the approximate solution from the exact one of not more than 10–12. 5.2. Cauchy problem for the system of IDEAs with a delayed argument; the delay value is τ = 1. t

∫

y1'''(t ) + [y3(u) + y 2(u − τ)]du +y3(t ) − sh t − t 2 + 2t − 1 = 0, 2

0

y 2(t ) + y1"(t ) + y 4(t ) − ch t − exp t − t = 0, t

∫

y3(t ) − y 4(u)du + y 4(t ) + sin[y1(t − τ)] − exp(−t ) − sin[ch(t − τ)] − 1 = 0, 0

t

∫

y 4(t ) − y 4(u)du − y3(t ) + exp(−t ) − 1 = 0. 0

The initial conditions are y10 = 1, y 20 = 0, y30 = 1, and y 40 = 1. The exact solution is y1(t ) = ch t, y 2(t ) = t, y3(t ) = exp(−t ), and y 4(t ) = exp(t ). The solution is to be found on the segment 0 ≤ t ≤ 10. As a result of solving, we obtained the following maximal deviation values Δ:

Δ h=0.2 = 9.1 × 10 −6, Δ h=0.1 = 1.8 × 10 −8, Δ h=0.05 = 6.0 × 10 −9, Δ h=0.025 = 7.6 × 10 −9. 5.3. Boundary value problem for the system of three IDEAs with unknown functions u(x), v(x), and w(x): x

∫

u "+ [v (s)w(s) + 0

π 2

π 2

∫ w(t)dt]ds − x

3

6 = 0, v +

0

x

∫ u(s)w(s)ds − 1.5 = 0,

∫

w − u' + u(s)ds + cos x − 1 = 0.

0

0

The boundaries of the integration segment are x 0 = 0, x N = π 2. The boundary conditions are π 2

u(0) +

∫ 0

π 2

u(s)ds − 1 = 0,

u(π 2) +

∫ u'(s)ds − 2 = 0,

v(0) = 0.

0

The exact solution is u(x) = sin x, v(x) = x, w(x) = cos x. The initial approximation is u0(x) = 0, v 0(x) = 0, w0(x) = 0. MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012

460

BANDURIN, GUREEVA

As a result of solving, for various numbers N, we obtained the following values of the maximal deviation Δ: of the functions of the approximate solutions from the exact one:

Δ N =5 = 3.6 × 10 −8, Δ N =11 = 2.4 × 10 −8, Δ N =51 = 8 × 10 −14, Δ N =151 = 7 × 10 −14. One can see that the precision does not increase under excessive refinement of the grid, which can be explained by the accumulation of the influence of the number of rounding errors with the increase of the number of arithmetic operations. 5.4. Boundary value problem for a system of nonlinear algebraic equations with unknown functions u(x), v ( x ), and w(x). The solution is to be obtained on the segment [0, 2] in two forms: as a result of solving the algebraic equations (AEs) or the system of differential equations (DEs) obtained by differentiating (13) with respect to x. On various sections of the integration segment, the structure of the equations is different, namely: for 0 ≤ x ≤ 0.5,

u + w − sin x − cos x + x − 1 = 0,

v + u 3 − (sin x − x)3 − x = 0,

w − cos v − 1 = 0;

(13)

for 0.5 ≤ x ≤ 2:

u + sin v − 2 sin x + x = 0, v + w − cos x − x − 1 = 0, w + sin w − sin(cos x + 1) − cos x − 1 = 0. After differentiation, we obtain the system of DEs: for 0 ≤ x ≤ 0.5, u' + w' − cos x + sin x + 1 = 0,

v' + 3u 2u' − 3(sin x − x)2(cos x − 1) − 1 = 0,

w' + sin vv' = 0; for 0.5 ≤ x ≤ 2:

u' + cos vv' − 2 cos x + 1 = 0,

(14)

v' + w' + sin x − 1 = 0, w' + cos ww' + sin x(cos(cos x) + 1) = 0. The structure of the system of equations at the junction point of two sections at x = 0.5 ensures the continuity of the functions and their derivatives. The exact solution is u(x) = sin x − x, v(x) = x, w(x) = cos x + 1. The initial approximation is u0(x) = 0, v 0(x) = 0, w0(x) = 0. As a result of solving system (13) as algebraic equations, we obtained the deviation values of the approx imate solution from the exact one, which practically did not depend on the grid step and were roughly equal to 10–12, which can be explained by the fact that at the junction border of two sections, at x = 0.5, the continuity of the functions and their derivatives is observed, whereas the error is only caused by the inaccuracy of the arithmetic operations. In solving the system of differential equations (14) for various val ues h of the grid step, the following deviations were obtained:

Δ h=0.25 = 6.6 × 10 −9, Δ h=0.2 = 1.7 × 10 −10, Δ h=0.125 = 1.1 × 10 −12, Δ h=0.1 = 5.5 × 10 −12. In solving this differential problem, the error depends on the grid step, because now the inaccuracy of the interpolation procedure exerts its influence, whereas some increasing of the deviation Δ under the grid nodes’ concentration is connected with the errors in performing a large number of arithmetic operations. 5.5. Boundary value problem for a system of IDEAs with unknown function u(x), v(x), w(x): π /2

x

∫

u" + [v '(t )w(t ) + 0

t

∫ w(z)dz ∫ v(z)dz]dt − x 0

3

6 = 0,

0

π 2

v' +

∫ u(t)w(t)dt − 2.5 = 0,

(15)

0

MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012

A METHOD AND A SOFTWARE PACKAGE

461

x

∫

w − u' + u(t )dt + cos x − 1 = 0. 0

π 2

The boundary conditions are u(0) +

π 2

∫ u(t)dt = 0, u(π 2) + ∫ u'(t)dt = 0, v(0) = 0. 0

0

The exact solution is u(x) = sin x, v(x) = x, w(x) = cos x. The initial approximation is u0(x) = 0, v 0(x) = 0, w0(x) = 0. The values of the maximal deviation Δ are Δ h=π 8 = 6.0 × 10 −4,

Δ h=π 20 = 2.4 × 10 −8,

Δ h=π 40 = 1.4 × 10 −11,

Δ h=π 60 = 2.7 × 10 −13.

An attempt at solving the system of equations, obtained as a result of differentiation of (15), ended in failure, since the iteration process for the new problem turned out to be highly unstable. This effect is caused by the specificities of the structure of equations and corresponds to a general tendency, according to which the difficulties in obtaining numerical solution of differential equations increase with the increasing of the order of derivatives. 5.6. Cauchy problem for a system of nonlinear IDAEs with a delayed argument; the delay value is τ = 1. y1'(t ) + sin[ y 3(t )] + exp[ y 3(t − τ)] − exp(t ) − exp[sin(t − τ)] − sin(sin t ) = 0, t

∫

y 2''(t ) + [y 2' (u)w3(u)]du − exp(−t )(1 + (sin t + cos t )) 2 − 0.5 = 0, 0

y3(t ) + y1(t ) + exp[y3(t )] − sin t − exp(t ) − exp(sin t ) = 0, t

∫

y 4(t ) + [y1(u) + y 2(u) + y 2(u − τ)]du − cos t − exp t + exp 2(exp(−t ) − 1) = 0. 0

The initial conditions are y10 = 1, y20 = 1, y30 = 0, y 40 = 1. The exact solution is y1(t ) = exp t, y 2(t ) = exp(−t ), y3(t ) = sin t, y 4(t ) = cos t. The solution is to be found on the interval 0 ≤ t ≤ 10. It should be noted that, to ensure the reliable start of the iteration process, the input is provided in the program as the initial data of the initial values of the functions, whose derivatives are absent in the system of equations. As a result of solving, the maximal values of deviation Δ of the approximate solution from the exact one were obtained for the first function. These deviations for various values h of the grid step turned out to be the following:

Δ h=0.2 = 6.8 × 10 −6, Δ h=0.1 = 1.1 × 10 −8, Δ h=0.05 = 2.8 × 10 −9, Δ h=0.025 = 3.7 × 10 −9. The iteration process of this system is stable, and convergence is observed in a sufficiently wide interval of h. 5.7. Cauchy problem for a rigid differential equation [28, 29]

εy' + y + t − 1 = 0, y(0) = 0.

(16)

The exact solution is y(t ) = (1 + ε − t ) − (1 + ε)exp(− t ε). The solution is to be found in the segment 0 ≤ t ≤ 1. We obtained the values of the maximal deviation Δ for ε = 0.001:

Δ h=0.002 = 7.4 × 10 −3,

Δ h=0.001 = 7.7 × 10 −5,

Δ h=0.005 = 3.0 × 10 −7,

Δ h=0.0001 = 2.5 × 10 −11.

The maximal value of the function y max = 0.993091 at t = 0.0069 is automatically printed. As a result of solving, we obtained Δ h=0.0001 = 7.8 × 10 −11 for ε = 0.0005. The method and algorithm also allow obtaining the solution for smaller values of ε . However, it will require further refining of the grid, which is impossible for this program, as, due to the limitations of the MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012

462

BANDURIN, GUREEVA

Delphi algorithmic language and the specificities of the algorithm, we are not supposed to use more than 10000 nodes on the t axis. As a result of integration and differentiation of (16) it is possible to get the following equivalent prob lems: t

∫

ε y + (y(s) + s − 1)ds = 0; 0

ε y "+ y' + 1 = 0, y(0) = 0, y'(0) − 1 ε = 0;

ε y'" + y" = 0, y(0) = 0, y'(0) − 1 ε = 0, y"(0) + (1 + ε) ε 2 = 0. For the last three problems, we again obtained the same values of the deviation Δ as for (16). 5.8. Boundary value problem for the rigid secondorder differential equation [30] ε y" + xy' + y = 0, y(0) = 1, y(1) − exp(− 0.5ε) = 0, ε = 1 300. The following problems are equivalent problems:

(17)

ε y'" + xy" + 2 y' = 0, y(0) = 1, y(1) − exp(− 0.5ε) = 0, y'(0)=0;

(18)

ε y (IV ) + xy'" + 3y" = 0, y(0) − 1 = 0, y(1) − exp(−0.5ε) = 0, y'(0) = 0; y′′(0) + 1 ε = 0.

(19)

The exact solution is y(t ) = exp(−t 2ε). As a result of solving problem (17), we obtained the values of the maximal deviation Δ for various num bers of nodes m on the segment [0,1]: 2

Δ m=25 = 0.043,

Δ m =101 = 9.3 × 10 −8,

Δ m=151 = 1.0 × 10 −9.

For problem (18), Δ m=25 = 0.24, Δ m =101 = 6.2 × 10 −7, and Δ m=151 = 1.5 × 10 −9. For (19), Δ m =25 = 4.86, Δ m=101 = 2.3 × 10 −5, and Δ m=151 = 1.6 × 10 −8. It can be seen that under a sufficiently fine grid, we obtained satisfactory results for all three versions of this problem and the previous one, which testifies to the reliability of the method, algorithm, and pro grams in solving a certain class of rigid problems as well. REFERENCES 1. C. Runge, Math. Ann. 46, 167–178 (1895). 2. W. Kutta, Z. Math. Phys. 46, 435–453 (1901). 3. The Modern Numerical Methods of Solving Ordinary Differential Equations, Ed. by J. Hall and J. Watt (Mir, Mos cow, 1979) [in Russian]. 4. Ts. Na, Computational Methods of Solving the Applied Boundary Value Problems (Mir, Moscow, 1982) [in Russian]. 5. A. M. Samoilenko and N. I. Ronto, NumericalAnalytical Methods of Studying Boundary Value Problems (Nauk ova Dumka, Kiev, 1986) [in Russian]. 6. H. H. Robertson, The Solution of a Set of Reaction Rate Equations / Numerical Analysis: an Introduction (Aca demic Press, 1966). 7. T. Alishenas and O. Olafsson, BIT 34, 455–483. 8. J. C. Butcher, Math. Comput. 18, 50 (1964). 9. C. W. Gear, “The Automatic Integration of Stiff Ordinary Differential Equations”, in Information Processing 68: Proceedings of the IFIR Congress 1968 (North Holland Publ. Co., 1969), pp.187–193. 10. A. F. Verlan’ and A. F. Sizikov, Integral Equations: the Methods, Algorithms and Programs. A Reference Manual (Naukova dumka, Kiev, 1986) [in Russian]. 11. A. B. Vasil’eva and N. A. Tikhonov, Integral Equations (Fizmatlit, Moscow, 1984) [in Russian]. 12. D. F. Davidenko, Dokl. Akad. Nauk SSSR 88 (4), 601 (1953). 13. A. N. Tikhonov, Mat. Sb. 31, 575 (1952). 14. M. O. Korpusov and A. G. Sveshnikov, Zh. Vychislit. Matem. Matem. Fiz. 45 (1), 145 (2005). 15. N. S. Bakhvalov, Zh. Vychislit. Matem. Matem. Fiz. 9 (4), 842 (1969). 16. N. V. Shirobokov, Zh. Vychislit. Matem. Matem. Fiz. 49 (6), 1080 (2009). 17. E. B. Kuznetsov and V. I. Shalashilin, Zh. Vychislit. Matem. Matem. Fiz. 33 (12), 1792 (1993). 18. S. D. Krasnikov and E. B. Kuznetsov, Zh. Vychislit. Matem. Matem. Fiz. 45 (12), 2148 (2005). MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012

A METHOD AND A SOFTWARE PACKAGE 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37.

463

E. B. Kuznetsov and V. I. Shalashilin, Izvestiya Vuzov. Matematika, No. 11, 56 (1998). E. B. Kuznetsov, Zh. Vychislit. Matem. Matem. Fiz. 44 (9), 1540 (1993). N. N. Kalitkin, Numerical Methods (Nauka, Moscow, 1978) [in Russian]. L. M. Skvortsov, Mat. Model. 14 (2), 3 (2002). V. F. Antyushin and A. V. Budanov, et al., Vestnik VGU, Ser. Fiz. Matem, No. 1, 51 (2004). A. R. Validi and E. Balolian, Int. Journal of Math Analysis 3 (36), 1769 (2009). E. Celir, M. Bayram, and T. Yelogly, Int. J. Pure and Appl. Math. Sciences 3 (1), 93 (2006). B. Batiha, M. Noorani, and I. Hashim, Int. J. Open Problems Comput. Math. 1 (1) (2008). E. Hairer and G. Wanner, Solving Ordinary Differential Equations. Rigid and DifferentialAlgebraic Problems (SpringerVerlag, Berlin, 1996). V. G. Zverev and V. D. Gol’din, Vychisl. Tekhnologii 13 (3), 54 (2008). B. M. Bagaev and V. V. Shaidurov, Grid Methods of Solving the Problems with Boundary Layer (Nauka, Novosi birsk, 1998) [in Russian]. E. Doolan, J. Miller and W. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layer (Boole Press, Dublin, 1980). S. S. Dmitriev and E. B. Kuznetsov, Zh. Vychislit. Matem. Matem. Fiz. 48 (3) 430 (2008). N. G. Bandurin, Vychisl. Tekhnologii 7 (2), 3 (2002). N. G. Bandurin and V. A. Ignat’ev, Mat. Model. 19 (2), 105 (2007). N. G. Bandurin, Vychisl. Tekhnologii, No. 3, 31 (2010). G. A. Corn and T. M. Corn, Mathematical Handbook for Scientists and Engineers (McGrawHill, New York, 1961). A. M. Nakhushev, Fractional Calculus and Its Application (Fizmatlit, Moscow, 2003) [in Russian]. N. G. Bandurin and V. A. Ignat’ev, “Software Complex for Solving the Systems of Nonlinear IntegroDifferen tial Equations (IDE), Involving Fractional Derivatives and Integrals”, in: Proceedings of the IX International Conference “Intellectual Systems and Computer Sciences” 2006, MGU, vol. 2, part 1, pp. 56–58.

MATHEMATICAL MODELS AND COMPUTER SIMULATIONS

Vol. 4

No. 5

2012