Численные методы для обыкновенных дифференциальных уравнений - Numerical methods for ordinary differential equations

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

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

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

Эта проблема

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

${\ Displaystyle у '(т) = е (т, у (т)), \ qquad у (т_ {0}) = у_ {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 необходимо использовать разные методы. Например, метод стрельбы (и его варианты) или глобальные методы, такие как конечные разности , методы Галеркина или методы коллокации , подходят для этого класса проблем.

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

Методы

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

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

Метод Эйлера

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

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

${\ Displaystyle у '(т) \ приблизительно {\ гидроразрыва {у (т + ч) -у (т)} {ч}},}$

( 2 )

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

${\ Displaystyle у (т + ч) \ приблизительно у (т) + хи '(т)}$

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

${\ Displaystyle у (т + ч) \ приблизительно у (т) + ВЧ (т, у (т)).}$

( 3 )

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

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

( 4 )

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

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

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

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

${\ Displaystyle у '(т) \ приблизительно {\ гидроразрыва {у (т) -у (й)} {ч}},}$

( 5 )

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

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

( 6 )

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

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

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

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

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

${\ Displaystyle у '(т) = - А \, у + {\ mathcal {N}} (у),}$

( 7 )

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

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

${\ 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}} (у (т_ {п} + \ тау))}$

${\ 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 {\ 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 [т_ {п}, т_ {п + 1}]}$

Расширенные возможности

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

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

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

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

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

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

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

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

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

Для приложений, требующих параллельных вычислений на суперкомпьютерах , степень параллелизма, предлагаемая численным методом, становится актуальной. Ввиду проблем, связанных с экзадачными вычислительными системами, изучаются численные методы для задач начального значения, которые могут обеспечить параллелизм во временном направлении. Parareal - относительно хорошо известный пример такого метода интеграции во времени , но его ранние идеи восходят к 1960-м годам.

Анализ

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

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

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

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

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

Устойчивость и жесткость

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

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

История

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

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

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

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

${\ 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 {align} & {} {\ frac {d ^ {2} u} {dx ^ {2}}} - u = 0, \\ & {} u (0) = 0, \\ & {} u (1) = 1. \ end {выровнено}}}$

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

${\ 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 и (1) = и_ {п}}$

Смотрите также

• Условие Куранта – Фридрихса – Леви.
• Энергетический дрейф
• Общие линейные методы
• Список тем численного анализа # Численные методы для обыкновенных дифференциальных уравнений
• Алгоритм обратимого распространения системы отсчета
• Язык 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)