Lorentz dynamic system and computational experiment

This post is a continuation of my article [1] on Habrahabr about the Lorenz attractor. Here we consider a method for constructing approximate solutions of the corresponding system, paying attention to the software implementation.

The dynamic Lorentz system is an autonomous system of ordinary differential equations of the third order



where , r and b are positive numbers. To study the behavior of system solutions, one usually takes the classical values ​​of the system parameters, i.e.



As noted in [1], in this case, instability of its solutions on the attractor takes place in system (1). In fact, this makes it incorrect to apply classical numerical methods over large periods of time (and the limit sets of dynamical systems are built on such segments). One of the ways to overcome this problem is the transition to high-precision computing, but this approach puts the researcher in a tight framework: firstly, a small degree of freedom to reduce errors (change in the step of integration and the accuracy of the representation of a real number to control the computing process), and second , a large amount of computation at very small .

Another solution to this problem may be the use of the power series method. In [2], a modification of this method for dynamic systems of the form (1) is described, which allows one to quickly determine the expansion coefficients



as



where the coefficients are and are determined by the initial conditions; the length (step) of the segment of convergence of the series (2) is also obtained :



if , then , otherwise ,

is any positive number.

Note that the scheme given in [2] for obtaining the length of the convergence segment of power series can be transferred by analogy to other third-order dynamical systems with nonlinearities of the form (1) (for example, a system describing the behavior of a self-developing market economy [3, p. 261]) .

Despite the fact that all trajectories of system (1) are bounded and its right-hand side is everywhere analytic, an initial computational experiment showed that the radius of convergence of series (2) is bounded and depends on the choice of initial conditions. Therefore, in the described way, we can get only a part of the trajectory. The procedure for constructing a trajectory arc on any time interval consists in stitching together the parts of the trajectory that make up the desired solution, on which the series (2) converge. The integration error accumulated during the transition from arc to arc of the trajectory due to the error in finding the current approximate solution can be controlled by varying the accuracy of expansion in a power series. The use of high-precision calculations allows one to continue the solutions for very large periods of time, sinceimpossible to make arbitrarily small due to the restriction on machine epsilon.

We calculate by formulas (3) , and up to a value of i at which the estimate takes place.



It is clear that we can’t avoid a significant accumulation of integration errors over long time periods, therefore we will implement high-precision floating-point calculations based on the GNU MPFR Library , or rather, C ++ - the MPFR library interface for working with real numbers of arbitrary precision - MPFR C ++ . It is convenient in that it has the mpreal class with overloaded arithmetic operations and friendly mathematical functions. To install the library on Ubuntu Linux, do

sudo apt-get install libmpfrc++-dev

In addition to the package libmpfr-dev, package managers will also be affected by libgmp-dev. This is the devel package of the GNU Multiple Precision Arithmetic Library or GMP (MPFR is an extension of it).

Consider an example C ++ code for computing the values ​​of phase coordinates at a finite moment in time, as well as checking the values ​​found.

#include 
#include 
using namespace std;
#include 
using namespace mpfr;
// Точность по степенному ряду
#define eps_c   "1e-50"
// Параметры системы уравнений Лоренца
#define sigma   10
#define r       28
#define b       8/(mpreal)3
// Количество бит под мантиссу вещественного числа
#define prec    180
// Как считать шаг: 1 - по оценке отрезка сходимости, 0 - заданная величина
#define FL_CALC 1
// Фиксированный шаг по времени
#define step_t  "0.02"
// Функция возвращает длину отрезка сходимости степенного ряда
mpreal get_delta_t(mpreal &alpha0, mpreal &beta0, mpreal &gamma0)
{
    mpreal h2 = (mpfr::max)((mpfr::max)(fabs(alpha0), fabs(beta0)), fabs(gamma0)),
           h1 = (mpfr::max)((mpfr::max)(2*sigma, r+2*h2+1), b+2*h2+1);
    mpreal h3 = h2 >= 1 ? h1 * h2 : (mpfr::max)((mpfr::max)(2*sigma, r+2), b+1);
    return 1/(h3 + "1e-10");
}
// Функция вычисления значений фазовых координат в конечный момент времени
// x, y и z - координаты начальной точки; T - длина отрезка интегрирования;
// way - направление поиска решений: 1 - вперед по времени, -1 - назад по времени
void calc(mpreal &x, mpreal &y, mpreal &z, mpreal T, int way = 1)
{
    cout << "\nВ начальный момент времени:\n\ndx/dt = " << sigma*(y-x) << "\ndy/dt = " <<
            r*x-y-x*z << "\ndz/dt = " << x*y-b*z << endl;
    mpreal t = 0, delta_t, L, p, s1, s2;
    bool fl_rp;
    do
    {
        if(FL_CALC)
            delta_t = get_delta_t(x, y, z);
        else
            delta_t = step_t;
        t += delta_t;
        if(t < T)
            fl_rp = true;
        else if(t > T)
        {
            delta_t -= t-T;
            fl_rp = false;
        }
        else
            fl_rp = false;
        vector alpha, beta, gamma;
        alpha.push_back(x);
        beta.push_back(y);
        gamma.push_back(z);
        int i = 0;
        L = sqrt(alpha[0]*alpha[0] + beta[0]*beta[0] + gamma[0]*gamma[0]);
        p = way * delta_t;
        while(L > eps_c)
        {
            // Вычисляем новые коэффициенты степенных рядов
            s1 = s2 = 0;
            for(int j = 0; j <= i; j++)
            {
                s1 += alpha[j] * gamma[i-j];
                s2 += alpha[j] * beta[i-j];
            }
            alpha.push_back(sigma*(beta[i] - alpha[i])/(i+1));
            beta.push_back((r*alpha[i] - beta[i] - s1)/(i+1));
            gamma.push_back((s2 - b*gamma[i])/(i+1));
            i++;
            x += alpha[i] * p;
            y += beta[i] * p;
            z += gamma[i] * p;
            L = fabs(p) * sqrt(alpha[i]*alpha[i] + beta[i]*beta[i] + gamma[i]*gamma[i]);
            p *= way * delta_t;
        }
    }
    while(fl_rp);
    cout << "\nКоординаты в конечный момент времени:\nx = " << x.toString() << "\ny = " <<
            y.toString() << "\nz = " << z.toString() << endl;
    cout << "\nЗначения производных:\n\ndx/dt = " << sigma*(y-x) << "\ndy/dt = " <<
            r*x-y-x*z << "\ndz/dt = " << x*y-b*z << endl;
}
int main()
{
    mpreal::set_default_prec(prec);
    cout << "Машинный эпсилон = " << machine_epsilon() << endl;
    mpreal T;
    cout << "\nВведите длину отрезка времени > ";
    cin >> T;
    mpreal x, y, z;
    cout << "\nx0 > ";
    cin >> x;
    cout << "y0 > ";
    cin >> y;
    cout << "z0 > ";
    cin >> z;
    cout << endl;
    calc(x, y, z, T);
    cout << "\n\n*** Проход назад ***\n";
    calc(x, y, z, T, -1);
    return 0;
}

To compile the lorenz.cpp file with this code, do the following

g++ lorenz.cpp -std=gnu++11 -lmpfr -lgmp -o lorenz

As can be seen from the program listing, vectors from the STL library are used to store the values ​​of the coefficients of power series. It was assumed that the number of bits under the mantissa of a real number is 180. The range of the exponent is the default from to (on 32-bit platforms). The machine_epsilon () function of machine epsilon definition was used to find out its value for a given representation of a real number. In this case, we present the



results of a computational experiment. As initial conditions, we take values ​​close to the Lorentz attractor:



The length of time at which the calculations were performed is - . Initially, the accuracy of the estimation of the general term of the series was taken equal . This is enough to get the result.

;

derivatives:



(decreasing the value , the coordinate values ​​will not change). In this case, you can take real numbers with less memory.

The values ​​of the derivatives, as well as the values ​​of the coordinates, are given in order to illustrate the fact of the return of the trajectory to the vicinity of the starting point, but non-closure, as one would expect from the hypothesis of the existence of cycles in system (1). After such a rapprochement, the point of the trajectory departs from its initial position, but then again returns to its vicinity. Such behavior is predicted by the qualitative theory of differential equations (Poisson stability of points of recurrent trajectories on the attractor; this was indicated in [1] and discussed at the conference [4]).

The indicated bit depth for a real number was chosen to track not only the return of the trajectory of the Lorentz system to the vicinity of the initial conditions [1, Fig. 1], but also for going backward in time from the end point to the initial trajectory along the arc (equality of the way parameter to calc () -1). Then in the calculations you need to take . The result values coinciding with the initial to the seventh decimal place (Method toString () class mpreal, called without parameter indicates the conversion of a string with a maximum number of decimal places):

x = 13.412656286837273085165416945301946328440634370684244
y = 13.4643000297481126631507883918720904312092673686014399
z = 33.4615641630148784946354299167181879731357599130041067

This small valueIt was chosen because when moving along the trajectory to negative values ​​of time, a strong instability of the solutions is observed - they immediately go to infinity from the attractor, since in the calculations we are close to it, and not directly on it.

The program also provides a fixed step value. It is verified that the number 0.02 can be used to construct approximate solutions near the Lorentz attractor. This value is much larger than that obtained from the above estimate from [2] (FL_CALC flag is 1) for any initial conditions, but when the starting point is removed a considerable distance from the attractor, the method ceases to work (the rows do not converge).

In [5, p. 90, 91] to study the trajectories of system (1), the Euler method with a variable step is usedintegration. The value is selected with error tracking at the current step (i.e., local control) using interval arithmetic. However, there is no verification of a common integration error. Passing back in time, described above, guarantees the correct construction of an approximate solution. By varying the accuracy of the estimation of the general term of the series, we can reduce the error at the step, which does not allow the Euler method. Also, in the modification of the power series method considered above, the advantage over the general scheme of the Taylor series method is the quick calculation of the expansion coefficients using formulas (3) compared to the procedure of symbolic differentiation of the right-hand sides of system equations (in addition, in the nonlinear case, symbolic expressions require a lot of memory to store them in the calculation of derivatives of higher orders).

Thus, knowing the state of the Lorentz system in the past, with a sufficient degree of accuracy, we can predict the behavior of its trajectories over long time intervals, and also go back. In fact, the formal definition of chaos is violated here [6, p. 118, 119].

Literature

1 . Pchelintsev A.N. A critical look at the Lorenz attractor, 2013. Habrahabr .
2 . Pchelintsev AN Numerical and physical modeling of the dynamics of the Lorenz system // Numerical Analysis and Applications, 7 (2): 159-167, 2014. DOI: 10.1134 / S1995423914020098 .
3 . Magnitsky N.A., Sidorov S.V. New methods of chaotic dynamics. - M .: Editorial URSS, 2004.
4 . Pchelintsev A.N. On modeling the dynamics of the Lorentz system // International Conference “Kolmogorov Readings VI. General management problems and their applications (OPU-2013) ”, Tambov State University named after G.R. Derzhavina (Tambov, 2013). Math-Net.Ru .
5 . Tucker W.Rigorous ODE Solver A and Smale's 14Th Problem // Foundations of Computational Mathematics, 2 (1): 53-117, 2002.
6 . Berger P., Pomo I., Vidal K. Order in chaos. About the deterministic approach to turbulence. - M.: Mir, 1991.