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

[Quiz]

Ждём комментария pchelintsev_an

[pchelintsev_an]

Таким образом, нельзя говорить ни о каком возвращении траектории в окрестность начальной точки — об этом говорится в статье. А делать выводы на основе расчетов необходимо всегда с оглядкой на погрешность вычислений (как методическую так и вычислительную).

Не я придумал, что предельная траектория возвращается в окрестность своего начального положения. Она является рекуррентной => любая точка на ней устойчива по Пуассону, в том числе и начальная. Это говорит о том, что траектория должна возвращаться в любую окрестность данной точки. Если моменты времени, которые этому соответствуют, образуют относительно плотное множество, то траектория является почти периодической. Всю эту «радость» можно прочитать в классических работах по качественной теории дифференциальных уравнений (например, книга В.В. Немыцкого и В.В. Степанова).
Теперь о радиусе окрестности. Если посмотрите числа, которые я получил в статье, то можно определить, что расстояние между начальной и конечной точкой на дуге траектории, взятое по Евклиду, будет меньше величины 0.13. Я не понимаю, зачем Вы привели проекции дуги траектории (они, на мой взгляд, не информативны). Вот, смотрите, что я имел в виду, когда говорил о возвращении траектории.

Если ее проматывать дальше, то должен найтись другой момент времени, когда мы вернёмся в окрестность меньшего радиуса. Еще раз — так говорит нам теория. Как видите, радиус окрестности и величина ошибки у меня отличаются на несколько порядков.
Не понимаю, почему у нас отличаются результаты вычислительного эксперимента. Как конкретно считали? Какой брали шаг, или он был переменным (как у меня в программе), ничего об этом Вы не написали, только кинули фразу

Несложно видеть, что результаты расчетов несколько отличаются от изложенных в статье.

и показали числа. Хорошо бы еще Вашу программу посмотреть.

[Surgun]

По поводу возвращения в окрестность начального положения — вероятно, я не прав. С работой В.В. Немыцкого и В.В. Степанова не знаком, попробую почитать.
Вот результаты моего расчета, пошагово: «Лоренц»
Давайте попробуем сравнить с Вашими аналогичными результатами.

Процедура выбора шага расчета основана на анализе сходимости числовых мажорант по признаку Вейерштрасса. Подробное описания процедуры поиска шага расчета выходит за рамки комментария.
По поводу результатов вычислительного эксперимента — попробуйте подставить результат Вашего расчета в исходную систему ОДУ — численные значения производных немного отличаются от представленных Вами в статье. Или я где-то ошибаюсь?

Исходный код моей программы занимает несколько тысяч строк кода (зачастую без комментариев) и с ее помощью можно рассчитывать различные нелинейные ОДУ, не только эту систему.

[pchelintsev_an]

Сейчас проверяю Ваш результат. Я взял в своей программе тип не mpreal, а double, и провел вычисления ещё раз. Получается следующий результат:

x(6.827) = 13.3595, y(6.827) = 13.3786, z(6.827) = 33.4297,

то есть Ваш результат до четвёртого знака после запятой. Это наводит на мысль о том, что внутри Ваша программа где-то преобразует промежуточный результат в double, тем самым начинает собираться ошибка. Проверяйте свой код.
Кстати, если Ваш код работает с любой системой, имеющей голоморфные решения, то, скорее всего, Вы используете символьные вычисления.

Исходный код моей программы занимает несколько тысяч строк кода (зачастую без комментариев) и с ее помощью можно рассчитывать различные нелинейные ОДУ, не только эту систему.

А почему бы не привести либо фрагменты, либо описать алгоритм. Я бы еще залил сырцы куда-нибудь, секретного в них ничего нет.

Подробное описания процедуры поиска шага расчета выходит за рамки комментария.

Дайте ссылку на статью, интересно было бы познакомиться с ней.

[Surgun]

В программе используется библиотека символьных вычислений GiNaC.
Исходный код ядра программы, в котором реализованы формализованные процедуры аналитически-численного метода, залил на гитхаб: anm
Преобразование в double используется только при построении графиков. При расчетах нигде не используется.
Статей, опубликованных в онлайн по этой теме нет.

[pchelintsev_an]

Значит, может быть неявное преобразование. Чтобы проверить результат, попробуйте пройти назад, как делал я. Если там действительно где-то есть конверт в double, то буквально через несколько шагов из-за сильной неустойчивости приближенное решение унесёт в бесконечность.

[Surgun]

Попробовал пройти назад, как делали Вы. Расчет прошел успешно от t = 6.827 до t = 0.
Что примечательно, если в качестве предначальных условий взять округленные до 6 знака после запятой значения, полученные при прямом расчете, то x(t) действительно становится неустойчивым в районе t = 5.4.
Но, если задавать предначальные условия не округленными, а точно такими, какими они были получены, расчет идет нормально до t = 0.

[pchelintsev_an]

Сморите, я поправил свой код на постоянный шаг 0.001

#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;
}

Отличие от того, что приведен в моей статье, в том, что я раньше выходил из цикла, когда момент времени, в который были произведены расчеты, был больше 6.827, т.е. на один шаг (не на очень большую величину) пересчитывал, чем Вы. Поэтому у нас и сходилось до первого знака после запятой.
Как видите, в данном случае я условие занёс в тело цикла, но здесь при переменном шаге (или шаге, не укладывающимся в отрезке интегрирования целое число раз) возможно то, что мы не пересчитаем (как было раньше), а не досчитаем.
Теперь результат этого кода:

x = 13.3595190692347808097837326351375661843913136556917914
y = 13.3786456324828156454170338272219075410817124310292579
z = 33.4297498991624334573543250465721240701422078412407597

Как видите, хорошо сошлось с Вашими результатами.
Теперь, почему с double выдавало тоже самое. У меня на компе две проги: одна с mpreal, а другая c double (старый вариант программы). Так вот во втором случае условие выхода из цикла стояло внутри тела, только сейчас заметил отличие.

Итак, по сути, наши результаты не противоречат друг другу.

[Surgun]

Рад что удалось разобраться в чем была причина расхождения результатов.

[pchelintsev_an]

Я немного поправил свой топик (внесены изменения в код, параметры и результаты вычислений).

[pchelintsev_an]

Можно код доделать до того, чтобы доходя до конечной точки, пересчитывался шаг как разность конечного момента времени и предыдущего, но это, как я считаю, мелочи.

[pchelintsev_an]

Забыл ещё поправить: там при выводе надо писать не

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

а

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

поскольку момент времени будет увеличен на величину шага.

[Surgun]

Забыл предупредить — в исходниках, в комментариях, может присутствовать ненормативная лексика.

[Surgun]

Провел расчеты в Maple 17. Результат аналогичен моему:
Digits := 180
l1 := diff(x(t), t) = sigma*(y(t)-x(t));
l2 := diff(y(t), t) = ro*x(t)-y(t)-x(t)*z(t);
l3 := diff(z(t), t) = x(t)*y(t)-beta*z(t);
sigma := 10;
ro := 28;
beta := 8/3;
ics := x(0) = 13.41265629, y(0) = 13.46430003, z(0) = 33.46156416;
sol := dsolve([l1, l2, l3, ics], [x(t), y(t), z(t)], numeric, method = taylorseries)
sol(6.827);
[t = 6.827,
x(t) = 13.359519...,
y(t) = 13.3786456...,
z(t) = 33.4297499...]

[pchelintsev_an]

Детали вычислений скрыты в пакете, что там внутри делается, мы не знаем. Многие мои знакомые математики пакетам не доверяют, сами пишут программы, где нужен контроль численного метода (и погрешности получаемого результата).

[Surgun]

Полностью согласен. И когда плотно занимался анализом даже находил примеры ОДУ, для которых пакеты находили некорректное решение.
Но для проверки собственных вычислений нужно же чем-то пользоваться.

[pchelintsev_an]

Quiz, спасибо, что дали знать, что появилась такая статья.

[paunch]

То, что аттрактор Лоренца не может стать переодическим, следует из результатов полученных еще Андроновым задолго до изобретения компьютеров и самого аттрактора. Автор верно указал на ошибку, но при наличии аналитического доказательства его следует предпочесть численному эксперименту.

Применяется для решения обыкновенных нелинейных неавтономных нестационарных интегродифференциальных уравнений

Просто с ходу не нашел книгу, указанную в источниках литературы, а есть ли обоснование для такого усложнения? Интегродифференциальные уровнения достаточно сложны для решения, не мог бы автор привести формулу после преобразования?

[Surgun]

К моему великому сожалению авторы книги не пожелали выкладывать ее в онлайн — опасаются нарушения своих прав, видимо. Поэтому ее можно найти только в аналоговых библиотеках.
Обоснование усложнения достаточно простое — используемый инструмент (ряды Тейлора, преобразования Лапласа, ряд Лорана) позволяет решать описанный класс задач. Он достаточно широк, но в плане используемого аппарата это один класс задач и один подход.
Столь точная формулировка необходима чтобы ограничить класс задач.
Какую именно формулу Вы хотите увидеть?

[paunch]

Просто хочу подробнее узнать про ваш численный метод: похож ли он на явные методы Рунге — Кутты, неявный, если свелся к интегродифференциальным уравнениям, то с какой точностью вычисляются интегралы на каждом шаге, какой у него порядок сходимости? Может вы знакомы с этими изданиями:
Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи
Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи
По-моему, в них описаны все лучшие на данный момент численные методы для ОДУ, было бы интересно узнать, как называется предложенный метод в ихней формулировке.

[Surgun]

Поддчеркну, что это не мой метод. Над ним работали и сейчас продолжают работать несколько ученых. Я лишь осуществил программную реализацию метода.
С указанными изданиями знаком поверхностно. Как оценить «похожесть» на методы Рунге-Кутты?
Например, аналитически-численный метод (он так и называется) снабжён процедурой выбора величины шага расчёта, которая принципиально отличается от всех других существующих методов. Процедура основана на анализе сходимости числовых мажорант по признаку Вейерштрасса.

[paunch]

Вроде я понял. Вы не делаете программу, которой просто скармливают систему ОДУ и она выдает частные решения, как это делают все пакеты, а предварительно разбираете систему аналитически. Подход не универсальный, но интересный, с радостью бы посмотрел на результаты и примеры, когда он дает более высокую точность или скорость.

[Surgun]

Совершенно верно. Метод включает в себя аналитическую часть и численную. Причем аналитическая часть выполняется единожды, а численная — на каждом шаге расчета.
Если сообществу интересно могу описать процедуру расчета на примере анализа системы Лоренца — это будет уже отдельная статья, поскольку описание слишком велико для комментария.

[Surgun]

Книга доступна во всех крупных библиотеках России. А также в библиотеке СПбГУ «ЛЭТИ» и ряде других вузов.

[pchelintsev_an]

А я рекомендую указанную выше пользователем paunch книгу Хайрер Э., Нёрсетт С. и Ваннер Г. «Решение обыкновенных дифференциальных уравнений. Нежесткие задачи». Это классика. Стр. 49-53, там же есть ссылки на зарубежные работы.

[pchelintsev_an]

Эта книга есть в сети.

[pchelintsev_an]

Также можно посмотреть и другую классику: Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Стр. 365-366.

[pchelintsev_an]

Вот еще статья.