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

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

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

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

Дифференциальное уравнение первого порядка — это задача начального значения (IVP) вида ^[2]

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

( 1 )

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

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

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

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

Численные методы решения задач первого порядка часто попадают в одну из двух больших категорий: ^[5] линейные многошаговые методы или методы Рунге–Кутты . Дальнейшее разделение может быть реализовано путем разделения методов на явные и неявные. Например, неявные линейные многошаговые методы включают методы Адамса–Моултона и методы обратного дифференцирования (BDF), тогда как неявные методы Рунге–Кутты ^[6] включают диагонально неявные методы Рунге–Кутты (DIRK), ^{[7] [8]} однократно диагонально неявные методы Рунге–Кутты (SDIRK), ^[9] и Гаусса–Радау ^[10] (основанные на квадратурных уравнениях Гаусса ^[11] ) численные методы. Явные примеры из семейства линейных многошаговых методов включают методы Адамса–Башфорта , и любой метод Рунге–Кутты с нижней диагональной таблицей Бутчера является явным . Нестрогое эмпирическое правило гласит, что жесткие дифференциальные уравнения требуют использования неявных схем, тогда как нежесткие задачи можно решить более эффективно с помощью явных схем.

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

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

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

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

( 2 )

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

${\displaystyle y(t+h)\approx y(t)+hy'(t)}$

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

${\displaystyle y(t+h)\approx y(t)+hf(t,y(t)).}$

( 3 )

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

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

( 4 )

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

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

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

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

( 5 )

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

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

( 6 )

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

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

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

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

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

( 7 )

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

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

${\displaystyle 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 .}$

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

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

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

( 8 )

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

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

Почти все практические многошаговые методы попадают в семейство линейных многошаговых методов , которые имеют вид$${\displaystyle {\begin{align}&{}\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{align}}}$$

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

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

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

Расширение этой идеи заключается в динамическом выборе между различными методами разных порядков (это называется методом переменного порядка ). Методы, основанные на экстраполяции Ричардсона ^[14] , такие как алгоритм Булирша–Штера ^{[15] [16],} часто используются для построения различных методов разных порядков.

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

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

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

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

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

Ранние методы PinT (самые ранние были предложены в 1960-х годах) ^[20] изначально не были замечены исследователями из-за того, что требуемые им параллельные вычислительные архитектуры еще не были широко доступны. С появлением большей вычислительной мощности интерес возобновился в начале 2000-х годов с разработкой Parareal , гибкого, простого в использовании алгоритма PinT, который подходит для решения широкого спектра IVP. Появление exascale вычислений означало, что алгоритмы PinT привлекают все большее внимание исследователей и разрабатываются таким образом, чтобы они могли использовать самые мощные суперкомпьютеры в мире . Наиболее популярные методы по состоянию на 2023 год включают Parareal, PFASST, ParaDiag и MGRIT. ^[21]

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

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

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

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

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

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

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

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

${\displaystyle \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}).}$

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

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

Метод имеет порядок , если${\displaystyle p}$

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

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

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

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

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

Ниже приведена хронология некоторых важных событий в этой области. ^{[26] [27]}

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

Краевые задачи (BVP) обычно решаются численно путем решения приблизительно эквивалентной матричной задачи, полученной путем дискретизации исходной BVP. ^[28] Наиболее часто используемый метод численного решения BVP в одном измерении называется методом конечных разностей . ^[3] Этот метод использует линейные комбинации точечных значений для построения коэффициентов конечных разностей , которые описывают производные функции. Например, центральное разностное приближение второго порядка к первой производной определяется как:

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

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

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

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

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

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

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

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

${\displaystyle {\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 есть член, включающий граничные значения и и поскольку эти два значения известны, можно просто подставить их в это уравнение и в результате получить неоднородную систему линейных уравнений , которая имеет нетривиальные решения.${\displaystyle u(0)=u_{0}}$${\displaystyle u(1)=u_{n}}$

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

• Брэди, Брайан (2006). Дружественное введение в численный анализ . Верхняя Сэддл-Ривер, Нью-Джерси: Pearson Prentice Hall. 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 . (Эта двухтомная монография систематически охватывает все аспекты этой области.) 
  
• Hochbruck, Marlis ; Ostermann, Alexander (май 2010 г.). "Экспоненциальные интеграторы". Acta Numerica . 19 : 209–286. Bibcode :2010AcNum..19..209H. CiteSeerX  10.1.1.187.6794 . doi :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 г.
• Доминик Турнес, «L'intégration approchée des équations différentielles ordinaires» (1671–1914) , докторская степень Парижского университета 7 – Дени Дидро, июнь 1996 г. Réimp. Вильнёв д'Аск: Press universitaires du Septentrion, 1997, 468 стр. (Обширный онлайн-материал по истории численного анализа ОДУ; англоязычные материалы по истории численного анализа ОДУ см., например, в цитируемых им бумажных книгах Чаберта и Голдстайна.)
• Пчелинцев, АН (2020). «Точный численный метод и алгоритм построения решений хаотических систем». Журнал прикладной нелинейной динамики . 9 (2): 207–221. arXiv : 2011.10664 . doi :10.5890/JAND.2020.06.004. S2CID  225853788.
• kv на GitHub ( библиотека C++ со строгими решателями ОДУ)
• INTLAB (библиотека, созданная MATLAB / GNU Octave , включающая строгие решатели ОДУ)