Жесткое уравнение

Дифференциальное уравнение, демонстрирующее необычную неустойчивость

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

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

Мотивирующий пример

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

Рассмотрим задачу начального значения

Точное решение (показано голубым цветом) —

Мы ищем численное решение , демонстрирующее такое же поведение.

Рисунок (справа) иллюстрирует численные результаты для различных числовых интеграторов, применяемых к уравнению.

1. Метод Эйлера с размером шага резко колеблется и быстро выходит за пределы графика (показан красным). ${\displaystyle h={\tfrac {1}{4}}}$
2. Метод Эйлера с шагом в два раза меньше, , дает решение в пределах графика, но колеблется около нуля (показано зеленым цветом). ${\displaystyle h={\tfrac {1}{8}}}$
3. Метод трапеций (то есть двухэтапный метод Адамса–Моултона ) задается формулой где . Применение этого метода вместо метода Эйлера дает гораздо лучший результат (синий). Численные результаты монотонно убывают до нуля, как и точное решение. ${\displaystyle y'=f(t,y)}$

Одним из наиболее ярких примеров жестких обыкновенных дифференциальных уравнений (ОДУ) является система, описывающая химическую реакцию Робертсона: ^[1]

Если рассматривать эту систему на коротком интервале, например, то не возникает никаких проблем с численным интегрированием. Однако если интервал очень большой ( скажем, ¹⁰¹¹ ), то многие стандартные коды не могут правильно его интегрировать. ${\displaystyle t\in [0,40]}$

Коэффициент жесткости

Рассмотрим линейную неоднородную систему с постоянным коэффициентом

где и — постоянная, диагонализируемая матрица с собственными значениями (предполагаемыми различными) и соответствующими собственными векторами . Общее решение ( 5 ) принимает вид ${\displaystyle \mathbf {y} ,\mathbf {f} \in \mathbb {R} ^{n}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle n\times n}$ ${\displaystyle \lambda _{t}\in \mathbb {C} ,t=1,2,\ldots ,n}$ ${\displaystyle \mathbf {c} _{t}\in \mathbb {C} ^{n},t=1,2,\ldots ,n}$

где - произвольные константы, а - частный интеграл. Теперь предположим, что ${\displaystyle \kappa _{t}}$ ${\displaystyle \mathbf {g} (x)}$

что подразумевает, что каждый из членов равен , так что решение приближается асимптотически к ; член будет монотонно убывать, если является действительным, и синусоидально, если является комплексным. ${\displaystyle e^{\lambda _{t}x}\mathbf {c} _{t}\to 0}$ ${\displaystyle x\to \infty }$ ${\displaystyle \mathbf {y} (x)}$ ${\displaystyle \mathbf {g} (x)}$ ${\displaystyle x\to \infty }$ ${\displaystyle e^{\lambda _{t}x}\mathbf {c} _{t}}$ ${\displaystyle \lambda _{t}}$ ${\displaystyle \lambda _{t}}$

Интерпретируя как время (как это часто бывает в физических задачах), называется переходным решением и стационарным решением . Если велико, то соответствующий член будет быстро затухать по мере увеличения и, таким образом, называется быстрым переходным процессом ; если мало, то соответствующий член затухает медленно и называется медленным переходным процессом . Пусть определяется как ${\displaystyle x}$ ${\textstyle \sum _{t=1}^{n}\kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}}$ ${\displaystyle \mathbf {g} (x)}$ ${\displaystyle \left|\operatorname {Re} (\lambda _{t})\right|}$ ${\displaystyle \kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}}$ ${\displaystyle x}$ ${\displaystyle \left|\operatorname {Re} (\lambda _{t})\right|}$ ${\displaystyle \kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}}$ ${\displaystyle {\overline {\lambda }},{\underline {\lambda }}\in \{\lambda _{t},t=1,2,\ldots ,n\}}$

так что это самый быстрый переходный процесс и самый медленный. Теперь мы определяем коэффициент жесткости как ^[2] ${\displaystyle \kappa _{t}e^{{\overline {\lambda }}x}\mathbf {c} _{t}}$ ${\displaystyle \kappa _{t}e^{{\underline {\lambda }}x}\mathbf {c} _{t}}$

Характеристика жесткости

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

Дж. Д. Ламберт определяет жесткость следующим образом:

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

Существуют и другие характеристики, которые демонстрируются многими примерами жестких проблем, но для каждого из них есть контрпримеры, поэтому эти характеристики не являются хорошими определениями жесткости. Тем не менее, определения, основанные на этих характеристиках, широко используются некоторыми авторами и являются хорошими подсказками относительно наличия жесткости. Ламберт называет их «утверждениями», а не определениями, по вышеупомянутым причинам. Вот некоторые из них:

1. Линейная система с постоянными коэффициентами является жесткой, если все ее собственные значения имеют отрицательную действительную часть, а коэффициент жесткости велик.
2. Жесткость возникает, когда требования устойчивости, а не точности, ограничивают длину шага.
3. Жесткость возникает, когда некоторые компоненты раствора распадаются гораздо быстрее, чем другие. ^[3]

Этимология

Происхождение термина «жесткость» точно не установлено. Согласно Джозефу Окленду Хиршфельдеру , термин «жесткий» используется потому, что такие системы соответствуют жесткой связи между водителем и ведомыми сервомеханизмами . [ 4 ^] Согласно Ричарду Л. Бердену и Дж. Дугласу Фейресу,

Значительные трудности могут возникнуть при применении стандартных численных методов для приближенного решения дифференциального уравнения , когда точное решение содержит члены вида , где — комплексное число с отрицательной действительной частью. ${\displaystyle e^{\lambda t}}$ ${\displaystyle \lambda }$

. . .

Проблемы, связанные с быстро затухающими переходными решениями, естественным образом возникают в самых разных приложениях, включая изучение пружинных и демпфирующих систем, анализ систем управления и проблемы химической кинетики . Все это примеры класса задач, называемых жесткими (математическая жесткость) системами дифференциальных уравнений, из-за их применения при анализе движения пружинных и массовых систем, имеющих большие константы пружины (физическая жесткость ). ^[5]

Например, задача начального значения

с , , , можно записать в виде ( 5 ) с и ${\displaystyle m=1}$ ${\displaystyle c=1001}$ ${\displaystyle k=1000}$ ${\displaystyle n=2}$

и имеет собственные значения . Оба собственных значения имеют отрицательную действительную часть, а коэффициент жесткости равен ${\displaystyle {\overline {\lambda }}=-1000,{\underline {\lambda }}=-1}$

что довольно велико. Система ( 10 ) тогда, безусловно, удовлетворяет утверждениям 1 и 3. Здесь константа пружины велика, а константа демпфирования еще больше. ^[6] (хотя «большая» — это не совсем определенный термин, но чем больше вышеуказанные величины, тем более выраженным будет эффект жесткости.) Точное решение ( 10 ) — ${\displaystyle k}$ ${\displaystyle c}$

Уравнение 13 ведет себя очень похоже на простую экспоненту , но наличие члена , даже с малым коэффициентом, достаточно, чтобы сделать численные вычисления очень чувствительными к размеру шага. Устойчивая интеграция ( 10 ) требует очень малого размера шага до тех пор, пока не достигнет гладкой части кривой решения, что приводит к ошибке, намного меньшей, чем требуется для точности. Таким образом, система также удовлетворяет утверждению 2 и определению Ламберта. ${\displaystyle x_{0}e^{-t}}$ ${\displaystyle e^{-1000t}}$

А-стабильность

Поведение численных методов на жестких задачах можно проанализировать, применив эти методы к тестовому уравнению с начальным условием с . Решение этого уравнения — . Это решение стремится к нулю, когда Если численный метод также демонстрирует такое поведение (для фиксированного размера шага), то говорят, что метод является A-устойчивым. ^[7] Численный метод, который является L-устойчивым (см. ниже), обладает более сильным свойством, заключающимся в том, что решение стремится к нулю за один шаг, когда размер шага стремится к бесконечности. A-устойчивые методы не демонстрируют проблем неустойчивости, описанных в мотивирующем примере. ${\displaystyle y'=ky}$ ${\displaystyle y(0)=1}$ ${\displaystyle k\in \mathbb {C} }$ ${\displaystyle y(t)=e^{kt}}$ ${\displaystyle t\to \infty }$ ${\displaystyle \operatorname {Re} (k)<0.}$

Методы Рунге–Кутты

Методы Рунге–Кутты, примененные к тестовому уравнению, принимают вид , и, по индукции, . Функция называется функцией устойчивости . Таким образом, условие, что как эквивалентно . Это мотивирует определение области абсолютной устойчивости (иногда называемой просто областью устойчивости ), которая является множеством . Метод является A-устойчивым, если область абсолютной устойчивости содержит множество , то есть левую полуплоскость. ${\displaystyle y'=k\cdot y}$ ${\displaystyle y_{n+1}=\phi (hk)\cdot y_{n}}$ ${\displaystyle y_{n}={\bigl (}\phi (hk){\bigr )}^{n}\cdot y_{0}}$ ${\displaystyle \phi }$ ${\displaystyle y_{n}\to 0}$ ${\displaystyle n\to \infty }$ ${\displaystyle |\phi (hk)|<1}$ ${\displaystyle {\bigl \{}z\in \mathbb {C} \,{\big |}\,|\phi (z)|<1{\bigr \}}}$ ${\displaystyle {\bigl \{}z\in \mathbb {C} \,{\big |}\,\operatorname {Re} (z)<0{\bigr \}}}$

Пример: Методы Эйлера

Розовый диск показывает область устойчивости метода Эйлера.

Рассмотрим методы Эйлера выше. Явный метод Эйлера, примененный к тестовому уравнению, имеет вид ${\displaystyle y'=k\cdot y}$

${\displaystyle y_{n+1}=y_{n}+h\cdot f(t_{n},y_{n})=y_{n}+h\cdot (ky_{n})=y_{n}+h\cdot k\cdot y_{n}=(1+h\cdot k)y_{n}.}$

Следовательно, при . Область абсолютной устойчивости для этого метода, таким образом, представляет собой круг, изображенный справа. Метод Эйлера не является A-устойчивым. ${\displaystyle y_{n}=(1+hk)^{n}\cdot y_{0}}$ ${\displaystyle \phi (z)=1+z}$ ${\displaystyle {\bigl \{}z\in \mathbb {C} \,{\big |}\,|1+z|<1{\bigr \}}}$

Мотивирующий пример имел . Значение z при выборе размера шага равно , что находится за пределами области устойчивости. Действительно, численные результаты не сходятся к нулю. Однако при размере шага , мы имеем , что находится как раз внутри области устойчивости, и численные результаты сходятся к нулю, хотя и довольно медленно. ${\displaystyle k=-15}$ ${\displaystyle h={\tfrac {1}{4}}}$ ${\displaystyle z=-15\times {\tfrac {1}{4}}=-3.75}$ ${\displaystyle h={\tfrac {1}{8}}}$ ${\displaystyle z=-1.875}$

Пример: метод трапеций

Розовая область — область устойчивости для метода трапеций.

Рассмотрим метод трапеций.

${\displaystyle y_{n+1}=y_{n}+{\tfrac {1}{2}}h\cdot {\bigl (}f(t_{n},y_{n})+f(t_{n+1},y_{n+1}){\bigr )},}$

при применении к тестовому уравнению , ${\displaystyle y'=k\cdot y}$

${\displaystyle y_{n+1}=y_{n}+{\tfrac {1}{2}}h\cdot \left(ky_{n}+ky_{n+1}\right).}$

Решение для урожайности ${\displaystyle y_{n+1}}$

${\displaystyle y_{n+1}={\frac {1+{\frac {1}{2}}hk}{1-{\frac {1}{2}}hk}}\cdot y_{n}.}$

Таким образом, функция устойчивости имеет вид

${\displaystyle \phi (z)={\frac {1+{\frac {1}{2}}z}{1-{\frac {1}{2}}z}}}$

и область абсолютной устойчивости

${\displaystyle \left\{z\in \mathbb {C} \ \left|\ \left|{\frac {1+{\frac {1}{2}}z}{1-{\frac {1}{2}}z}}\right|<1\right.\right\}.}$

Эта область содержит левую полуплоскость, поэтому трапециевидный метод является A-устойчивым. Фактически, область устойчивости идентична левой полуплоскости, и, таким образом, численное решение сходится к нулю тогда и только тогда, когда сходится точное решение. Тем не менее, трапециевидный метод не имеет идеального поведения: он гасит все затухающие компоненты, но быстро затухающие компоненты гасятся только очень слабо, потому что при . Это привело к концепции L-устойчивости : метод является L-устойчивым, если он является A-устойчивым и при . Трапециевидный метод является A-устойчивым, но не L-устойчивым. Неявный метод Эйлера является примером L-устойчивого метода. ^[8] ${\displaystyle y'=k\cdot y}$ ${\displaystyle \phi (z)\to 1}$ ${\displaystyle z\to -\infty }$ ${\displaystyle |\phi (z)|\to 0}$ ${\displaystyle z\to \infty }$

Общая теория

Функция устойчивости метода Рунге–Кутты с коэффициентами и определяется выражением ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {b} }$

${\displaystyle \phi (z)={\frac {\det \left(\mathbf {I} -z\mathbf {A} +z\mathbf {e} \mathbf {b} ^{\mathsf {T}}\right)}{\det(\mathbf {I} -z\mathbf {A} )}},}$

где обозначает вектор со всеми единицами. Это рациональная функция (один многочлен , деленный на другой). ${\displaystyle \mathbf {e} }$

Явные методы Рунге–Кутты имеют строго нижнюю треугольную матрицу коэффициентов , и, таким образом, их функция устойчивости является полиномом. Из этого следует, что явные методы Рунге–Кутты не могут быть A-устойчивыми. ${\displaystyle \mathbf {A} }$

Функция устойчивости неявных методов Рунге–Кутты часто анализируется с использованием звезд порядка. Звезда порядка для метода с функцией устойчивости определяется как множество . Метод является A-устойчивым тогда и только тогда, когда его функция устойчивости не имеет полюсов в левой плоскости, а его звезда порядка не содержит чисто мнимых чисел. ^[9] ${\displaystyle \phi }$ ${\displaystyle {\bigl \{}z\in \mathbb {C} \,{\big |}\,|\phi (z)|>|e^{z}|{\bigr \}}}$

Многошаговые методы

Линейные многошаговые методы имеют вид

${\displaystyle y_{n+1}=\sum _{i=0}^{s}a_{i}y_{n-i}+h\sum _{j=-1}^{s}b_{j}f\left(t_{n-j},y_{n-j}\right).}$

Применительно к тестовому уравнению они становятся

${\displaystyle y_{n+1}=\sum _{i=0}^{s}a_{i}y_{n-i}+hk\sum _{j=-1}^{s}b_{j}y_{n-j},}$

что можно упростить до

${\displaystyle \left(1-b_{-1}z\right)y_{n+1}-\sum _{j=0}^{s}\left(a_{j}+b_{j}z\right)y_{n-j}=0}$

где . Это линейное рекуррентное соотношение . Метод является A-устойчивым, если все решения рекуррентного соотношения сходятся к нулю при . Характеристический многочлен равен ${\displaystyle z=hk}$ ${\displaystyle \{y_{n}\}}$ ${\displaystyle \operatorname {Re} (z)<0}$

${\displaystyle \Phi (z,w)=w^{s+1}-\sum _{i=0}^{s}a_{i}w^{s-i}-z\sum _{j=-1}^{s}b_{j}w^{s-j}.}$

Все решения сходятся к нулю для заданного значения , если все решения лежат в единичной окружности. ${\displaystyle z}$ ${\displaystyle w}$ ${\displaystyle \Phi (z,w)=0}$

Область абсолютной устойчивости для многошагового метода вышеуказанной формы — это множество всех, для которых все такие, что удовлетворяют . Опять же, если это множество содержит левую полуплоскость, многошаговый метод называется A-устойчивым. ${\displaystyle z\in \mathbb {C} }$ ${\displaystyle w}$ ${\displaystyle \Phi (z,w)=0}$ ${\displaystyle |w|<1}$

Пример: Метод Адамса–Башфорта второго порядка

Розовая область — область устойчивости метода Адамса–Башфорта второго порядка.

Определим область абсолютной устойчивости для двухшагового метода Адамса–Башфорта

${\displaystyle y_{n+1}=y_{n}+h\left({\tfrac {3}{2}}f(t_{n},y_{n})-{\tfrac {1}{2}}f(t_{n-1},y_{n-1})\right).}$

Характеристический многочлен равен

${\displaystyle \Phi (w,z)=w^{2}-\left(1+{\tfrac {3}{2}}z\right)w+{\tfrac {1}{2}}z=0}$

который имеет корни

${\displaystyle w={\tfrac {1}{2}}\left(1+{\tfrac {3}{2}}z\pm {\sqrt {1+z+{\tfrac {9}{4}}z^{2}}}\right),}$

таким образом, область абсолютной устойчивости

${\displaystyle \left\{z\in \mathbb {C} \ \left|\ \left|{\tfrac {1}{2}}\left(1+{\tfrac {3}{2}}z\pm {\sqrt {1+z+{\tfrac {9}{4}}z^{2}}}\right)\right|<1\right.\right\}.}$

Эта область показана справа. Она не включает всю левую полуплоскость (фактически она включает только действительную ось между ), поэтому метод Адамса–Башфорта не является A-устойчивым. ${\displaystyle -1\leq z\leq 0}$

Общая теория

Явные многошаговые методы никогда не могут быть A-устойчивыми, как и явные методы Рунге–Кутты. Неявные многошаговые методы могут быть A-устойчивыми только в том случае, если их порядок не превышает 2. Последний результат известен как второй барьер Дальквиста ; он ограничивает полезность линейных многошаговых методов для жестких уравнений. Примером метода второго порядка A-устойчивого является упомянутое выше правило трапеций, которое также можно рассматривать как линейный многошаговый метод. ^[10]

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

• Формула обратной дифференциации , семейство неявных методов, используемых в основном для решения жестких дифференциальных уравнений.
• Номер состояния
• Дифференциальное включение , расширение понятия дифференциального уравнения, допускающее разрывы, отчасти как способ обойти некоторые проблемы жесткости.
• Явные и неявные методы

Примечания

1. Робертсон, ХХ (1966). «Решение набора уравнений скорости реакции». Численный анализ: введение . Academic Press. С. 178–182.
2. Ламберт (1992, стр. 216–217)
3. Ламберт (1992, стр. 217–220)
4. Хиршфельдер (1963)
5. Бремя и ярмарки (1993, стр. 314)
6. Крейциг (1972, стр. 62–68)
7. Это определение принадлежит Далквисту (1963).
8. Определение L-стабильности принадлежит Эле (1969).
9. Определение принадлежит Ваннеру, Хайреру и Нёрсетту (1978); см. также Исерлес и Нёрсетт (1991).
10. См. Далквист (1963).

Ссылки

• Берден, Ричард Л.; Фейрес, Дж. Дуглас (1993), Численный анализ (5-е изд.), Бостон: Prindle, Weber and Schmidt, ISBN 0-534-93219-3.
• Дальквист, Джермунд (1963), «Специальная проблема устойчивости для линейных многошаговых методов», BIT , 3 (1): 27–43, doi : 10.1007/BF01963532, hdl : 10338.dmlcz/103497 , S2CID  120241743.
• Эберли, Дэвид (2008), Анализ устойчивости систем дифференциальных уравнений (PDF).
• Эле, Б. Л. (1969), О приближениях Паде к экспоненциальной функции и A-устойчивых методах для численного решения задач с начальными значениями (PDF) , Университет Ватерлоо.
• Gear, CW (1971), Численные задачи начального значения в обыкновенных дифференциальных уравнениях , Englewood Cliffs: Prentice Hall , Bibcode : 1971nivp.book.....G.
• Gear, CW (1981), «Численное решение обыкновенных дифференциальных уравнений: осталось ли что-то еще сделать?», SIAM Review , 23 (1): 10–24, doi :10.1137/1023002.
• Хайрер, Эрнст; Ваннер, Герхард (1996), Решение обыкновенных дифференциальных уравнений II: Жесткие и дифференциально-алгебраические проблемы (второе изд.), Берлин: Springer-Verlag , ISBN 978-3-540-60452-5.
• Хиршфельдер, Дж. О. (1963), «Прикладная математика, используемая в теоретической химии», Симпозиум Американского математического общества : 367–376.
• Изерлес, Арье; Норсетт, Сиверт (1991), Звезды Ордена , Чепмен и Холл , ISBN 978-0-412-35260-7.
• Крейциг, Эрвин (1972), Advanced Engineering Mathematics (3-е изд.), Нью-Йорк: Wiley , ISBN 0-471-50728-8.
• Ламберт, Дж. Д. (1977), Д. Якобс (ред.), «Задача начального значения для обыкновенных дифференциальных уравнений», The State of the Art in Numerical Analysis , Нью-Йорк: Academic Press : 451–501.
• Ламберт, Дж. Д. (1992), Численные методы для обыкновенных дифференциальных систем , Нью-Йорк: Wiley , ISBN 978-0-471-92990-1.
• Мэтьюз, Джон; Финк, Куртис (1992), Численные методы с использованием MATLAB.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Раздел 17.5. Жесткие наборы уравнений". Numerical Recipes: The Art of Scientific Computing (3-е изд.). Нью-Йорк: Cambridge University Press. ISBN 978-0-521-88068-8. Архивировано из оригинала 2011-08-11 . Получено 2011-08-17 .
• Шэмпайн, Л. Ф.; Гир, К. В. (1979), «Взгляд пользователя на решение жестких обыкновенных дифференциальных уравнений», SIAM Review , 21 (1): 1–17, doi :10.1137/1021001.
• Ваннер, Герхард; Хайрер, Эрнст; Норсетт, Сиверт (1978), «Звезды порядка и теория стабильности», BIT , 18 (4): 475–489, doi : 10.1007/BF01932026, S2CID  8824105.
• Устойчивость методов Рунге-Кутты [1]

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

• Введение в физически обоснованное моделирование: энергетические функции и жесткость
• Жесткие системы Лоуренс Ф. Шампайн и Скип Томпсон Scholarpedia , 2(3):2855. doi:10.4249/scholarpedia.2855