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

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

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

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

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

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

Проблема

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

${ Displaystyle у '(т) = е (т, у (т)), qquad у (т_ {0}) = у_ {0}, qquad qquad (1)}$

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

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

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

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

Методы

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

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

Метод Эйлера

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

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

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

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

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

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

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

Эта формула обычно применяется следующим образом. Выбираем размер шага час, и построим последовательность т₀, т₁ = т₀ + час, т₂ = т₀ + 2час,… Обозначим через у_п численная оценка точного решения у(т_п). Руководствуясь (3), мы вычисляем эти оценки с помощью следующих рекурсивный схема

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Экспоненциальные интеграторы строятся умножением (7) на ${ displaystyle 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 })) . qquad qquad (8)}$

Обобщения

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

Одна из возможностей - использовать не только ранее вычисленное значение. у_п определить у_п+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].}$

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

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

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

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

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

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

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

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

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

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

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

Для приложений, требующих параллельные вычисления на суперкомпьютеры степень параллелизма, предлагаемая численным методом, становится актуальной. Учитывая вызовы со стороны эксафлопсные вычисления системы, численные методы для проблемы начального значения которые могут обеспечить параллелизм во временном направлении, изучаются.^[20]Парареальный является относительно известным примером такого параллельно во времени метод интеграции, но ранние идеи восходят к 1960-м годам.^[21]

Анализ

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

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

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

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

${ displaystyle lim _ {h to 0 +} max _ {n = 0,1, dots, lfloor t ^ {*} / h rfloor} | y_ {n, h} -y (t_ {n}) | = 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, поэтому они согласованы. Большинство используемых на практике методов достигают более высокого порядка. Согласованность - необходимое условие сходимости^{[нужна цитата ]}, но не достаточно; чтобы метод был конвергентным, он должен быть как согласованным, так и нулевой стабильный.

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

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

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

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

История

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

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

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

Краевые задачи (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}).}$

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

${ displaystyle { frac {d ^ {2} u} {dx ^ {2}}} - u = 0,}$

${ displaystyle u (0) = 0,}$

${ displaystyle u (1) = 1.}$

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

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

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

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

• Условие Куранта – Фридрихса – Леви.
• Энергетический дрейф
• Общие линейные методы
• Список тем численного анализа # Численные методы для обыкновенных дифференциальных уравнений
• Алгоритм обратимого распространения системы отсчета
• 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 г.). «Экспоненциальные интеграторы». Acta Numerica. 19: 209–286. Bibcode:2010AcNum..19..209H. CiteSeerX  10.1.1.187.6794. Дои:10.1017 / S0962492910000048.
• Арие Изерлес, Первый курс численного анализа дифференциальных уравнений, Издательство Кембриджского университета, 1996. ISBN  0-521-55376-8 (переплет), ISBN  0-521-55655-4 (мягкая обложка).
  (Учебник, предназначенный для студентов и аспирантов по математике, в котором также обсуждается числовые уравнения в частных производных.)
• Джон Денхольм Ламберт, Численные методы для обыкновенных дифференциальных систем, Джон Уайли и сыновья, Чичестер, 1991. ISBN  0-471-92990-5.
  (Учебник, чуть более требовательный, чем книга Изерлеса.)

внешняя ссылка

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