Исследование схем на устойчивость

Канонический вид двухслойной схемы

Пусть имеем произвольную задачу Коши с любой старшей производной. Из курса диффуров знаем, что её всегда можно переписать в виде системы:

\[ \frac{d \mathbf{u}}{d t}=\mathbf{f}(t, \mathbf{u}), \quad t>0, \quad \mathbf{u}(0)=\mathbf{u_0} , \quad \mathbf{u}: \mathbb{R} \rightarrow \mathbb{R}^N , \quad \mathbf{f} : \mathbb{R}^{N+1} \rightarrow \mathbb{R}^N \]

которую формально можно записать в операторном виде:

\[ \hat{L} \; \mathbf{u}(t)=\mathbf{f}(t, \mathbf{u}) \]

и поставить ей в соответствие разностную схему (дискретизация):

\[ \hat{L}_{h} \mathbf{y}^{(h)}=\mathbf{F}^{(h)} \]

Пусть эта схема двухслойная (т.е. \(\mathbf{y}_{n+1}\) может зависеть только от \(\mathbf{y}_{n}\)).

Определение. Каноническая форма двухслойной схемы - это схема в виде:

\[ \mathbf{y}_{n+1}=\hat{R}_h \cdot \mathbf{y}_n+h \cdot \mathbf{\rho}_n \]

где явно выделен так называемый оператор послойного перехода \(\hat{R}_h\). Естественно, этот оператор должен быть линеен, и действует только на \(\mathbf{y}_n\).

Каноническую форму можно всегда вывести из обычного вида численной схемы.

Для решения разностной задачи каноническая форма, в общем-то, бесполезна, т.к. в общем случае для её нахождения необходимо линеаризировать истинное решение диффура. Тем не менее, оператор оператор перехода \(\hat{R}_h\) не будет зависеть в явном виде от истинного решения, что даёт нам возможность анализировать через него систему на устойчивость.

Пример. Найдём канонический вид неявной схемы Эйлера.

Пусть имеем ОДУ:

\[u^{\prime}=f(t, u), \quad u(0)=u_0\]

Неявная схема Эйлера имеет вид:

\[ \frac{y_{n+1}-y_n}{\tau}=f\left(t_{n+1}, y_{n+1}\right), \quad y_0=u_0 . \]

Перепишем разностную схему в виде

\[ y_{n+1}=y_n+\tau f\left(t_{n+1}, y_{n+1}\right) \]

Хотим избавиться от \(y_{n+1}\) в нелинейной в общем случае \(f\). Пусть \(\varphi(t)=y(t)\) - истинное (математическое) решение. Для малых \(\tau\) отклонение от истинной траектории будет невилико. Таким образом, раскладывая вблизи истинной траектории:

\[ f\left(t_{n+1}, y_{n+1}\right) \approx f\left(t_{n+1}, \; \varphi \left(t_{n} \right)\right)+\frac{\partial f}{\partial u}\left(y_{n+1}-\varphi(t_{n})\right) \]

А тогда канонический вид:

\[ y_{n+1}=\hat{R}_\tau y_n+\tau \rho_n \]

\[ \hat{R}_\tau=\left(1-\tau \left.\frac{\partial f}{\partial u}\right|_{(y_n, \, t_n)}\right)^{-1}, \quad \rho_n=f\left(t_{n+1}, \varphi(t_{n})\right)-\left.\frac{\partial f}{\partial u}\right|_{(y_n,\, t_n)} \varphi(t_{n}) \]

Как и было обещано, \(\hat{R}_\tau\) не зависит от истинного решения \(\varphi(t)\), но зависит от локального положения траектории.

Также отметим, что в общем случае (системы ОДУ), \(\frac{\partial f}{\partial u}\) - матрица Якоби и \(\hat{R}_\tau\) является именно матрицей.

Исследование устойчивости. Спектральный признак Фон Неймана

Теорема (достаточное условие устойчивости, спектральный признак Фон Неймана). Пусть имеем записанную в каноническом виде систему \(\mathbf{y}_{n+1}=\hat{R}_h \cdot \mathbf{y}_n+h \cdot \mathbf{\rho}_n\). Тогда для её корректности (а значит и устойчивости) на любом отрезке интегрирования достаточно выполнения в каждой точке траектории

\[ \exists C_2: \forall n \; \hookrightarrow \;\left\|\mathbf{\rho}_n\right\| \leq C_2\|\mathbf{f}\| \; \; \text{и} \; \; \left\|\mathbf{y}_0\right\| \leq C_2\|\mathbf{f}\| \]

\[ \forall \lambda_i \hookrightarrow \left\|\lambda_i\right\| \leq 1 \]

где \(\lambda_i\) - собственные числа матрицы послойного перехода \(\hat{R}_h\).

Примечание. Условие в первой строчке должны быть выполнены в каждой точке траектории, т.е. на практике проверяются пост-фактум.

Доказательство спектрального признака Фон Неймана

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

Лемма. Пусть разностная схема \(\hat{L}_\tau y^{(\tau)}=f^{(\tau)}\) приведена к каноническому виду \(y_{n+1}=\hat{R}_\tau \cdot y_n+\tau \cdot \rho_n\), и выполнены неравенства:

\[ \left\|y_0\right\| \leq C_2\|\mathbf{f}\| \]

\[ \left\|\rho_n\right\| \leq C_2\|\mathbf{f}\| \]

Тогда для корректности \(\left\|y^{(\tau)}\right\| \leq C\|f^{(\tau)}\|\) (а значит и для устойчивости) достаточно, чтобы нормы степеней оператора \(\left\|\hat{R}_\tau^m\right\|\) были равномерно по \(\tau\) ограничены, т.е. чтобы выполнялась оценка

\[ \exists C_3 > 0: \; \forall \tau \in [0, \; \tau_0], \;\; \forall m \geq 1 \; \hookrightarrow \; \left\|\hat{R}_\tau^m\right\| \leq C_3 \]

При этом в качестве числа С может быть взята величина (\(T\) - общая длина отрезка вычисляемой траектории):

\[ C=(1+T) \cdot C_2 C_3 \]

Откуда видим, что чем ближе \(C_3\) к 1, тем лучше устойчивость.

Доказательство спектрального признака Фон Неймана

Для проверки выполнения требования

\[\left\|\hat{R}_\tau^m\right\| \leq C_3\]

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

Для нормы оператора справедливо неравенство (см. семинар 1.4):

\[ \left\|\hat{R}_\tau\right\| \leq \max _i\left|\lambda_i\right| \]

Тогда для выполнения требования \(\left\|\hat{R}_\tau^m\right\| \leq C_3\) необходимо, чтобы все собственные значения оператора послойного перехода (которые являются функциями чисел \(t, u\), т.к. матрица \(\hat{R}_\tau\) сама зависит от \(t, u\)) лежали в круге (см. семинар 3.3):

\[ \exists c\geq0: \; \forall i, t, u \; \hookrightarrow \; \left|\lambda_i\right| \leq 1+c \tau \]

при этом постоянная \(c\) не должна зависеть от сеточных параметров. Математически, можем определить \(c\) как:

\[ c = \max{\left(\cfrac{\max_{i, t, \mathbf{u}} \left|\lambda_i\right| - 1}{\tau}, \;\; 0\right)} \]

А тогда при выполнении условия выше имеем (\(N = \frac{T}{\tau}\) - общее количество шагов интегрирования):

\[ \left\|\hat{R}_\tau^m\right\| \leq \max _i\left|\lambda_i\right|^m \leq |1+c \tau|^m \leq e^{c \tau \cdot m} \leq e^{c \tau N} = e^{cT} = C_3 \]

Определение. Если получается выбрать только \(c>0\) (что эквивалентно \(\exists \lambda_i : \left\| \lambda_i \right\| >1 \)), то устойчивость называется нестрогой, так как можем обеспечить устойчивость только на на коротком интервале \(T\) из условия малости \(c \cdot T \approx 0\). Более того, устремление \(\tau \rightarrow 0\) даёт \(c \rightarrow +\infty\), что отвечает бесконечному росту \(\left\|\hat{R}_\tau^m\right\|\), что делает нашу оценку бессмысленной - надо делать дополнительный анализ.

Определение. Если можно выбрать \(c=0\) (что эквивалентно \(\forall \lambda_i \hookrightarrow \left\| \lambda_i \right\| \leq1 \)), то устойчивость называется строгой. В этом случае устойчивость формально сохраняется при бесконечном \(T\). Это необходимо при расчетах диффуров на неограниченных (заранее неизвестных какой длины) интервалах и возможности устремления \(\tau \rightarrow 0\). Более того, только в этом случае при \(\tau \rightarrow 0\) оценка выше доказывает устойчивость. Это и есть достаточное условие устойчивости.

Явные методы, как правило, условно устойчивы - т.е. мы не можем до бесконечности понижать шаг \(\tau\). Неявные методы лишены этого недостатка, однако при их использовании приходится решать систему алгебраических уравнений, вообще говоря, нелинейную.

Вывод. Если со строгой устойчивостью всё понятно - тут всё хорошо, то с нестрогой нужно делать дополнительный анализ. Главное из этого пункта уяснить, что собственные числа матрицы послойного перехода \(\hat{R}_\tau\) связаны со сходимостью.

Исследование устойчивости. Алгоритм Далквиста

Иной подход к анализу устойчивости даёт метод Далквиста. Он основан на исследовании нашего метода дискретизации для модельного семейства систем ОДУ.

Пусть имеем произвольную задачу Коши (записанную в виде системы):

\[ \frac{d \mathbf{u}}{d x}=\mathbf{f}(x, \mathbf{u}), \quad x>0, \quad \mathbf{u}(0)=\mathbf{u_0} , \quad \mathbf{u}: \mathbb{R} \rightarrow \mathbb{R}^N , \quad \mathbf{f} : \mathbb{R^{N+1}} \rightarrow \mathbb{R}^N \]

Поставим ей в соответствие разностную схему по какому-то алгоритму дискретизации:

\[ \hat{L}_{h} \mathbf{y}^{(h)}=\mathbf{F}^{(h)} \]

Хотим проверить, устойчив ли наш алгоритм дискретизации. Один из способов - использовать уравнение Далквиста (явно подставили \(\mathbf{f}(x, \mathbf{u}) =\lambda \mathbf{u}\)):

\[ \mathbf{u}^{\prime}=\lambda \mathbf{u} \]

Здесь константа \(\lambda\) играет роль какого-либо собственного значения матрицы Якоби исходной системы ОДУ и может принимать комплексные значения.

Дискритизованная версия по ранее определенному алгоритму, записанная в каноническом виде, будет выглядеть как:

\[ \mathbf{y}_{n+1}=\hat{R}_h \cdot \mathbf{y}_n+h \cdot {\mathbf{\rho}_n}=\hat{R}_h \cdot \mathbf{y}_n \]

\[ \mathbf{y}_{n+1}=\hat{R}_h \cdot \mathbf{y}_n \]

Т.е. в случае рассмотрения уравнения Далквиста, часть \(\mathbf{\rho}_n\) всегда зануляется из-за линейности \(\mathbf{f}(x, \mathbf{u}) =\lambda \mathbf{u}\). Поэтому уравнение Далквиста и классное.

Более того, оператор послойного перехода будет явной функцией от \(z = \lambda \cdot h \in \mathbb{C}\) (т.к. \(\frac{\mathbf{du}}{dx}=\lambda \mathbf{u} \Leftrightarrow \frac{\mathbf{du}}{d(\lambda x)}= \mathbf{u}\)). То есть можем обозначить \(\hat{R}_h = R(z)\). Тогда в итоге можем записать,

\[ \mathbf{y}_{n+1}=R(z) \mathbf{y}_n, \quad z=\lambda h \]

где \(R(z)\) теперь называется разрешающим оператором и одновременно функцией устойчивости. Стоит различать функцию учтойчивости и оператор послойного перехода, т.к. первое вводится только для уравнения Далквиста. В духе спектрального признака Фон Неймана разумно ввести следующее определение:

Определение 1. Множество \(S=\{z=h \lambda \in \mathbb{C}, \left\|R(z)\right\| \leq 1\}\) называется областью (абсолютной) устойчивости данного метода. Для понимания, почему эта область важна, см. доказательство спектрального признака Фон Неймана выше.

Примечание. На практике, устойчивость в конкретной точке траектории определяется из условия \(h \cdot \lambda (y_n, t_n) \in S\), где \(\lambda (y_n, t_n)\) - все собственные числа матрицы Якоби в данной точке.

Найдём область убсолютной устойчивости какого-нибудь метода.

Пример 1. Область устойчивости явного метода Эйлера.

Напишем явный метод Эйлера:

\[ \frac{y_{n+1}-y_n}{\tau}=f(t_n, y_n), \quad y_0=a, \quad n=0,1, \ldots, \]

В контексте уравнения Далквиста (\(f=\lambda y\))

\[ \frac{y_{n+1}-y_n}{\tau}=\lambda y_n, \quad y_0=a, \quad n=0,1, \ldots, \]

откуда следует

\[ y_{n+1}=(1+\tau \lambda) y_n \]

Т.е. функция устойчивости имеет вид

\[ R(z)=1+z \]

Нарисуем область устойчивости:

from sympy import symbols, Abs
from sympy import I
from sympy.plotting import plot_implicit

# Define the variable z
z = symbols('z', complex=True)
x, y = symbols('x y', real=True)
# Define the function R
R = 1 + z

# Define the absolute value of R
abs_R = Abs(R)

# Define the condition for the set |R(z)|<1
cond = abs_R < 1
display(cond)

z_number = x + I*y
cond = cond.subs(z, z_number)

display(cond)
# Plot the set where |R(z)|<1

plot_implicit(cond, adaptive = False,
              title='Область устойчивости явной схемы Эйлера',
              xlabel='Real part of z', ylabel='Imaginary part of z')

\[\displaystyle \left|{z + 1}\right| < 1\]

\[\displaystyle \sqrt{x^{2} + 2 x + y^{2} + 1} < 1\]

<sympy.plotting.plot.Plot at 0x7feb48d4e430>

Т.е. явная схема Эйлера устойчива только для не слишком быстро затухающих процессов с «несильными» осцилляциями.

Пример 2. Область устойчивости неявного метода Эйлера.

Напишем неявный метод Эйлера:

\[ \frac{y_{n+1}-y_n}{\tau}=f(t_{n+1}, y_{n+1}), \quad y_0=a, \quad n=0,1, \ldots, \]

В контексте уравнения Далквиста (\(f=\lambda y\))

\[ \frac{y_{n+1}-y_n}{\tau}=\lambda y_{n+1}, \quad y_0=a, \quad n=0,1, \ldots, \]

откуда следует

\[ y_{n+1}=\cfrac{y_n}{(1-\tau \lambda)} \]

Т.е. функция устойчивости имеет вид

\[ R(z)=\cfrac{1}{1-z} \]

Нарисуем область устойчивости:

from sympy import symbols, Abs
from sympy import I
from sympy.plotting import plot_implicit

# Define the variable z
z = symbols('z', complex=True)
x, y = symbols('x y', real=True)
# Define the function R
R = 1/(1 - z)

# Define the absolute value of R
abs_R = Abs(R)

# Define the condition for the set |R(z)|<1
cond = abs_R < 1
display(cond)

z_number = x + I*y
cond = cond.subs(z, z_number)

display(cond.doit())
# Plot the set where |R(z)|<1

plot_implicit(cond, adaptive = False,
              title='Область устойчивости неявной схемы Эйлера',
              xlabel='Real part of z', ylabel='Imaginary part of z')

\[\displaystyle \frac{1}{\left|{z - 1}\right|} < 1\]

\[\displaystyle \frac{1}{\sqrt{x^{2} - 2 x + y^{2} + 1}} < 1\]

<sympy.plotting.plot.Plot at 0x7feb48c880d0>

То есть неявная схема Эйлера устойчива для всего, кроме слабо растущих процессов с несильными осцилляциями.

Классификация функций устойчивости

Определение 2 (Далквист, 1963). Если область абсолютной устойчивости \(S\) включает в себя левую полуплоскость комплексной плоскости \(\operatorname{Re}(z) \leq 0\), то метод называется А-устойчивым.

Множество \(\operatorname{Re}(z) \leq 0\) является существенно важным в вычислительных методах, потому что именно ему отвечают экспонциально затухающие процессы, с которыми мы как правило и работаем (т.к с экспонциально растущими процессами всё как правило и так плохо, то если ещё и с быстрозатухающими процессами всё плохо, то всё совсем плохо-плохо).

Определение 3. Если область абсолютной устойчивости \(S\) включает в себя угол с полураствором \(\alpha\), отсчитываемым от отрицательного направления действительной оси в плоскости \(z\), то метод называется \(\mathrm{A}(\alpha)\)-устойчивым.

Метод является \(\mathrm{A}(0)-\) устойчивым, если условие верно только при \(\alpha \rightarrow 0\).

Определение 4. Численный метод называется L-устойчивым (или асимптотически устойчивым), если он А-устойчив и выполнено условие \(|R(z)| \rightarrow 0\) при \(\operatorname{Re} z \rightarrow-\infty\).

Иногда в литературе можно встретить определение асимптотической устойчивости, не включающее в себя А-устойчивости. В этом случае численный метод должен быть, например, жестко устойчивым. Понятие жесткой устойчивости строго не определяется, но подразумевается, что для жесткой устойчивости метод будет абсолютно устойчивым на всем спектре рассматриваемой задачи. Кроме того. существует и альтернативная теория устойчивости численных методов для решения ЖС ОДУ.

Определение 5. Численный метод называется \(Lp\)-устойчивым, если он А-устойчив и функция устойчивости затухает как \(z^{-p}\) при \(\operatorname{Re} z \rightarrow-\infty\).

Решения уравнения Далквиста, полученные L-устойчивыми методами, будут затухающими.

Еще одним важным свойством разностных схем является их монотонность. Говорить о монотонности численного метода для решения систем ОДУ можно только в том случае, когда решение дифференциальной задачи также монотонно, т.е. только для действительных отрицательных \(z\) и, таким образом, для строго убывающих решений.

Определение 6. Если функция устойчивости положительна на отрицательной части действительной оси (т.е. \(R(z)>0\) при \(\operatorname{Im} z=0, \operatorname{Re} z<0\) ), то метод называется монотонным.

Теорема (барьер Далквиста).

Теорема (барьер Далквиста). Не существует A-устойчивых многошаговых разностных схем порядка аппроксимации выше второго.

То есть, самым важным классом численных схем являются именно одношаговые методы, т.к. для них возможна (но вовсе необязательна) А-устойчивость.

Для повышения порядка аппроксимации с использованием многошаговых схем приходится мириться с отсутствием А-устойчивости и использовать вспомогательные понятия, как A(\(\alpha\))-устойчивость. Такие алгоритмы могут иметь порядок аппроксимации выше 2-го, но очень чувствительны к малым осцилляциям определённых частот.