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

Иллюстрация численного интегрирования для дифференциального уравнения синего света: метод Эйлера , Зеленый: метод средней точки , красный: точное решения, . Размер шага .${\ displaystyle y '= y, y (0) = 1.}$${\textstyle y=e^{t}}$${\displaystyle h=1.0}$

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

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

Многие дифференциальные уравнения не могут быть решены с помощью символьных вычислений («анализа»). Однако для практических целей, например, в инженерии, часто бывает достаточно численного приближения к решению. В алгоритмах изученных здесь могут быть использованы для вычисления такого приближения. Альтернативный метод - использовать методы исчисления для получения разложения решения в ряд .

Обычные дифференциальные уравнения встречаются во многих научных дисциплинах, включая физику , химию , биологию и экономику . ^[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] или методы коллокации , подходят для этого класса проблем.

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

Методы [ править ]

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

Так называемые общие линейные методы (GLM) являются обобщением двух вышеупомянутых больших классов методов. ^[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 и строим последовательность t ₀ , t ₁ = t ₀ + h , t ₂ = t ₀ + 2 h ,… Обозначим через y _n численную оценку точного решения y ( t _n ) . Руководствуясь ( 3 ), мы вычисляем эти оценки по следующей рекурсивной схеме

${\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(t-h)}{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 -Ay}$${\displaystyle {\mathcal {N}}(y)}$

Экспоненциальные интеграторы строятся путем умножения ( 7 ) на и точного интегрирования результата за интервал времени :${\textstyle e^{At}}$${\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} , но и сделать решение зависимым от большего количества прошлых значений. Это дает так называемый многоступенчатый метод . Возможно, самым простым является метод чехарды второго порядка, который (грубо говоря) основан на двух значениях времени.

Практически все практические многоступенчатые методы относятся к семейству линейных многоступенчатых методов , которые имеют вид

${\displaystyle \alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}}$

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

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

Расширенные функции [ править ]

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

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

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

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

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

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

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

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

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

Для приложений, требующих параллельных вычислений на суперкомпьютерах , степень параллелизма, предлагаемая численным методом, становится актуальной. Ввиду проблем, связанных с экзадачными вычислительными системами, изучаются численные методы для задач начального значения, которые могут обеспечить параллелизм во временном направлении. ^[20] Parareal - относительно хорошо известный пример такого метода интеграции параллельно во времени , но его ранние идеи восходят к 1960-м годам. ^[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, поэтому они согласованы. Большинство применяемых на практике методов достигают более высокого порядка. Согласованность - необходимое условие конвергенции ^{[ необходима цитата ]} , но ее недостаточно; Чтобы метод был сходимым, он должен быть как согласованным, так и устойчивым к нулю .

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

Устойчивость и жесткость [ править ]

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

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

История [ править ]

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

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

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

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

Внешние ссылки [ править ]

• Джозеф В. Рудмин, Применение метода Паркера – Сохацкого в небесной механике , 1998.
• Доминик Турнес, L'intégration Approchée des équations différentielles ordinaires (1671-1914) , докторская диссертация Парижского университета 7 - Дени Дидро, июнь 1996. Réimp. Вильнёв д'Аск: Presses Universitaires du Septentrion, 1997, 468 с. (Обширный онлайн-материал по истории численного анализа ODE, англоязычный материал по истории численного анализа ODE, см., Например, цитируемые им бумажные книги Чабера и Голдстайна.)
• Пчелинцев, АН (2020). «Точный численный метод и алгоритм построения решений хаотических систем» (PDF) . Журнал прикладной нелинейной динамики . 9 (2): 207–221. DOI : 10.5890 / JAND.2020.06.004 .
• kv на GitHub ( библиотека C ++ со строгими решателями ODE)
• INTLAB (библиотека, созданная MATLAB / GNU Octave, которая включает строгие решатели ODE)