Ищите интегратора / решателя ODE с непринужденным отношением к производной точности

Question

У меня есть система (первого порядка) ODE с довольно дорогими для вычисления производных.

Однако, производные могут быть вычислены значительно дешевле в пределах заданных границ ошибки, либо потому, что производные вычисляются из сходящегося ряда и границы могут быть наложены на максимальный вклад от пропущенных терминов, либо путем использования предварительно вычисленной информации о диапазоне, хранящейся в таблицы поиска kd-tree / octree.

К сожалению, я не смог найти каких-либо общих решателей ODE, которые могли бы извлечь из этого пользу; кажется, все они просто дают вам координаты и хотят получить точный результат. (Имейте в виду, я не эксперт по ODE; я знаком с Рунге-Куттой, материалом из книги «Числовые рецепты», LSODE и решателя научной библиотеки Gnu).

то есть для всех решателей, которые я видел, вы предоставляете функцию обратного вызова derivs, принимающую t и массив x и возвращающую массив dx/dt назад; но в идеале я ищу тот, который дает обратный вызов t, x s, и массив допустимых ошибок и получает обратно массивы dx/dt_min и dx/dt_max с производным диапазоном гарантированно в пределах требуемой точности. (Возможно, существует множество не менее полезных вариантов).

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

Answer 1

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

Answer 2

Подумав об этом еще немного, мне пришло в голову, что интервальная арифметика, вероятно, является ключевой. Моя derivs функция в основном возвращает интервалы. Интегратор, использующий арифметику интервалов, будет поддерживать x как интервалы. Все, что меня интересует, - это получение достаточно маленькой ошибки, связанной с x s в конечном t. Очевидный подход заключается в итеративной реинтеграции, улучшающей качество выборки, вводящей наибольшее количество ошибок на каждой итерации, пока мы, наконец, не получим результат с приемлемыми границами (хотя это звучит так, как если бы это было «лекарство хуже, чем болезнь» в отношении для общей эффективности). Я подозреваю, что адаптивное управление размером шага могло бы хорошо вписаться в такую ​​схему, при этом размер шага был выбран таким образом, чтобы «неявная» ошибка дискретизации была сопоставима с «явной ошибкой», т. Е. Интервалом).

В любом случае, прибегая к помощи "арифметики интервала решения ode" или просто "интервальной оды", вы обнаруживаете массу новых интересных и актуальных вещей ( VNODE и, в частности, ссылки).

Answer 3

Я не уверен, что это правильный вопрос.

Во многих алгоритмах, например, решение нелинейных уравнений, f (x) = 0, оценка производной f '(x) - это все, что требуется для использования в чем-то подобном методу Ньютона, поскольку вы только нужно идти в «общем направлении» ответа.

Однако в этом случае производная является основной частью уравнения (ODE), которое вы решаете - получите производную неправильно, и вы просто получите неправильный ответ; это все равно что пытаться решить f (x) = 0 только с приближением для f (x).

Как показал другой ответ, если вы настроите свой ODE как примененный f (x) + g (x), где g (x) - это термин ошибки, вы должны быть в состоянии связать ошибки в ваших производных с ошибками в вашем входы.

Answer 4

Грубо говоря, если вы знаете f 'до абсолютной ошибки eps и интегрируете от x0 до x1, погрешность интеграла, возникающая из ошибки в производной, будет <= eps * (x1 - x0). Существует также ошибка дискретизации, исходящая от вашего ODE-решателя. Подумайте, насколько большим может быть eps * (x1 - x0) для вас, и задайте для солдера ODE значения f ', вычисленные с ошибкой <= eps. </p>

Answer 5

Проверка на метод конечных элементов с линейными базисными функциями и квадратурой средней точки. Решение следующего ODE требует только одного вычисления каждого из f (x), k (x) и b (x) на элемент:

-k (x) u '' (x) + b (x) u '(x) = f (x)

Ответ будет иметь точечную ошибку, пропорциональную ошибке в ваших оценках.

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

Answer 6

Не уверен, что я много помогаю, но в мире фарм-моделирования мы используем LSODE, DVERK и DGPADM. DVERK - это хороший быстрый простой 5/6 Runge-Kutta решатель. DGPADM - хороший решатель показателей матрицы. Если ваши ODE линейны, матричный показатель является наилучшим. Но ваша проблема немного другая.

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

Возможно, вы выходите на новую теоретическую территорию. Удачи!

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

Answer 7

Вы изучали использование odeset ? Он позволяет вам устанавливать параметры для солвера ODE, а затем вы передаете структуру опций в качестве четвертого аргумента для любого солвера, который вы вызываете. Свойства контроля ошибок (RelTol, AbsTol, NormControl) могут представлять для вас наибольший интерес. Не уверен, что это именно та помощь, которая вам нужна, но это лучшее предложение, которое я мог бы предложить, когда последний раз использовал функции MATLAB ODE много лет назад.

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