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

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

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

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

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

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

${\displaystyle \,y'(t)=-15y(t),\quad t\geq 0,y(0)=1.}$

(1)

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

${\displaystyle y(t)=e^{-15t}\,}$ с ${\displaystyle y(t)\to 0}$ в качестве ${\displaystyle t\to \infty .}$

(2)

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

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

1. Метод Эйлера с размером шага час = 1/4 сильно колеблется и быстро выходит за пределы графика (показано красным).
2. Метод Эйлера с половиной шага, час = 1/8, дает решение в границах графика, но колеблется около нуля (показано зеленым).
3. В трапециевидный метод (то есть двухступенчатый Метод Адамса – Моултона ) дан кем-то

   ${\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}h\left(f(t_{n},y_{n})+f(t_{n+1},y_{n+1})\right),}$

   (3)

   куда ${\displaystyle \textstyle y'=f(t,y)}$. Применение этого метода вместо метода Эйлера дает гораздо лучший результат (синий). Численные результаты монотонно убывают до нуля, как и точное решение.

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

${\displaystyle {\dot {x}}=-0.04x+10^{4}y\cdot z}$ ${\displaystyle {\dot {y}}=0.04x-10^{4}y\cdot z-3\cdot 10^{7}y^{2}}$ ${\displaystyle {\dot {z}}=3\cdot 10^{7}y^{2}}$

(4)

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

Дополнительными примерами являются наборы ODE, полученные в результате временной интеграции крупных механизмов химических реакций. Здесь жесткость возникает из-за сосуществования очень медленных и очень быстрых реакций.^{[нужна цитата ]} Для их решения программные пакеты КПП и Автохим может быть использован.

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

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

${\displaystyle \mathbf {y} '=\mathbf {A} \mathbf {y} +\mathbf {f} (x),}$

(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}$. Общее решение (5) принимает вид

${\displaystyle \mathbf {y} (x)=\sum _{t=1}^{n}\kappa _{t}\exp(\lambda _{t}x)\mathbf {c} _{t}+\mathbf {g} (x),}$

(6)

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

${\displaystyle Re(\lambda _{t})<0,\qquad t=1,2,\ldots ,n,}$

(7)

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

${\displaystyle |Re({\overline {\lambda }})|\geq |Re(\lambda _{t})|\geq |Re({\underline {\lambda }})|,\qquad t=1,2,\ldots ,n}$

(8)

так что ${\displaystyle \kappa _{t}\exp({\overline {\lambda }}x)\mathbf {c} _{t}}$ самый быстрый переходный и ${\displaystyle \kappa _{t}\exp({\underline {\lambda }}x)\mathbf {c} _{t}}$ самый медленный. Теперь определим коэффициент жесткости в качестве

${\displaystyle {\frac {|Re({\overline {\lambda }})|}{|Re({\underline {\lambda }})|}}.}$[2]

(9)

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

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

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

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

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

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

Этимология

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

Значительные трудности могут возникнуть при стандартном численные методы применяются для аппроксимации решения дифференциальное уравнение когда точное решение содержит члены вида е^λt, где λ - комплексное число с отрицательной действительной частью.

...

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

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

${\displaystyle m{\ddot {x}}+c{\dot {x}}+kx=0,\qquad x(0)=x_{0},\qquad {\dot {x}}(0)=0,}$

(10)

с м = 1, c = 1001, k = 1000, можно записать в виде (5) с п = 2 и

${\displaystyle \mathbf {A} =\left({\begin{array}{rr}0&1\\-1000&-1001\end{array}}\right),}$

(11)

${\displaystyle \mathbf {f} (t)=\left({\begin{array}{c}0\\0\end{array}}\right),}$

(12)

${\displaystyle \mathbf {x} (0)=\left({\begin{array}{c}x_{0}\\0\end{array}}\right),}$

(13)

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

${\displaystyle {\frac {|-1000|}{|-1|}}=1000,}$

(14)

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

${\displaystyle x(t)=x_{0}\left(-{\frac {1}{999}}e^{-1000t}+{\frac {1000}{999}}e^{-t}\right)\approx x_{0}e^{-t}.}$

(15)

Обратите внимание, что (15) ведет себя почти как простая экспоненциальная Икс₀е^−т, но наличие е^−1000т Даже с небольшим коэффициентом достаточно, чтобы числовые вычисления были очень чувствительны к размеру шага. Стабильная интеграция (10) требует очень небольшого шага, пока он не войдет в гладкую часть кривой решения, что приведет к ошибке, намного меньшей, чем требуется для точности. Таким образом, система также удовлетворяет утверждению 2 и определению Ламберта.

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

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

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

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

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

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

Рассмотрим методы Эйлера выше. Явный Метод Эйлера применительно к испытательному уравнению ${\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}.}$

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

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

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

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

Рассмотрим трапециевидный метод

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

применительно к испытательному уравнению ${\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)={1+{1 \over 2}z \over 1-{1 \over 2}z}}$

а область абсолютной стабильности

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

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

Общая теория

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

${\displaystyle \phi (z)={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},}$

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

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

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

Многоступенчатые методы

Линейные многоступенчатые методы иметь форму

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

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

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

который можно упростить до

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

куда z = гонконгский. Это линейный отношение повторения. Метод является A-устойчивым, если все решения {у_п} рекуррентного соотношения сходятся к нулю, когда 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}.}$

Все решения сходятся к нулю при заданном значении z если все решения ш из Φ (z,ш) = 0 лежат в единичной окружности.

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

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

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

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

${\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}-(1+{\tfrac {3}{2}}z)w+{\tfrac {1}{2}}z=0}$

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

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

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

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

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

Общая теория

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

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

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

Примечания

1. Робертсон, Х. Х. (1966). «Решение системы уравнений скорости реакции». Численный анализ: введение. Академическая пресса. С. 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-е изд.), Бостон: Принл, Вебер и Шмидт, ISBN  0-534-93219-3.
• Дальквист, Гермунд (1963), «Специальная задача устойчивости для линейных многоступенчатых методов», КУСОЧЕК, 3 (1): 27–43, Дои:10.1007 / BF01963532, HDL:10338.dmlcz / 103497.
• Эберли, Дэвид (2008), Анализ устойчивости систем дифференциальных уравнений (PDF).
• Эле, Б. Л. (1969), Об аппроксимациях Паде экспоненциальной функции и A-устойчивых методах численного решения начальных задач (PDF), Университет Ватерлоо.
• Гир, К. У. (1971), Численные начальные задачи в обыкновенных дифференциальных уравнениях, Энглвудские скалы: Prentice Hall.
• Gear, C. W. (1981), "Численное решение обыкновенных дифференциальных уравнений: что-нибудь еще нужно сделать?", SIAM Обзор, 23 (1): 10–24, Дои:10.1137/1023002.
• Хайрер, Эрнст; Ваннер, Герхард (1996), Решение обыкновенных дифференциальных уравнений II: жесткие и дифференциально-алгебраические задачи (второе изд.), Берлин: Springer-Verlag, ISBN  978-3-540-60452-5.
• Хиршфельдер, Дж. О. (1963), «Прикладная математика, используемая в теоретической химии», Симпозиум Американского математического общества: 367–376.
• Изерлес, Арье; Норсетт, Syvert (1991), Заказать звезды, Чепмен и Холл, ISBN  978-0-412-35260-7.
• Крейсциг, Эрвин (1972), Высшая инженерная математика (3-е изд.), Нью-Йорк: Wiley, ISBN  0-471-50728-8.
• Ламберт, Дж. Д. (1977), Д. Джейкобс (ред.), "Задача начального значения для обыкновенных дифференциальных уравнений", Современное состояние численного анализа, Нью-Йорк: Академическая пресса: 451–501.
• Ламберт, Дж. Д. (1992), Численные методы для обыкновенных дифференциальных систем., Нью-Йорк: Wiley, ISBN  978-0-471-92990-1.
• Мэтьюз, Джон; Финк, Куртис (1992), Численные методы с использованием MATLAB.
• Нажмите, WH; Теукольский С.А.; Феттерлинг, штат Вашингтон; Фланнери, ВР (2007). «Раздел 17.5. Жесткие множества уравнений». Числовые рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк: Издательство Кембриджского университета. ISBN  978-0-521-88068-8.
• Шампин, Л. Ф .; Gear, C. W. (1979), «Взгляд пользователя на решение жестких обыкновенных дифференциальных уравнений», SIAM Обзор, 21 (1): 1–17, Дои:10.1137/1021001.
• Ваннер, Герхард; Хайрер, Эрнст; Норсетт, Сиверт (1978), "Звезды порядка и теория устойчивости", КУСОЧЕК, 18 (4): 475–489, Дои:10.1007 / BF01932026.
• Устойчивость методов Рунге-Кутты. [1]

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

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