Как стать автором
Обновить

Комментарии 12

Не увидел автоматического выбора шага интегрирования. Собственно, возможность этой процедуры обеспечивается одновременным получением решений разного порядка точности и является основным преимуществом схем методов типа Рунге-Кутта подобного рода.
В приведенных модельных задачах E=sum((1/30)*k1-(128/4275)*k3-(2197/75240)*k4+(1/50)*k5+(1/55)*k6) меняется в пределах порядка и шаг остаётся равным заданному.
При необходимости можно добавить
if abs(E)>e:
tau=tau/2
elif abs(E)<e/32:
tau=2*tau
elif e/32<abs(E)<e:
tau=tau
В целом не зависимо от метода это повышает устойчивость!!! (о погрешности понятно не говорим)
За замечание спасибо!!! Приятно встретить практика по численным методам
В случае непродолжаемых решений (как в приведённом примере с уравнением x' = 1 + x^2) РК-методы с постоянным шагом ненадёжены — не зная о вертикальных асимптотах, мы можем просто проскочить их (нужно менять шаг при приближении к ним, при этом контролируя погрешность численного решения). Опять же для нелинейных динамических систем, где имеются странные аттракторы, при численном интегрировании на больших отрезках времени может получиться значительная ошибка по фазовым координатам. Выходом здесь является переход, например, к методу степенных рядов (для систем общего вида с квадратичными нелинейностями метод описан здесь; программа на C++). Вот пример, где РК-метод явно проиграет, — возьмите отсюда начальную точку для третьего цикла в системе Лоренца и проинтегрируйте на периоде (сильно мельчить шагом придётся).
В приведенном модельном примере с du/dt = 1 + u^2 есть точное решение u(t)=tan(t) для сравнительного анализа погрешностей указанных методов в установленном диапазоне. Это корректная постановка задачи. Применяемые численные методы не предполагают «чёрного ящика» цель моделирования всегда известна.Как реализовать переменный шаг опять не проблема. Всё что Вы написали есть справедливым для конкретной задачи.
Спасибо за комментарий.
Простите за вопрос, а зачем тогда численно решать уравнение, когда известно точное решение? Определить погрешность? Для реальных систем точное решение мы можем не знать.
Я сформулировал задачу достаточно точно. Есть и другие способы оценки точности, но я выбрал этот. В ответе речь шла не об известности решения, а об известности задачи, например, из нелинейного дифференциального уравнения движения тела, брошенного под углом к горизонту при сопротивлении пропорциональном квадрату скорости, найти численные значения максимальной высоты подъёма и времени движения. Есть начальные условия, есть представление о форме траектории, нужно найти численные значения времени и максимальной высоты полёта формы траектории в зависимости от сопротивления среды. Эта простая задача, но и более сложные формулируются по такому же принципу. Всегда известно почему конкретное нелинейное ДУ не решается аналитическим методом, из этих знаний можно выбрать метод численного решения.
Всегда известно почему конкретное нелинейное ДУ не решается аналитическим методом, из этих знаний можно выбрать метод численного решения

Разве? Я то думал, что из анализа устойчивости решений (сделать как можно полный качественный анализ системы) и жёсткости системы.
Для одного дифференциального уравнения n – го порядка, задача Коши — СДУ-метод численного анализа -анализ устойчивости решений. У Вас другая последовательность?
1. Продолжаемость решений. Есть теоремы, позволяющие установить нелокальную продолжаемость решений (здесь также можно попробовать построить функцию Ляпунова, см. пункт 2). Если проверить это нельзя, значит, могут быть так называемые сингулярности в решении, т.е. наличие асимптот. Тогда нужно использовать численные методы, которые «чувствуют» асимптоты. Есть такие. Примеры — модифицированный метод Эйлера и упоминаемая работа в предыдущем моём комментарии.
2. Ограниченность решений на оси (мы строим приближённое решение на отрезке, однако, такие сведения тоже нужны). Если решение неорганичено, то осложняется работа с ним в плане хранения вещественных чисел для узлов сетки, а также устойчивости.
3. Если решение ограничено (как в случае с аттракторами), посмотреть, есть ли положения равновесия и их устойчивость по первому приближению. Если есть сёдла, то предельные решения будут неустойчивыми, и применение РК-методов на больших отрезках времени будет некорректным.
4. Если с устойчивостью всё в порядке, то исследуем жёсткость системы, чтобы выбрать явный или неявный метод решения СДУ.
5. В случае нежёсткой задачи подходят РК-методы с автоматическим выбором шага. Вообще, при выборе РК-метода для жёсткой и нежёсткой системы рекомендую воспользоваться двумя монографиями Э. Хайрера и др. «Решение обыкновенных дифференциальных уравнений».
6. Также надо обращать внимание на наличие разрывов и модулей (а также разрывов в производных) в правой части системы по фазовым координатам. Вообще, строго говоря, при переходе через них нарушается теорема существования и единственности решения задачи Коши. В плане численных методов может произойти скачок в погрешности численного решения.
Благодарю Вас за содержательный комментарий приведенного в моём предыдущем ответе четвёртого этапа (в основном применительно к методу Эйлера). Я обязательно воспользуюсь Вашим опытом особенно п.6, вернее уже воспользовался.Но это технология а не этапность. Думаю Вам хорошо известно, что стройной теории нелинейных моделей нет и все решают частные задачи.
Тяжеловато делать какие-либо выводы из статьи.
Вторая функция ode из этого модуля реализует несколько методов, в том числе и упомянутый пятиранговый метод Рунге-Кутта-Фельберга, но, вследствие универсальности, имеет ограниченное быстродействие.
Вот что это значит — «вследствие универсальности»? Универсальность не может приводить к потере быстродействия, это просто разные характеристики.
Необходимо провести исследование быстродействия.
Исследование быстродействия, по моему скромному мнению, отвечает на вопрос «а почему быстродействие такое и как это исправить». Вы просто сравнили два конкретных примера без какого-либо объяснения, почему именно эти примеры и где же именно происходит просадка быстродействия. Т.е. это не исследование.
Целью настоящей публикации является сравнительный анализ
Сравнение — есть, но где же анализ?
Другими словами, мне вообще непонятно как (и можно ли) из статьи сделать PR в репо scipy.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Истории