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

Прочитав статью «Динамическая система Лоренца и вычислительный эксперимент», проверил расчеты с помощью аналитически-численного метода [1].

Результаты расчета на фазовой плоскости z(x):


И y(x):


Кажется, что кривые замкнуты, но давайте рассмотрим результат поподробнее.

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

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



где — приближенное решение (i-я фазовая координата);
— неизвестное точное решение;
— верхняя оценка предельной полной погрешности расчета приближенного решения;


Взяв параметры для расчета из статьи «Динамическая система Лоренца и вычислительный эксперимент»:
Предначальные условия, параметры динамической системы, точность математических операций — 180 знаков после запятой, точность по степенному ряду 1e-9, получим следующий результат в точке t = 6.827:





Значения производных:



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



Отмечу, что повышение точности расчетов (количество учитываемых знаков после запятой и точность по степенному ряду) приводит лишь к сужению области, содержащей точные решения. Например, при задании точности 1e-55, область в точке t = 6.827 сужается до .

Далее, я решил продолжить расчет до точки t = 12.827 и рассмотреть график результатов расчета на фазовых плоскостях z(x):



И y(x):



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

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

Литература:
1. Бычков Ю., Щербаков С. Аналитически-численный метод расчета динамических систем. — Санкт-Петербург: Энергоатомиздат,2001.

Similar posts

Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 31

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

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

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

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

        Исходный код моей программы занимает несколько тысяч строк кода (зачастую без комментариев) и с ее помощью можно рассчитывать различные нелинейные ОДУ, не только эту систему.
          0
          Сейчас проверяю Ваш результат. Я взял в своей программе тип не mpreal, а double, и провел вычисления ещё раз. Получается следующий результат:

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

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

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

          Дайте ссылку на статью, интересно было бы познакомиться с ней.
            0
            В программе используется библиотека символьных вычислений GiNaC.
            Исходный код ядра программы, в котором реализованы формализованные процедуры аналитически-численного метода, залил на гитхаб: anm
            Преобразование в double используется только при построении графиков. При расчетах нигде не используется.
            Статей, опубликованных в онлайн по этой теме нет.
              0
              Значит, может быть неявное преобразование. Чтобы проверить результат, попробуйте пройти назад, как делал я. Если там действительно где-то есть конверт в double, то буквально через несколько шагов из-за сильной неустойчивости приближенное решение унесёт в бесконечность.
                0
                Попробовал пройти назад, как делали Вы. Расчет прошел успешно от t = 6.827 до t = 0.
                Что примечательно, если в качестве предначальных условий взять округленные до 6 знака после запятой значения, полученные при прямом расчете, то x(t) действительно становится неустойчивым в районе t = 5.4.
                Но, если задавать предначальные условия не округленными, а точно такими, какими они были получены, расчет идет нормально до t = 0.
                  0
                  Сморите, я поправил свой код на постоянный шаг 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 (старый вариант программы). Так вот во втором случае условие выхода из цикла стояло внутри тела, только сейчас заметил отличие.

                  Итак, по сути, наши результаты не противоречат друг другу.
                    0
                    Рад что удалось разобраться в чем была причина расхождения результатов.
                      0
                      Я немного поправил свой топик (внесены изменения в код, параметры и результаты вычислений).
                    +1
                    Можно код доделать до того, чтобы доходя до конечной точки, пересчитывался шаг как разность конечного момента времени и предыдущего, но это, как я считаю, мелочи.
                      0
                      Забыл ещё поправить: там при выводе надо писать не
                      cout << "\nКоординаты в конечный момент времени "  << t << ...
                      

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

                      поскольку момент времени будет увеличен на величину шага.
                  +1
                  Забыл предупредить — в исходниках, в комментариях, может присутствовать ненормативная лексика.
                    0
                    Провел расчеты в 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...]
                      0
                      Детали вычислений скрыты в пакете, что там внутри делается, мы не знаем. Многие мои знакомые математики пакетам не доверяют, сами пишут программы, где нужен контроль численного метода (и погрешности получаемого результата).
                        0
                        Полностью согласен. И когда плотно занимался анализом даже находил примеры ОДУ, для которых пакеты находили некорректное решение.
                        Но для проверки собственных вычислений нужно же чем-то пользоваться.
                +2
                Quiz, спасибо, что дали знать, что появилась такая статья.
                +5
                То, что аттрактор Лоренца не может стать переодическим, следует из результатов полученных еще Андроновым задолго до изобретения компьютеров и самого аттрактора. Автор верно указал на ошибку, но при наличии аналитического доказательства его следует предпочесть численному эксперименту.

                Применяется для решения обыкновенных нелинейных неавтономных нестационарных интегродифференциальных уравнений
                Просто с ходу не нашел книгу, указанную в источниках литературы, а есть ли обоснование для такого усложнения? Интегродифференциальные уровнения достаточно сложны для решения, не мог бы автор привести формулу после преобразования?
                  +1
                  К моему великому сожалению авторы книги не пожелали выкладывать ее в онлайн — опасаются нарушения своих прав, видимо. Поэтому ее можно найти только в аналоговых библиотеках.
                  Обоснование усложнения достаточно простое — используемый инструмент (ряды Тейлора, преобразования Лапласа, ряд Лорана) позволяет решать описанный класс задач. Он достаточно широк, но в плане используемого аппарата это один класс задач и один подход.
                  Столь точная формулировка необходима чтобы ограничить класс задач.
                  Какую именно формулу Вы хотите увидеть?
                    0
                    Просто хочу подробнее узнать про ваш численный метод: похож ли он на явные методы Рунге — Кутты, неявный, если свелся к интегродифференциальным уравнениям, то с какой точностью вычисляются интегралы на каждом шаге, какой у него порядок сходимости? Может вы знакомы с этими изданиями:
                    Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи
                    Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи
                    По-моему, в них описаны все лучшие на данный момент численные методы для ОДУ, было бы интересно узнать, как называется предложенный метод в ихней формулировке.
                      0
                      Поддчеркну, что это не мой метод. Над ним работали и сейчас продолжают работать несколько ученых. Я лишь осуществил программную реализацию метода.
                      С указанными изданиями знаком поверхностно. Как оценить «похожесть» на методы Рунге-Кутты?
                      Например, аналитически-численный метод (он так и называется) снабжён процедурой выбора величины шага расчёта, которая принципиально отличается от всех других существующих методов. Процедура основана на анализе сходимости числовых мажорант по признаку Вейерштрасса.
                        0
                        Вроде я понял. Вы не делаете программу, которой просто скармливают систему ОДУ и она выдает частные решения, как это делают все пакеты, а предварительно разбираете систему аналитически. Подход не универсальный, но интересный, с радостью бы посмотрел на результаты и примеры, когда он дает более высокую точность или скорость.
                          +2
                          Совершенно верно. Метод включает в себя аналитическую часть и численную. Причем аналитическая часть выполняется единожды, а численная — на каждом шаге расчета.
                          Если сообществу интересно могу описать процедуру расчета на примере анализа системы Лоренца — это будет уже отдельная статья, поскольку описание слишком велико для комментария.
                • UFO just landed and posted this here
                    –1
                    Книга доступна во всех крупных библиотеках России. А также в библиотеке СПбГУ «ЛЭТИ» и ряде других вузов.
                    • UFO just landed and posted this here
                      +1
                      А я рекомендую указанную выше пользователем paunch книгу Хайрер Э., Нёрсетт С. и Ваннер Г. «Решение обыкновенных дифференциальных уравнений. Нежесткие задачи». Это классика. Стр. 49-53, там же есть ссылки на зарубежные работы.
                        0
                        Эта книга есть в сети.
                          0
                          Также можно посмотреть и другую классику: Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Стр. 365-366.
                            0
                            Вот еще статья.
                            • UFO just landed and posted this here

                          Only users with full accounts can post comments. Log in, please.