Horner’s Scheme for investigation of solutions of differential equations with polynomial right-hand side Текст научной статьи по специальности «Математика»

Horner's Scheme for investigation of solutions of differential equations with polynomial right-hand side1

Alexander P. Afanas'ev

Head of Laboratory of Distributed Computational Systems

Institute for Information Transmission Problems RAS;

Professor, Head of Department of High-performance Computations

of Institute for Information Transmission Problems RAS

National Research University Higher School of Economics

Address: 19, build.1, Bolshoy Karetny Per., Moscow, 127051, Russian Federation

E-mail: apa@isa.ru

Sergei M. Dzyuba

Professor, Department of Information Systems Tver State Technical University

Address: 22, Afanasiy Nikitin Embankment, Tver, 170026, Russian Federation E-mail: sdzyuba@mail.ru

Irina I. Emelyanova

Senior Lecturer, Department of Information Systems Tver State Technical University

Address: 22, Afanasiy Nikitin Embankment, Tver, 170026, Russian Federation E-mail: emelyanova-123@yandex.ru

Abstract

We present a method for investigating solutions of systems of ordinary differential equations with polynomial right-hand side. Similar systems are of long-term interest for applications, because many process models have different physical, biological and economical natures described by these systems. The standard methods of numerical analysis are usually applied obtaining system solutions with the polynomial right-hand side, disregarding the specific form of the right-hand side. We suggest a different method starting from the fact that the right side of the equation appears to be a multidimensional polynomial. The relative simplicity of the right-hand side of the system under consideration made it possible to construct by this method approximate analytic solutions in the form of functions not only of time but of the initial conditions as well. In contrast to the majority of known methods, the latter made it possible in many cases to directly trace the systematic computational error. The implementation of the method is based on the construction of a discrete dynamical system along the solutions of the original system with subsequent use of the generalized Horner's Scheme. The computation peculiarity of Horner's Scheme lies in the fact that in many cases the scheme allows us to reduce the number of machine operations required for computation of the polynomial in comparison with the ordinary computation process. The second peculiarity of the generalized Horner's Scheme lies in the fact that there is a good decomposition of computation process that allows us to make calculations in parallel on independent nodes. According to computational experiments, this enables us to reduce computation time hugely even in the simplest cases while retaining required accuracy.

Key words: differential equations with polynomial right-hand side, Horner's Scheme, interactive system along solutions, MathCloud system.

Citation: Afanas'ev A.P., Dzyuba S.M., Emelyanova I.I. (2017) Horner's Scheme for investigation of solutions of differential equations with polynomial right-hand side. Business Informatics, no. 2 (40), pp. 33—39. DOI: 10.17323/1998-0663.2017.2.33.39.

1 The research was supported by the Russian Scientific Foundation (project No. 16-11-10352)

C

Introduction

onsider the system of ordinary differential equations written in the vector form

(1)

where x = (x1.....x") is a real vector function of the real

variable i and /=(/',•••,/") is a real vector function in which each element f' is a multidimensional polynomial in the variables x1, ..., x".

The system (1) has a special meaning in our days. This is explained by numerous models of processes of various physical, economic and other nature as described by such systems (e.g., [1—4]). The special significance of the system is that the strange attractors are quite common.

As is well known, each attractor is a compact invariant set (e.g., [1]). The classical general theory of dynamical systems says that every compact invariant set contains a compact minimal set (e.g., [5]). Therefore, for understanding the structure of strange attractors, constructing the minimal sets is extremely important.

The problem of constructing minimal sets has been posed in the general case only in [6]. As is clear from the results [6], the construction of minimum sets supposes the construction and investigation of the system (1) solutions on infinite horizon for which purpose we need to implement the construction of a certain discrete-time system in the line of primary solutions (see also [7]).

The aim of the present paper is to develop a method for constructing and investigating solutions of system (1) in the MathCloud system. The method is based on a special construction of local solutions of the system taking into account the fact that f is a multidimensional polynomial [7, 8]. In contrast to the majority of the known methods (e.g. [9—13]), the latter made it possible in many cases to directly trace the systematic computational error. Such tracing is extremely important, e.g., for studying the equations of chaotic dynamics, because computation accuracy in this case plays the decisive role and computations should be performed in rather long time intervals.

1. Construct local solutions of the system (1)

First of all, consider the problem of constructing a local solution x(t) of system (1) with the initial condition

*(0) = x0.

(2)

Acting in the standard way, replace system (1) by the integral equation

i

X(Í) = X0+J/(X(t))Í/T

(3)

To find a solution x(f) of Equation (3), we use the Picard successive approximation method and write

i

xn+l (t) = Xo+jf{XniT))dT-

(4)

Assume that

x1(t) = x0.

Then

I

*2 (0 = *0 + J/(*0 )dT = X0+f(Xo) t-0

It follows that equality

t e2

x3(t) = x0+jf(x0+ f(x0) t) dT=x0 + Yjp2i (x0 ) t

o ¡=i

holds, and 02 is a positive number depending on the polynomial f and the ip2 i are the corresponding real vector functions defined as multidimensional polynomial at the variable x0.

Acting formally by induction, using equation

= (5)

/=i

where dN is a positive integer depending only on N and the polynomial f and the ipNt are the corresponding real vector functions defined as multidimensional polynomial at the variable x0.

The question of the convergence of the sequence (xn)n<= n to the solution x(t) requires additional investigation. The answer can be found in Theorem 1 (the proof of Theorem 1 is contained in paper [7]).

Theorem 1. Let x0eR be a given point, and let a be some positive number. Let

m = max |/(x)|. Then, for all values T > 0 satisfying relation

M

Hold the following equality

lim

jv—>+<*>

= 0.

According to Theorem 1, formula (5) gives approxi-

mate analytical solutions of system (1) with initial condition (2) defined on [— T, 7] . This approximate analytical solution is a multidimensional polynomial of time t and vector jc0 which determine the initial state of the system. If the function f is nonlinear, then an instinctive feature of this representation consists in equality

V N ~

lim — = 0.

iv-m- eN

More than all, equality

Hm (0w+1-0w)=+~

is always executed.

2. Interactive system along solution (1)

Generally speaking, investigation attractors of system (1), if, of course, they exist, suggests the continuation of local solutions to the right on a sufficiently large time interval. In paper [7] one was invited to implement the construction of discrete dynamical system along solution jc(i).

Let x(t) be the solution of system (1) with the initial condition (2) defined for all t > 0 and is bounded for these t. Fixed real number n satisfies the relation

a = 2iMsup|x(i)-x0|.

Let

2m

(6)

(7)

where m is the number satisfying the conditions of Theorem 1.

By gt we denote the phase flow for which the field f if the phase velocity field. Then the relation t > 0

x(t + kT) = g'x(kT), £ = 0,1,...

holds for all t > 0. One can readily see that, by virtue of (6), (4), and (7), for te [0, T], the points g'x(kT) lie in the closed ball

2 1 and hence in the closed ball

Therefore, by Theorem 1 we see that approximation g'N to the operator g t satisfies equality

<>N

Now acting as the paper [7] it is not difficult to prove that it holds.

Theorem 2. Let conditions (6), (7) be satisfied for some n > 1, and number m satisfies conditions of Theorem 1. Then, for each positive number e, there exists a positive integer Nc such that, for N > Nc, the inequality

holds for all g0 e Ba (x0).

2

Therefore, for all sufficiently small positive numbers T the operator gT can always be approximated by an operator gTN with prescribed accuracy s. However, note that the replacement of gT by gTN is justified only if the sequence

(gTN... gTNxQ),k = 0,1,... k

lies in the ball Ba (jc0 ). This means that number fi is a parameter that permits one to control the systematic error of calculations.

3. Construct of the minimal sets of the Lorenz system

Consider the Lorenz system

x = a(y-x), y = rx-y-xz, z = xy-bz,

(8)

where a, r and b are positive numbers, which play the role of system parameters.

For simplicity, we restrict considerations to the following case of classical values of parameters:

Set

and

(7 = 10, r = 28, ¿ = 8/3.

= ^ —1,2,..., where c = (x0, j0, ^0). Then, obviously, the relations t

and

~ xo> yi ~ yo ' — ^o

x'2 = 10/y0 -10ùc0 +x0, y'2 = ~ ÍVo + 28öc0, xz'2 = -2.667tZo +Z0 +tx0y0.

Here have

x\ = -5t2x0z0 -55t2y0 + lOiVo +190*2x0-10fit0+x0,

y[ = S.mt3y0z0-5t2y0z0 - 8.889i3x0^0 + 6.833i2x0z0 --tx0z0 -3.333fxayl + 3.333i3x02j>0 -0.5t2x%y0 + 140.5i2j>0 -

-ty0 + y0-l54t2x0 + 2Stx0, zl =-3.333i3x0j;0z0 +3.333i3x02z0 -0.5/2x02z0 +3.555i2z0 --2.667feo + Z0 - 3.333?Vo + Si^o + 96.667t3x0y0 --6.S33t2x0y0 +tx0y0 -93.333i3x„ +14?2x„.

It is very difficult to construct the operator g'N for N> 3 without special computer programs. For this purpose, we have used a specially developed software (see Section 4). We cannot cite the regret g'N with N> 3. For example, the operators y'l0 and z{0 contain about 30 000 terms.

The minimal sets of system (8) were constructed for a wide range of initial conditions. Here, in particular, consider the construction of the minimal set lying in the tw-limit set Q of the solution with the initial condition

x0 = -15.720831, y0= -16.587193, z0=36.091132. (9)

The Figure 1 represents the projection of an arc of the orbit L onto the plane xOy, constructed for T = 0.001 on the basis of points

t t

SNC> S M ■

t

with tremendous accuracy

tt t t

&n—&n c ~ &n+\ ■■■Sn+1 c

<10"

Note that solutions of system (8) are defined and bounded for t > 0. Since <£ is the unique attractor of system (8), it follows that there exists a limit

lim inf \g'q-p |= 0,

for almost all çeM3[1].

Therefore, it is fair to say that the path L is set by the desired minimal set to computation error.

4. Horner's Scheme

As was already mentioned, even in the simplest situation of Lorenz system we require multiple computations of a rather complicated multidimensional polynomial for the local solutions construction. With respect to the efficient implementation of this procedure, using special methods is required. For this purpose we have used the generalized Horner's Scheme integrated into the Math-Cloud system [14, 15].

Generally speaking, the description of any version of the generalized Horner's Scheme is extremely awkward. Thus, we will limit ourselves to the situation of four variables polynomial (see Section 4).

Let's consider polynomial

P(t,x,y,z) = Y,aijHt'xJyk

U,k,I

(10)

where a— are some real numbers and t, x, y and z — are real variables. The generalized Horner's Scheme for the polynomial (10) can be depicted as:

25 20 15 10 5

y 0

-5 -10 -15 -20 -25

Fig. 1. The projection onto the plane xOy of an arc of the orbit L near the attractor <L of system (8)

x

P (t,x,y,z) = (...({bm(x,y,z)t+ôm_! (x,y,z) t + + K-2 {x,y,z)t + ...)t + b0 {x,y,z),

where

b, ix>y>z) = (■ ■ -((c™ {y>z)x+cim_ 1 (y,z))x + +Ci,m-2{y,z))x + ..)x + +cm (y,z)

c¡¡ M = (• ■ ■■ i(dVm («) y+m_! (*)) y + +dj/,m-2(z))y + ---)y + +djf0(z)

(11)

and

d,]k (i) = {■■■({aiJkmz+aijkm_l)z +

+ aijkm-2 +.. •) + +aijk0

(12)

(13)

(14)

It is obvious that (12)—(14) scheme is good at parallelizing. According to simulation experiments, the use of this scheme (in comparison to direct scheme (11)) at the same level of accuracy makes it possible to essentially reduce the number of computations even in the simplest situations.

For the construction of the Lorenz system minimum set precisely (12)—(14) version of the generalized Horner's Scheme have been used. A schema of the computation process, using services from the MathCloud system is performed in the Figure 2. As a result, the time it takes for the construction of the arc described at the Figure 1 falls by tens.

Conclusion

Let's discuss computation aspects of the suggested method implementation. Let's recall that in accordance with the method any solution x( t) of the system (1) is restored in the form of variables t, x0 multidimensional polynomial where x(0) = jc0 and t is time. Thus, the storage of the solution obtained comes down to the storage of multidimensional matrix determining this polynomial and matrix determining analytical entry of the expression for the error. One can perform this storage by going through initially all matrix multiplications, or, if more compact, by remembering a parent matrix and set of rules for further conversions. In any case, this storage parallelizes fairly well.

In case of further work with the solution obtained, one for example may need to calculate a trajectory path at given xQ value at t = rt, t2, . . . , r, points. It is obvious that each evaluation at t. instant is independent of others and these evaluations can be done at independent nodes.

As for the evaluation of the multidimensional polynomial, then in the general case this is an ambitious task. It is reasonable to use variations of the generalized Horner's method to solve the task. There are good grounds to believe that the appropriate procedure will be good at parallelizing.

In relatively simple cases, it is possible to overcome the computational complexity using standard methods and to study complex systems behavior. In the paper, there is an approximate construction of Lorentz system minimum sets as an example of implementation of the results obtained. ■

Initial data

Construction of discrete dynamic system

Í

Horner's Scheme

Construction of polynomials

MathCloud

Construction of minimal sets

Construction of local solutions

J

Fig. 2. Computer process organization in the MathCloud system

References

1. Lorenz E.N. (1963) Deterministic nonperiodic flow. Journal of Atmospheric Sciences, no. 20, pp. 130—141.

2. Vallis G.K. (1988) Conceptual models of El Nffio and the southern oscillation. Journal of Geophysical Research, vol. 98, no. 11, pp. 13979— 13991.

3. Bakasov A.A. (1991) Dynamical model of single-mode laser. I. Regime of stable stationary generation. Theoretical and Mathematical Physics, vol. 89, no. 2, pp. 278-292.

4. Shapovalov V.I., Kablov V.F., Bashmakov V.A., Avakumov V.E. (2004) Sinergeticheskaya model' ustoychivosti sredney firmy [Synergetic model of stability of a medium company]. Sinergetika iproblemy teorii upravleniya [Synergetics and problems of control theory]. Moscow, Fizmatlit, pp. 454-464 (in Russian).

5. Nemytskii V.V., Stepanov V.V. (2004) Kachestvennaya teoriyadifferentsial'nykhuravneniy [Qualitative theory of differential equations]. Moscow, 2004 (in Russian).

6. Afanas'ev A.P., Dzyuba S.M. (2015) Method for contraction minimal sets of dynamical systems. Differential Equations, vol. 51, no. 7, pp. 831-837.

7. Afanas'ev A.P., Dzyuba S.M. (2015) Construction of the minimal sets of differential equations with polynomial right-hand side. Differential Equations, vol. 51, no. 11, pp. 1403-1412.

8. Afanas'ev A.P., Dzyuba S.M. (2015) Method for constructing approximate analytic solutions of differential equations with a polynomial right-hand side. Computational Mathematics and Mathematical Physics, vol. 55, no. 10, pp. 1665-1673.

9. Hall G., Watt J.M., eds. (1976) Modern numerical methods for ordinary differential equations. Oxford: Clarendon Press.

10. Hairer E., Norsett S.P., Wanner G. (1987) Solving ordinary differential equations. I: Nonstiff problems. Berlin: Springer-Verlag.

11. Fedorenko R.P. (2008) Vvedenie v vychislitel'nuyufiziku [Introduction to computational physics]. Dolgoprudnyi: Intellect, 2008 (in Russian).

12. Demidovich B.P., Maron I.A., Shuvalova E.Z. (2010) Chislennye metody analiza. Priblizhenie funkttsiy, differentsial'nye i integral'nye uravneni-ya [Numerical methods for analysis: Approximation offunctions, differential and integral equations. St. Petersburg: Lan', 2010 (in Russian).

13. Bakhvalov N.S., Zhidkov N.P., Kobelkov G.M. Chislennye metody [Numerical methods]. Moscow: Binom. Laboratory of Knowledge, 2011 (in Russian).

14. Suhoroslov O.V. (2009) Unifitsirovannyy interfeys dostupa k algoritmicheskim servisam v Web [The unified interface to access algorithmic services in Web]. Problemy vychisleniy v raspredelennoy srede. Trudy ISA RAN [Problems of computations in a distributed environment. Transactions of ISA RAS], vol. 46, pp. 60-82 (in Russian).

15. Afanas'ev A.P., Dzyuba S.M., Emelyanova I.I. (2015) Analytical and numerical investigation for the problem of optimal control of nonlinear system class via quadratic criteria. Procedia Computer Science, no. 66, pp. 23-32.

Схема Горнера для исследования решений дифференциальных уравнений

у у О

с полиномиальной правой частью2

А.П. Афанасьев

доктор физико-математических наук, профессор заведующий лабораторией распределенных вычислительных систем Институт проблем передачи информации им. А.А. Харкевича РАН; профессор, заведующий базовой кафедрой высокопроизводительных вычислений Института проблем передачи информации им. А.А. Харкевича РАН Адрес: 127051, г. Москва, Большой Каретный пер., д. 19 стр. 1 E-mail: apa@isa.ru

С.М. Дзюба

доктор физико-математических наук, профессор кафедры информационных систем Тверской государственный технический университет Адрес: 170026, г. Тверь, наб. Афанасия Никитина, д. 22 E-mail: sdzyuba@mail.ru

И.И. Емельянова

старший преподаватель кафедры информационных систем Тверской государственный технический университет Адрес: 170026, г. Тверь, наб. Афанасия Никитина, д. 22 E-mail: emelyanova-123@yandex.ru

2 Работа выполнена при поддержке Российского научного фонда (проект 16-11-10352)

МАТЕМАТИЧЕСКИЕ МЕТОДЫ И АЛГОРИТМЫ БИЗНЕС-ИНФОРМАТИКИ

Аннотация

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

Ключевые слова: дифференциальные уравнения с полиномиальной правой частью, схема Горнера, интерактивная система вдоль решений исходной системы, система MathCloud.

Цитирование: Afanas'ev A.P., Dzyuba S.M., Emelyanova I.I. Horner's Scheme for investigation of solutions of differential equations with polynomial right-hand side // Business Informatics. 2017. No. 2 (40). P. 33—39. DOI: 10.17323/1998-0663.2017.2.33.39.

Литература

1. Lorenz E.N. Deterministic nonperiodic flow // Journal ofAtmospheric Sciences. 1963. No. 20. P. 130—141.

2. Vallis G.K. Conceptual models of El Nuio and the southern oscillation // Journal of Geophysical Research. 1988. Vol. 98. No. 11. P. 13979— 13991.

3. Bakasov A.A. Dynamical model of single-mode laser. I. Regime of stable stationary generation // Theoretical and Mathematical Physics. 1991. Vol. 89. No. 2. P. 278-292.

4. Шаповалов В.И., Каблов В.Ф., Башмаков В.А., Авакумов В.Е. Синергетическая модель устойчивости средней фирмы // Синергетика и проблемы теории управления / Под ред. А.А. Колесникова. — М.: Физматлит, 2004. С. 454-464.

5. Немыцкий В.В., Степанов В.В. Качественная теория дифференциальных уравнений. М., 2004.

6. Afanas'ev A.P., Dzyuba S.M. Method for contraction minimal sets of dynamical systems // Differential Equations. 2015. Vol. 51. No. 7. P. 831—837.

7. Afanas'ev A.P., Dzyuba S.M. Construction of the minimal sets of differential equations with polynomial right-hand side // Differential Equations. 2015. Vol. 51. No. 11. P. 1403—1412.

8. Afanas'ev A.P., Dzyuba S.M. Method for constructing approximate analytic solutions of differential equations with a polynomial right-hand side // Computational Mathematics and Mathematical Physics. 2015. Vol. 55. No. 10. P. 1665—1673.

9. Hall G., Watt J.M., eds. Modern numerical methods for ordinary differential equations. Oxford: Clarendon Press, 1976.

10. Hairer E., Norsett S.P., Wanner G. Solving ordinary differential equations. I: Nonstiff problems. Berlin: Springer-Verlag, 1987.

11. Федоренко Р.П. Введение в вычислительную физику. Долгопрудный: Интеллект, 2008.

12. Демидович Б.П., Марон И.А., Шувалова Е.З. Численные методы анализа. Приближение функций, дифференциальные и интегральные уравнения. СПб: Лань, 2010.

13. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Бином. Лаборатория знаний, 2011.

14. Сухорослов О.В. Унифицированный интерфейс доступа к алгоритмическим сервисам в Web // Проблемы вычислений в распределенной среде. Труды ИСА РАН / Под ред. С.В. Емельянова, А.П. Афанасьева. М., 2009. Т. 46. С. 60—82.

15. Afanas'ev A.P., Dzyuba S.M., Emelyanova I.I. (2015) Analytical and numerical investigation for the problem of optimal control of nonlinear system class via quadratic criteria // Procedia Computer Science. 2015. No. 66. P. 23—32.

БИЗНЕС-ИНФОРМАТИКА № 2(40) — 2017