Сходимость разностных схем

Сходимость = Аппроксимация + Устойчивость

Теорема. (Лакса-Рябенького-Филиппова). Пусть разностная задача аппроксимирует дифференциальную задачу на ее решении и устойчива. Тогда решение разностной задачи сходится к решению дифференциальной. Если аппроксимация имеет порядок \(k\) по времени и порядок р по пространству, то и сходимость имеет те же порядки.

Анализ аппроксимации

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

Пример 1. Для разнообразия рассмотрим УРЧП первого порядка.

\[ \frac{\partial u}{\partial t} = \frac{\partial u}{\partial x} + f(x, t) \]

относительно следующей разностной задачи (явная четырёхточечная схема Эйлера):

\[ \frac{u^{n+1}_j-u^{n}_j}{\tau} = \frac{u^n_{j+1}-u^n_{j-1}}{2\cdot h} + f^n_j \]

Разложение в ряд Тейлора удобнее всего делать с помощью Wolfram или sympy. Будьте внимательны, в sympy нет надежного способа разложить в ряд Тейлора функцию от нескольких переменных, так что он не подойдёт для всех случаев.

import sympy as sp

# Define the symbolic variables
f, x, t, h, tau = sp.symbols('f x t h tau')

# Define the symbolic function u(x,t)
u = sp.Function('u')(x,t)
f = sp.Function('f')(x,t)


# Define the finite difference scheme for u_t = u_x + f
ut = (u.subs(t, t+tau) - u) / tau
ux = (u.subs(x, x+h) - u.subs(x, x-h)) / (2*h)

# Define the Taylor series expansion of u(x+h,t+tau) and u(x-h,t+tau)
ut_series = sp.series(ut, tau, x0=0, n=4)
ux_series = sp.series(ux, h, x0=0, n=4)
ut_series - ux_series - f

\[\displaystyle \left. \frac{\partial}{\partial \xi_{2}} u{\left(x,\xi_{2} \right)} \right|_{\substack{ \xi_{2}=t }} - \left. \frac{\partial}{\partial \xi_{1}} u{\left(\xi_{1},t \right)} \right|_{\substack{ \xi_{1}=x }} - f{\left(x,t \right)} + \frac{\tau \left. \frac{\partial^{2}}{\partial \xi_{2}^{2}} u{\left(x,\xi_{2} \right)} \right|_{\substack{ \xi_{2}=t }}}{2} + \frac{\tau^{2} \left. \frac{\partial^{3}}{\partial \xi_{2}^{3}} u{\left(x,\xi_{2} \right)} \right|_{\substack{ \xi_{2}=t }}}{6} + \frac{\tau^{3} \left. \frac{\partial^{4}}{\partial \xi_{2}^{4}} u{\left(x,\xi_{2} \right)} \right|_{\substack{ \xi_{2}=t }}}{24} - \frac{h^{2} \left. \frac{\partial^{3}}{\partial \xi_{1}^{3}} u{\left(\xi_{1},t \right)} \right|_{\substack{ \xi_{1}=x }}}{6} + O\left(\tau^{4}\right) + O\left(h^{4}\right)\]

Видим, что порядок аппроксимации \(O(\tau, h^2)\) - неплохо, но к сожалению, данная схема безусловно неустойчива.

Аналогичным образом можно находить порядки аппроксимации любых схем (или граничных условий) - будьте разве что внимательны к разложению точек типа \(u_{m+1}^{n+1}\), где надо раскладывать одновременно по двум и более параметрам, - в sympy нет встроенных методов для вычисления ряда Тейлора по нескольким переменным.

Пример 2. Исследования на аппроксимацию схемы расщепления по направлению

Для последней схемы покажем, как определить порядок аппроксимации.

В этой схеме последовательно решаем задачи сначала на ~, потом на \(n+1\).

\[ \frac{\tilde{u}_{m l}-u_{m l}^n}{\tau}=\mathbf{\Lambda}_1 \tilde{u}_{m l} \]

\[ \frac{u_{m l}^{n+1}-\tilde{u}_{m l}}{\tau}=\boldsymbol{\Lambda}_2 u_{m l}^{n+1} \]

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

Запишем эти уравнения в следующем виде:

\[ \left(\mathbf{E}-\tau \mathbf{\Lambda}_1\right) \tilde{u}_{m l}=u_{m l}^n \]

\[ \left(\mathbf{E}-\tau \boldsymbol{\Lambda}_2\right) u^{n+1}=\tilde{u}_{m l} \]

Подействуем на обе части второго уравнения оператором \(\left(\mathbf{E}-\tau \mathbf{\Lambda}_1\right)\): $\( \begin{aligned} & \quad\left(\mathbf{E}-\tau \mathbf{\Lambda}_1\right)\left(\mathbf{E}-\tau \mathbf{\Lambda}_2\right) u_{m l}^{n+1}=\left(\mathbf{E}-\tau \mathbf{\Lambda}_1\right) \tilde{u}_{m l} \end{aligned} \)$

Так как \(\left(\mathbf{E}-\tau \mathbf{\Lambda}_1\right) \tilde{u}_{m l}=u_{m l}^n\), то

\[ \left(\mathbf{E}-\tau \boldsymbol{\Lambda}_1\right)\left(\mathbf{E}-\tau \boldsymbol{\Lambda}_2\right) u_{m l}^{n+1}=u_{m l}^n \]

Раскрывая скобки,

\[ \frac{u_{m l}^{n+1}-u_{m l}^n}{\tau}=\boldsymbol{\Lambda}_1 u_{m l}^{n+1}+\boldsymbol{\Lambda}_2 u_{m l}^{n+1}-\tau^2 \boldsymbol{\Lambda}_1 \boldsymbol{\Lambda}_2 u_{m l}^{n+1} \]

И здесь уже всё разлагаем в Тейлора. Порядок аппроксимации этой схемы \(O\left(\tau, h_x^2, h_y^2\right)\).

Примечание. Не забываем об аппроксимации граничных условий. Если они заданы в хитром виде (с производными), то их также надо аппроксимировать с нужным порядком.

Анализ устойчивости

Признак фон Неймана

• Самый простой в использовании

• Применяется в эволюционных задачах (когда есть явная зависимость от времени). Т.е. для параболических и гиперболических уравнений.

Введем вектор разностного решения для данного временного слоя

\[ \mathbf{y}^n=\left(y_0^n, y_1^n, \ldots, y_M^n\right)^{\mathrm{T}} \]

Это может быть как одномерный, так и k-мерный массив по пространственным координатам.

Пусть мы сумели (возможно, неединственным образом) представить решение на следующем слое в виде

\[ \mathbf{y}^{n+1}=\mathbf{R}_\tau \mathbf{y}^n+\boldsymbol{\sigma}^n . \]

Оператор \(\mathbf{R}_\tau\) называется оператором послойного перехода.

Теорема. (спектральный признак устойчивости фон Неймана). Необходимым условием устойчивости по начальным данным является выполнение следующего неравенства для всех собственных значений оператора перехода \(\mathbf{R}_\tau\):

\[ |\lambda| \leq 1 \]

Обычно устойчивость исследуется отдельно для каждой фурьегармоники решения $\( y_m^n=\lambda^n e^{i \alpha m} . \)$

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

Явный левый уголок

Рассмотрим схему, называемую явный левый уголок для решения одномерного линейного уравнения переноса:

\[\begin{split} \begin{aligned} & u_t^{\prime}+c u_x^{\prime}=f(t, x) \text { при } c>0: \\ & \frac{y_m^{n+1}-y_m^n}{\tau}+c \frac{y_m^n-y_{m-1}^n}{h}=f_m^n . \end{aligned} \end{split}\]

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

Подставим \(y_m^n=\lambda^n e^{i \alpha m}\) в однородное разностное уравнение, так как из устойчивости по начальным данным следует устойчивость по правой части. В результате получим

\[ \frac{\lambda^{n+1} e^{i \alpha m}-\lambda^n e^{i \alpha m}}{\tau}+c \frac{\lambda^n e^{i \alpha m}-\lambda^n e^{i \alpha(m-1)}}{h}=0 . \]

Нас интересует нетривиальное решение, поэтому сокращая на общий множитель всех членов \(\lambda^n e^{i \alpha}\), получим спектр оператора перехода

\[ \lambda=1-c \tau / h+e^{-i \alpha} c \tau / h . \]

При различных значениях параметра \(\alpha\) получим на комплексной плоскости окружность с центром в точке \(\quad(1-c \tau / h, 0)\) радиуса \(c \tau / h\). При числе Куранта \(c \tau / h \leq 1\) эта окружность лежит внутри единичной (рис. 2.1), тем самым выполнен спектральный признак устойчивости. Таким образом, данная разностная схема является условно устойчивой - она устойчива при (гиперболическом) числе Куранта меньше единицы.

Явная схема

\[ \frac{u_{m l}^{n+1}-u_{m l}^n}{\tau}=\frac{u_{m-1, l}^n-2 u_{m l}^n+u_{m+1, l}^n}{h_x^2}+\frac{u_{m, l-1}^n-2 u_{m l}^n+u_{m, l+1}^n}{h_y^2} \]

или, в операторном виде,

\[ \frac{u_{m l}^{n+1}-u_{m l}^n}{\tau}=\boldsymbol{\Lambda}_1 u_{m l}^n+\boldsymbol{\Lambda}_2 u_{m l}^n \]

Подставим в это уравнение следующую штуку

\[ u_{m l}^n=\lambda^n e^{i \alpha m+i \beta l}, \quad \alpha, \beta \in[0,2 \pi] \]

Сокращая всё, что сокращается, выражаем \(\lambda\):

\[ \lambda\left(\alpha, \beta, \tau, h_x, h_y\right)=1-4 \frac{\tau}{h_x^2} \sin ^2 \frac{\alpha}{2}-4 \frac{\tau}{h_y^2} \sin ^2 \frac{\beta}{2} \]

Дальше надо посмотреть, при каких параметрах \(\tau, h_x, h_y\) все наши \(\lambda\) (тут получилось только одна лямбда) по модулю не превосходят единицу равномерно по \(\alpha\) и \(\beta\). Здесь, очевидно,

\[ \tau \leq \frac{1}{2\left(h_x^{-2}+h_y^{-2}\right)} \]

Неявная схема

\[ \frac{u_{m l}^{n+1}-u_{m l}^n}{\tau}=\frac{u_{m-1, l}^{n+1}-2 u_{m l}^{n+1}+u_{m+1, l}^{n+1}}{h_x^2}+\frac{u_{m, l-1}^{n+1}-2 u_{m l}^{n+1}+u_{m, l+1}^{n+1}}{h_y^2} \]

или, в операторном виде,

\[ \frac{u_{m l}^{n+1}-u_{m l}^n}{\tau}=\boldsymbol{\Lambda}_1 u_{m l}^{n+1}+\boldsymbol{\Lambda}_2 u_{m l}^{n+1} \]

Подставим сюда

\[ u_{m l}^n=\lambda^n e^{i \alpha m+i \beta l}, \quad \alpha, \beta \in[0,2 \pi] \]

И выразим лямбду

\[ \lambda=\frac{1}{1+4 \frac{\tau}{h_x^2} \sin ^2 \frac{\alpha}{2}+4 \frac{\tau}{h_y^2} \sin ^2 \frac{\beta}{2}} \]

Видим, что эта лямбда по модулю меньше единицы равномерно по \(\alpha\) и \(\beta\) при любых \(\tau, h_x, h_y\), т.е. она абсолютно устойчива.

Схема расщепления по направлению

В этой схеме последовательно решаем задачи сначала на ~, потом на \(n+1\).

\[ \frac{\tilde{u}_{m l}-u_{m l}^n}{\tau}=\mathbf{\Lambda}_1 \tilde{u}_{m l} \]

\[ \frac{u_{m l}^{n+1}-\tilde{u}_{m l}}{\tau}=\boldsymbol{\Lambda}_2 u_{m l}^{n+1} \]

По алгоритму применения спектрального признака, надо подставить

\[ u_{m l}^n=\lambda^n e^{i \alpha m+i \beta l}, \quad \alpha, \beta \in[0,2 \pi] \]

где явно записан множитель временной эволюции \(\lambda^n\).

При переходе на слой ~ тоже аля происходит временная эволюция, поэтому

\[ \tilde{u}_{m l}=\lambda_1 u_{m l}^n \]

То же верно при переходе со слоя ~ на \(n+1\):

\[ u_{m l}^{n+1}=\lambda_2 \tilde{u}_{m l}=\left(\lambda_1 \lambda_2\right) u_{m l}^n \]

Откуда \(\lambda=\lambda_1 \lambda_2\).

Подставим формулы выше в разностные схемы и сократим всё, что сокращается. Итого,

\[ \lambda_1=\frac{1-4 \frac{\tau}{h_x^2} \sin ^2 \frac{\alpha}{2}}{1+4 \frac{\tau}{h_y^2} \sin ^2 \frac{\beta}{2}}, \quad \lambda_2=\frac{1-4 \frac{\tau}{h_x^2} \sin ^2 \frac{\beta}{2}}{1+4 \frac{\tau}{h_y^2} \sin ^2 \frac{\alpha}{2}} \]

Окончательный спектр оператора послойного перехода

\[ \lambda\left(\alpha, \beta, \tau, h_x, h_y\right)=\lambda_1 \lambda_2=\frac{\left(1-4 \frac{\tau}{h_x^2} \sin ^2 \frac{\alpha}{2}\right)\left(1-4 \frac{\tau}{h_x^2} \sin ^2 \frac{\beta}{2}\right)}{\left(1+4 \frac{\tau}{h_y^2} \sin ^2 \frac{\alpha}{2}\right)\left(1+4 \frac{\tau}{h_y^2} \sin ^2 \frac{\beta}{2}\right)} \]

Легко видим, что эта штука всегда меньше единицы. Т.е. схема безусловно устойчива.

Примечание. В случае переменных коэффициентов можно использовать принцип замороженных коэффициентов - когда мы считаем коэффициенты константными даже при их непостоянности. Часто приводит к неверным результатам, но существенно упрощает выкладки. В остальном же придётся делать оценку устойчивости равномерно по всей области допустимых параметров \(t, x, y\).