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

Иллюстрация численного интегрирования дифференциального уравнения $y'=y,y(0)=1.$

  Синий: метод Эйлера

  Зеленый: метод средней точки

  Красный: Точное решение: ${\textstyle y=e^{t}}$.

Размер шага $h=1.0$.

Та же иллюстрация для $h=0.25.$ Метод средней точки сходится быстрее, чем метод Эйлера, так как $h\to 0$.

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

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

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

Проблема

Дифференциальное уравнение первого порядка представляет собой задачу начального значения (IVP) вида:

$y'(t)=f(t,y(t)),\qquad y(t_{0})=y_{0},$ ( 1 )

где $f$ это функция $f:[t_{0},\infty )\times \mathbb {R} ^{d}\to \mathbb {R} ^{d}$, а начальное условие $y_{0}\in \mathbb {R} ^{d}$ является заданным вектором. Первый порядок означает, что в уравнении присутствует только первая производная от y , а более высокие производные отсутствуют.

Не ограничивая общности систем высшего порядка, мы ограничимся дифференциальными уравнениями первого порядка , поскольку ОДУ высшего порядка можно преобразовать в более крупную систему уравнений первого порядка введением дополнительных переменных. Например, уравнение второго порядка y " = − y можно переписать в виде двух уравнений первого порядка: y ′ = z и z ′ = − y .

В этом разделе мы описываем численные методы решения IVP и отмечаем, что краевые задачи (BVP) требуют другого набора инструментов. В BVP значения или компоненты решения y определяются более чем в одной точке. По этой причине для решения БВП необходимо использовать разные методы. Например, метод стрельбы (и его варианты) или глобальные методы вроде конечных разностей , методы Галеркина , или методы коллокации подходят для этого класса проблем.

Теорема Пикара –Линделёфа утверждает, что существует единственное решение, если f липшицево -непрерывно .

Методы

Численные методы решения ИВП первого порядка часто попадают в одну из двух больших категорий: линейные многошаговые методы , или методы Рунге–Кутты . Дальнейшее разделение можно осуществить, разделив методы на явные и неявные. Например, к неявным линейным многошаговым методам относятся методы Адамса-Моултона и методы обратного дифференцирования (BDF), тогда как неявные методы Рунге-Кутты включают диагонально неявный Рунге – Кутта (DIRK), однодиагонально неявный Рунге – Кутта (SDIRK), и Гаусс-Радау (на основе квадратуры Гаусса ) численные методы. Явные примеры из линейного многошагового семейства включают методы Адамса-Бэшфорта , и любой метод Рунге-Кутты с нижней диагональной таблицей Батчера является явным . Общее эмпирическое правило гласит, что жесткие дифференциальные уравнения требуют использования неявных схем, тогда как нежесткие задачи можно решать более эффективно с помощью явных схем.

Так называемые общие линейные методы (ОМЛМ) являются обобщением двух вышеуказанных больших классов методов.

метод Эйлера

Дополнительная информация: метод Эйлера.

Из любой точки кривой можно найти приближение к ближайшей точке кривой, переместившись на небольшое расстояние по касательной к кривой.

Начиная с дифференциального уравнения ( 1 ), заменим производную y ′ конечно-разностной аппроксимацией

$y'(t)\approx {\frac {y(t+h)-y(t)}{h}},$ ( 2 )

которая при перестановке дает следующую формулу

$y(t+h)\approx y(t)+hy'(t)$

и использование ( 1 ) дает:

$y(t+h)\approx y(t)+hf(t,y(t)).$ ( 3 )

Эта формула обычно применяется следующим образом. Выбираем размер шага h и строим последовательность $t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h,...$ Обозначим через $y_{n}$ численная оценка точного решения $y(t_{n})$. Руководствуясь ( 3 ), мы вычисляем эти оценки по следующей рекурсивной схеме

$y_{n+1}=y_{n}+hf(t_{n},y_{n}).$ ( 4 )

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

Метод Эйлера является примером явного метода. Это означает, что новое значение y _{n +1} определяется с точки зрения уже известных вещей, например y _n .

Обратный метод Эйлера

Дополнительная информация: Обратный метод Эйлера.

Если вместо ( 2 ) использовать приближение

$y'(t)\approx {\frac {y(t)-y(t-h)}{h}},$ ( 5 )

мы получаем обратный метод Эйлера :

$y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1}).$ ( 6 )

Обратный метод Эйлера — это неявный метод, означающий, что нам нужно решить уравнение, чтобы найти y _{n +1} . часто используют итерацию с фиксированной точкой или (некоторую модификацию) метода Ньютона – Рафсона Для этого .

Решение этого уравнения требует больше времени, чем явные методы; эту стоимость необходимо учитывать при выборе метода, который будет использоваться. Преимущество неявных методов, таких как ( 6 ), заключается в том, что они обычно более стабильны для решения жесткого уравнения больший размер шага h , а это означает, что можно использовать .

Метод экспоненциального интегратора первого порядка

Дополнительная информация: Экспоненциальный интегратор.

Экспоненциальные интеграторы описывают большой класс интеграторов, которые в последнее время получили большое развитие. Они датируются как минимум 1960-ми годами.

Вместо ( 1 ) мы предполагаем, что дифференциальное уравнение имеет вид

$y'(t)=-A\,y+{\mathcal {N}}(y),$ ( 7 )

или он был локально линеаризован относительно фонового состояния для получения линейного члена $-Ay$ и нелинейный член ${\mathcal {N}}(y)$.

Экспоненциальные интеграторы строятся путем умножения ( 7 ) на ${\textstyle e^{At}}$и точно интегрируя результат повременной интервал $[t_{n},t_{n+1}=t_{n}+h]$:

$y_{n+1}=e^{-Ah}y_{n}+\int _{0}^{h}e^{-(h-\tau )A}{\mathcal {N}}\left(y\left(t_{n}+\tau \right)\right)\,d\tau .$

Это интегральное уравнение является точным, но оно не определяет интеграл.

Экспоненциальный интегратор первого порядка можно реализовать, если ${\mathcal {N}}(y(t_{n}+\tau ))$ постоянная на всем интервале:

$y_{n+1}=e^{-Ah}y_{n}+A^{-1}(1-e^{-Ah}){\mathcal {N}}(y(t_{n}))\ .$ ( 8 )

Обобщения

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

Одна из возможностей состоит в том, чтобы использовать не только ранее вычисленное значение y _n для определения y _{n +1} , но и сделать решение зависящим от большего количества прошлых значений. Это дает так называемый многошаговый метод . Возможно, самым простым является метод чехарды , который является вторым порядком и (грубо говоря) опирается на два значения времени.

Почти все практические многошаговые методы относятся к семейству линейных многошаговых методов , которые имеют вид $${\begin{aligned}&{}\alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}\\&{}\quad =h\left[\beta _{k}f(t_{n+k},y_{n+k})+\beta _{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +\beta _{0}f(t_{n},y_{n})\right].\end{aligned}}$$

Другая возможность — использовать больше точек в интервале. $[t_{n},t_{n+1}]$. Это приводит к семейству методов Рунге-Кутты , названному в честь Карла Рунге и Мартина Кутты . Особой популярностью пользуется один из их методов четвертого порядка.

Расширенные функции

Хорошая реализация одного из этих методов решения ОДУ влечет за собой нечто большее, чем просто формулу шага по времени.

Зачастую неэффективно постоянно использовать один и тот же размер шага, поэтому методы с переменным размером шага были разработаны . Обычно размер шага выбирается таким, чтобы (локальная) ошибка на шаг была ниже некоторого уровня допуска. Это означает, что методы также должны вычислять индикатор ошибки — оценку локальной ошибки.

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

Другие желательные функции включают в себя:

• плотный вывод : дешевые численные аппроксимации для всего интервала интегрирования, а не только в точках t ₀ , t ₁ , t ₂ ,...
• Местоположение события : определение моментов, когда, скажем, конкретная функция исчезает. Обычно для этого требуется использование алгоритма поиска корня .
• поддержка параллельных вычислений .
• при использовании для интегрирования по времени обратимость времени

Альтернативные методы

Многие методы не подпадают под обсуждаемые здесь рамки. Некоторые классы альтернативных методов:

• мультипроизводные методы , в которых используется не только функция f , но и ее производные. В этот класс входят методы Эрмита-Обрешкова и методы Фельберга , а также такие методы, как метод Паркера-Сохацкого. или метод Бычкова–Щербакова, который рекурсивно вычисляет коэффициенты ряда Тейлора решения y .
• методы для ОДУ второго порядка . Мы сказали, что все ОДУ более высокого порядка можно преобразовать к ОДУ первого порядка вида (1). Хотя это, безусловно, правда, возможно, это не лучший способ действовать. В частности, методы Нистрема работают непосредственно с уравнениями второго порядка.
• методы геометрического интегрирования специально разработаны для специальных классов ОДУ (например, симплектических интеграторов для решения гамильтоновых уравнений ). Они заботятся о том, чтобы численное решение учитывало базовую структуру или геометрию этих классов.
• Методы квантовых систем состояний представляют собой семейство методов интеграции ОДУ, основанных на идее квантования состояний. Они эффективны при моделировании разреженных систем с частыми разрывами.

Параллельные во времени методы

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

Ранние методы PinT (самые ранние из них были предложены в 1960-х годах) изначально были упущены из виду исследователями из-за того, что требуемые им архитектуры параллельных вычислений еще не были широко доступны. Благодаря увеличению доступной вычислительной мощности интерес возобновился в начале 2000-х годов с разработкой Parareal , гибкого и простого в использовании алгоритма PinT, который подходит для решения широкого спектра задач IVP. Появление экзафлопсных вычислений привело к тому, что алгоритмы PinT привлекают все большее внимание исследователей и разрабатываются таким образом, чтобы они могли использовать самые мощные в мире суперкомпьютеры . К наиболее популярным методам по состоянию на 2023 год относятся Parareal, PFASST, ParaDiag и MGRIT.

Анализ

Численный анализ – это не только разработка численных методов, но и их анализ. Тремя центральными концепциями в этом анализе являются:

• сходимость : приближает ли метод решение,
• порядок : насколько хорошо оно приближает решение, и
• стабильность : устраняются ли ошибки.

Конвергенция

Основные статьи: Последовательность , Предел (математика) и Предел последовательности.

Численный метод называется сходящимся, если численное решение приближается к точному решению при стремлении размера шага h к 0. Точнее, мы требуем, чтобы для любого ОДУ (1) с липшицевой функцией f и любого t ^* > 0,

$\lim _{h\to 0^{+}}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\left\|y_{n,h}-y(t_{n})\right\|=0.$

Все упомянутые выше методы являются конвергентными.

Последовательность и порядок

Дополнительная информация: Ошибка усечения (численное интегрирование).

Предположим, что численный метод

$y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,$

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

$\delta _{n+k}^{h}=\Psi \left(t_{n+k};y(t_{n}),y(t_{n+1}),\dots ,y(t_{n+k-1});h\right)-y(t_{n+k}).$

Метод называется состоятельным, если

$\lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.$

Метод имеет порядок $p$ если

$\delta _{n+k}^{h}=O(h^{p+1})\quad {\mbox{as }}h\to 0.$

Следовательно, метод является непротиворечивым, если он имеет порядок больше 0. (Прямой) метод Эйлера (4) и обратный метод Эйлера (6), представленные выше, имеют порядок 1, поэтому они непротиворечивы. Большинство методов, используемых на практике, имеют более высокий порядок. Согласованность – необходимое условие конвергенции , но недостаточно; Чтобы метод был сходящимся, он должен быть непротиворечивым и нулевым .

Связанное с этим понятие — глобальная ошибка (ошибка усечения) , ошибка, возникающая на всех этапах, необходимых для достижения фиксированного времени. $t$. Явно глобальная ошибка во времени $t$ является $y_{N}-y(t)$ где $N=(t-t_{0})/h$. Глобальная ошибка $p$Одношаговый метод порядка $O(h^{p})$; в частности, такой метод является сходящимся. Это утверждение не обязательно верно для многошаговых методов.

Стабильность и жесткость

Дополнительная информация: Жесткое уравнение.

Для некоторых дифференциальных уравнений применение стандартных методов, таких как метод Эйлера, явные методы Рунге-Кутты или многошаговые методы (например, методы Адамса-Бэшфорта), демонстрирует нестабильность решений, хотя другие методы могут давать устойчивые решения. Это «трудное поведение» в уравнении (которое само по себе не обязательно может быть сложным) описывается как жесткость и часто вызвано наличием различных временных масштабов в основной проблеме. Например, столкновение в механической системе, такой как ударный генератор, обычно происходит в гораздо меньшем временном масштабе, чем время движения объектов; это несоответствие приводит к очень «крутым поворотам» кривых параметров состояния.

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

История

Ниже представлен график некоторых важных событий в этой области.

• 1768 г. – Леонард Эйлер публикует свой метод.
• 1824 г. — Огюстен Луи Коши доказывает сходимость метода Эйлера. В этом доказательстве Коши использует неявный метод Эйлера.
• 1855 — Первое упоминание о многошаговых методах Джона Кауча Адамса в письме Фрэнсиса Башфорта .
• 1895 — Карл Рунге публикует первый метод Рунге-Кутты .
• 1901 — Мартин Кутта четвертого порядка описывает популярный метод Рунге–Кутты .
• 1910 — Льюис Фрай Ричардсон объявляет о своем методе экстраполяции — экстраполяции Ричардсона .
• 1952 – Чарльз Ф. Кертисс и Джозеф Окленд Хиршфельдер вводят термин «жесткие уравнения» .
• 1963 - Гермунд Далквист представляет A-стабильность методов интегрирования.

Численные решения одномерных краевых задач второго порядка

Краевые задачи (БВП) обычно решаются численно путем решения приблизительно эквивалентной матричной задачи, полученной путем дискретизации исходной БВП. Наиболее часто используемый метод численного решения BVP в одном измерении называется методом конечных разностей . Этот метод использует преимущества линейных комбинаций точечных значений для построения коэффициентов конечной разности , которые описывают производные функции. Например, аппроксимация центральной разностью второго порядка для первой производной определяется следующим образом:

${\frac {u_{i+1}-u_{i-1}}{2h}}=u'(x_{i})+{\mathcal {O}}(h^{2}),$

второго порядка а центральная разность для второй производной определяется выражением:

${\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}=u''(x_{i})+{\mathcal {O}}(h^{2}).$

В обеих этих формулах $h=x_{i}-x_{i-1}$ — расстояние между соседними значениями x в дискретизированной области. Затем строится линейная система, которую затем можно решить стандартными матричными методами . Например, предположим, что уравнение, которое необходимо решить, имеет следующий вид:

${\begin{aligned}&{}{\frac {d^{2}u}{dx^{2}}}-u=0,\\&{}u(0)=0,\\&{}u(1)=1.\end{aligned}}$

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

$u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}$

и решить полученную систему линейных уравнений. Это приведет к таким уравнениям, как:

${\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}-u_{i}=0,\quad \forall i={1,2,3,...,n-1}.$

На первый взгляд кажется, что эта система уравнений имеет трудности, связанные с тем, что в уравнении нет членов, которые не умножаются на переменные, но на самом деле это неверно. При i = 1 и n − 1 существует член, включающий граничные значения $u(0)=u_{0}$ и $u(1)=u_{n}$ а так как эти две величины известны, то можно просто подставить их в это уравнение и в результате получить неоднородную систему линейных уравнений , имеющую нетривиальные решения.

См. также

• Условие Куранта – Фридрихса – Леви.
• Энергетический дрейф
• Общие линейные методы
• Список тем численного анализа#Численные методы для обыкновенных дифференциальных уравнений
• Обратимый алгоритм распространения системы отсчета
• Язык Modelica и OpenModelica программное обеспечение

Ссылки

• Брэди, Брайан (2006). Дружественное введение в численный анализ . Река Аппер-Сэддл, Нью-Джерси: Пирсон Прентис Холл. ISBN  978-0-13-013054-9 .
• Дж. К. Батчер , Численные методы решения обыкновенных дифференциальных уравнений , ISBN   0-471-96758-0
• Эрнст Хайрер, Сиверт Пауль Норсетт и Герхард Ваннер, Решение обыкновенных дифференциальных уравнений I: нежесткие задачи, второе издание, Springer Verlag, Берлин, 1993. ISBN   3-540-56670-8 .
• Эрнст Хайрер и Герхард Ваннер, Решение обыкновенных дифференциальных уравнений II: жесткие и дифференциально-алгебраические задачи, второе издание, Springer Verlag, Берлин, 1996. ISBN   3-540-60452-9 .
  (Эта двухтомная монография систематически охватывает все аспекты этой области.)
• Хохбрук, Марлис ; Остерманн, Александр (май 2010 г.). «Экспоненциальные интеграторы». Акта Нумерика . 19 : 209–286. Бибкод : 2010AcNum..19..209H . CiteSeerX   10.1.1.187.6794 . дои : 10.1017/S0962492910000048 . S2CID   4841957 .
• Арье Изерлес, Первый курс численного анализа дифференциальных уравнений, Cambridge University Press, 1996. ISBN   0-521-55376-8 (в твердом переплете), ISBN   0-521-55655-4 (мягкая обложка).
  (Учебник, предназначенный для студентов и аспирантов, изучающих математику, в котором также обсуждаются численные уравнения в частных производных .)
• Джон Денхольм Ламберт, Численные методы для обыкновенных дифференциальных систем, John Wiley & Sons, Чичестер, 1991. ISBN   0-471-92990-5 .
  (Учебник, несколько более требовательный, чем книга Изерлеса.)

Внешние ссылки

• Джозеф В. Рудмин, Применение метода Паркера-Сочацкого к небесной механике. Архивировано 16 мая 2016 г. в Португальском веб-архиве , 1998 г.
• Доминик Турне, Приблизительное интегрирование обыкновенных дифференциальных уравнений (1671–1914) , докторская диссертация Парижского университета 7 - Дени Дидро, июнь 1996 г. Réimp. Вильнев д'Аск: Press universitaire du Septentrion, 1997, 468 стр. (Обширный онлайн-материал по истории численного анализа ОДУ; англоязычные материалы по истории численного анализа ОДУ см., например, в цитируемых им бумажных книгах Чаберта и Голдстайна.)
• Пчелинцев, АН (2020). «Точный численный метод и алгоритм построения решений хаотических систем». Журнал прикладной нелинейной динамики . 9 (2): 207–221. arXiv : 2011.10664 . дои : 10.5890/JAND.2020.06.004 . S2CID   225853788 .
• kv на GitHub ( библиотека C++ со строгими решателями ОДУ)
• INTLAB (библиотека, созданная MATLAB / GNU Octave , которая включает в себя строгие решатели ОДУ)