Ищем интегратора/решателя ОДУ со спокойным отношением к точности производных

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

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

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

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

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

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

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

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

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

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

Не уверен, что внес большой вклад, но в мире фармацевтического моделирования мы используем LSODE, DVERK и DGPADM. DVERK — хороший быстрый простой решатель Рунге-Кутты порядка 5/6. DGPADM — хороший решатель матрицы-экспоненты. Если ваши ОДУ линейны, матричный экспонент на сегодняшний день лучше. Но ваша проблема немного в другом.

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

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

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

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

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

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

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