Гарантии получения корректного результата при расчете динамических систем

#include <iostream>
#include <vector>
using namespace std;

#include <mpreal.h>
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 0

// Фиксированный шаг по времени
#define step_t  "1e-3"

// Функция возвращает длину отрезка сходимости степенного ряда
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Координаты в начальный момент времени:\nx = " << x.toString() << "\ny = " << y.toString() <<
            "\nz = " << z.toString() << endl;
    mpreal t = 0;
    mpreal delta_t;
    while (1)
    {
        if(FL_CALC)
            delta_t = get_delta_t(x, y, z);
        else
            delta_t = step_t;

        t += delta_t;
	if(t > T)
	    break;

        vector<mpreal> alpha, beta, gamma;
        alpha.push_back(x);
        beta.push_back(y);
        gamma.push_back(z);

        int i = 0;
        mpreal L = sqrt(alpha[0]*alpha[0] + beta[0]*beta[0] +
                        gamma[0]*gamma[0]);
        mpreal p = way * delta_t;
        while(L > eps_c)
        {
            // Вычисляем новые коэффициенты степенных рядов
            mpreal s1 = 0, 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;
        }
    }

    cout << "\nКоординаты в конечный момент времени "  << t << ":\nx = " << x.toString() << "\ny = " << y.toString() << "\nz = " << z.toString() << 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;
}

cout << "\nКоординаты в конечный момент времени "  << t << ...

cout << "\nКоординаты в конечный момент времени "  << t-delta_t << ...