Численные методы решения уравнений математической физики

Методы расщепления

Разбить на страницы
Показывать лекцию целиком

8.1. Понятие о методах расщепления

Рассмотрим дифференциальную задачу для уравнения в частных производных с постоянными коэффициентами:

$$$ {\frac{{\partial}u}{{\partial}t} + {\mathbf{A}}u = 0; x \in \Omega_x , t \in \Omega_t , \left. u\right|_\Gamma = u_\Gamma , u(t_0 ) = u_0 .} $$$

Здесь оператор $${\mathbf{A}} \ge 0$$ — положительный дифференциальный оператор с постоянными коэффициентами. В запись оператора $$\mathbf{A}$$ входят производные по пространственным переменным. Для любого ненулевого элемента выполнено $$(\mathbf{A} \varphi, \varphi) \ge 0$$. $$\Gamma$$ — граница области интегрирования $$\Omega _{x}$$ ; $$\Lambda$$ — разностный оператор, аппроксимирующий $$\mathbf{A}$$. Можно проверить, что разностное уравнение

$$$ {\frac{{u^{{n} + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \frac{{u^{{n} + 1} + u^{n}}}{2} = 0, u^0 = u_0} $$$

аппроксимирует (8.1) со вторым порядком по $$\tau$$ ( схема Кранка - Никольсон ). Заметим, что (8.2) можно трактовать как результат попеременного применения явной и неявной схем первого порядка аппроксимации, записанных на интервалах [tn, t + 1/2], [tn + 1/2, tn + 1]

$$\begin{gather*} \frac{{u^{{n} + 1/2} - u^{n}}}{{{\tau}/2}} + {\mathbf{\Lambda}}u^{n} = 0, \\ {\frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{{{\tau}/2}} + {\mathbf{\Lambda}}u^{{n} + 1} = 0.} \end{gather*}$$

Исключая из уравнений (8.3) значения функции на промежуточном слое по времени (с полуцелым индексом), получим (8.2). Если $$\mathbf{A} = \mathbf{A}(t)$$, то

$$$ {\frac{{u^{{n} + 1} - u^{{n}}}}{\tau} + {\mathbf{\Lambda}}^{n} \frac{{u^{{n} + 1} + u^{{n}}}}{{2}} = 0} $$$

при этом разностный оператор также является положительным:

$$(\Lambda_nu, u) \ge 0,$$

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

$$$ u^{{n} + 1} = ({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}^{n} )^{- 1}({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}^{n})u^{n}, $$$

или

$$u^{{n} + 1} = {\mathbf{T}}^{n} u^{n},$$

где $$$ {\mathbf{T}}^{n} = ({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}^{n} )^{- 1}({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}^{n} ). $$$

Для доказательства устойчивости полученного разностного уравнения умножим скалярно (8.4) на (un + un + 1)/2, получим

$$$ {\frac{(u^{{n} + 1}, u^{{n} + 1}) - (u^{{n}}, u^{{n}})} {2 \tau} + \left({\mathbf{\Lambda}}^{n} \frac{{u^{{n} + 1} + u^{{n}}}}{{2}}, \frac{{u^{{n} + 1} + u^{{n}}}}{{2}}\right) = 0} $$$

Так как в силу положительности разностного оператора $${(\Lambda^{n}u, u) \ge 0}$$, то из (8.5}) следует, что $$$ \left\| {u^{{n} + 1}\right\| \le \left\| {u^{n}}\right\|, $$$ чем и обеспечена устойчивость схемы. Если разностный оператор $$\Lambda$$ (пространственные разности) выбран в виде полусуммы разностных операторов на верхнем и нижнем слоях по времени $$$ \frac{1}{2}({\mathbf{\Lambda}}^{{n} + 1} + {\mathbf{\Lambda}}^{n} ) = {\mathbf{\Lambda}}^{{n} + 1/2}, $$$ то схема имеет второй порядок аппроксимации по $$\tau.$$

8.2. Метод расщепления первого и второго порядка точности по tau

8.2.1. Локально - одномерные схемы

Положим, что дифференциальный оператор $$\mathbf{A}$$ и соответствующий ему разностный оператор $$\Lambda$$ можно представить в виде суммы операторов, каждый из которых включает производные лишь по одной пространственной переменной и разности лишь вдоль одного направления соответственно. Всего пространственных направлений N. Такие дифференциальные и разностные операторы будем называть локально - одномерными. И дифференциальный, и разностный операторы записываются в виде суммы локально - одномерных:

$${\mathbf{A}} = \sum\limits_{{i} = 1}^{N}{{\mathbf{A}}_i }, {\mathbf{\Lambda}} = \sum\limits_i^{N}{{\mathbf{\Lambda}}_i }.$$

Для однородной задачи можно выписать схему расщепления по направлениям:

$$\begin{gather*} \frac{u^{n + 1/N} - u^n}{\tau} + {\mathbf{\Lambda}}_1 u^{n + 1/N} = 0, \\ \frac{u^{n + 2/N} - u^{n + 1/N}}{\tau} + {\mathbf{\Lambda}}_2 u^{n + 2/N} = 0, \\ \ldots \\ \frac{u^{n + 1} - u^{n + (N - 1)/N}}{\tau} + {\mathbf{\Lambda}}_n u^{{n} + 1} = 0. \end{gather*}$$

Получена система разностных уравнений, каждое из которых не аппроксимирует исходное дифференциальное, но может быть легко решено (методом прогонки вдоль соответствующего направления, если разностные операторы содержат лишь первые и вторые разности). Тем не менее, последовательно примененные друг за другом, они дают на следующем слое по времени решение с разумной точностью. Говорят, что имеет место суммарная аппроксимация — результирующий оператор послойного перехода получился аппроксимирующим. Описанный выше способ называется иногда методом дробных шагов, и уже встречался при решении многомерного уравнения теплопроводности.

Для неоднородной задачи один из возможных вариантов схемы расщепления имеет вид

$$\begin{gather*} \frac{u^{n + 1/N} - u^n}{\tau} + {\mathbf{\Lambda}}_1 u^{{n} + 1/{N}} = 0, \\ \frac{u^{n + 2/N} - u^{n + 1/N}}{\tau} + {\mathbf{\Lambda}}_2 u^{n + 2/N} = 0, \\ \ldots \\ \frac{u^{n + 1} - u^{n + (N - 1) /N}}{\tau} + {\mathbf{\Lambda}}_n u^{{n} + 1} = f^{n} . \end{gather*}$$

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

Приведенные выше схемы расщепления по направлениям абсолютно устойчивы.

8.2.2. Схемы Кранка - Никольсон

Рассмотрим обобщение схемы Кранка - Никольсон на случай многомерных уравнений с локально - одномерными операторами. Положим, как и ранее, $${\mathbf{\Lambda}} = \sum\limits_i^{N}{{\mathbf{\Lambda}}_i }$$. Если коэффициенты разностного оператора явно зависят от времени, они берутся на промежуточном временном слое $${\mathbf{\Lambda}} = {\mathbf{\Lambda}}(t^{{n} + 1/2} )$$. Для простоты изложения рассмотрим двумерный случай.

Схему расщепления по направлениям представим в виде

$$\begin{gather*} \frac{{u^{{n} + 1/2} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{{n} + 1/2} + u^{n}}}{2} = 0 \\ \frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{{n} + 1} + u^{{n} + 1/2}}}{2} = 0, \end{gather*}$$

а решение на следующем слое по времени в операторной форме выписывается как $$u^{{n} + 1} = {\mathbf{T}}^{n} u^{n}$$. Для оператора послойного перехода получается следующая формула:

$$$ {\mathbf{T}}^{n} = ({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2 )^{- 1} ({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2)({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1 )^{- 1} ({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1). $$$

При выполнении условия

$$$ \frac{\tau}{2} \left\|{{\mathbf{\Lambda}}_i }\right\| < 1 $$$

схема устойчива, обладает вторым порядком аппроксимации по времени, если операторы $$\Lambda _{1}, \Lambda _{2}$$ коммутативны, и первым — если нет.

8.2.3. Общая формулировка методов расщепления

Заменим локально - одномерные дифференциальные операторы $$\mathbf{A}_i$$ разностными операторами на каждом шаге по времени $$t_n \le t \le t_{{n} + 1}$$. Представим схему расщепления в следующем общем виде:

$$\begin{gather*} \frac{{u^{{n} + 1/{N}} - u^{n}}}{\tau} + {{{\Lambda}}}_{10}u^{n} + {{{\Lambda}}}_{11} u^{{n} + 1/{N}} = 0, \\ \frac{{u^{{n} + 2/{N}} - u^{{n} + 1/{N}}}}{\tau } + {{{\Lambda}}}_{20} u^{n} + {{{\Lambda}}}_{21} u^{{n} + 1/{N}} + {{{\Lambda}}}_{22} u^{{n} + 2/{N}} = 0, \\ \ldots \\ \frac{u^{n + 1} - u^{n + \frac{N - 1}{N}}}{\tau} + {{{\Lambda}}}_{{N}0} u^{n} + {{{\Lambda}}}_{{N}1} u^{{n} + 1/2} + \ldots + {{{\Lambda}}}_{{N}{N}} u^{{n} + 1} = 0. \end{gather*}$$

Условие устойчивости такой схемы расщепления будет

$$\left\|{C_i C_{{i} - 1} \ldots C_1}\right\| \le 1 + c {\tau}, \quad c = const,$$

где $$C_i = {(\mathbf{E} + {\tau}\Lambda_{ii})}^{- 1} (\mathbf{E} + \tau \Lambda_{i, i - 1}), i = 1, 2, \ldots , N$$ Двухслойная схема расщепления с весовыми коэффициентами представлена в виде

$$\begin{gather*} \frac{{u^{{n} + 1/{N}} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \left[{(1 - {\sigma})u^{n} + {\sigma}u^{{n} + 1/{N}}}\right] = 0, \\ \frac{{u^{{n} + 2/{N}} - u^{{n} + 1/{N}}}}{\tau } + {\mathbf{\Lambda}}_2 \left[{(1 - \sigma )u^{{n} + 1/{N}} + \sigma u^{{n} + 2/{N}}}\right] = 0, \\ \ldots \\ \frac{{u^{{n} + 1} - u^{{n} + \frac{{{N} - 1}}{{N}}}}}{\tau} + {\mathbf{\Lambda}}_n \left[ {(1 - {\sigma})u^{{n} + \frac{{{N} - 1}}{{N}}} + {\sigma}u^{{n} + 1}}\right] = 0. \end{gather*}$$

Если в этой схеме расщепления положить веса верхнего и нижнего слоев по времени равными, $$\gamma = 0, 5$$, то в случае коммутирующих операторов $$\Lambda _{i}$$ (каждый такой разностный оператор аппроксимирует со вторым порядком соответствующий локально - одномерный дифференциальный оператор) схема будет иметь второй порядок аппроксимации и по времени. Если при этом каждый оператор $${\mathbf{\Lambda}}_i \ge 0$$, то схема будет абсолютно устойчивой.

8.2.4. Схемы расщепления для уравнения теплопроводности

Рассматриваем нестационарное уравнение теплопроводности

$$$ \frac{{\partial}u}{{\partial}t} - \Delta u = 0, t \in \Omega_t , \left\{{x, y, z}\right\} \in \Omega . $$$

Здесь оператор Лапласа определен как $$\Delta = {\partial}^2 / {\partial}x^2 + {\partial}^2 / {\partial}y^2 + {\partial}^2 / {\partial}z^2$$. Его также можно записать в виде суммы трех локально - одномерных операторов $$\mathbf{A} = \mathbf{A}_x + \mathbf{A}_y + \mathbf{A}_{z}$$. Соответствующие разностные операторы будут $${\mathbf{\Lambda}} = {\mathbf{\Lambda}}_{xx} + {\mathbf{\Lambda}}_{yy} + {\mathbf{\Lambda}}_{zz}$$, где $$$ \Lambda_{xx} = \frac{u_{m - 1, jk} - 2u_{mjk} + u_{m + 1, jk}}{h_x^2} $$$ аналогично определяются операторы вычисления второй разностной производной и по остальным направлениям $$\Lambda _{yy}, \Lambda _{zz}$$ Локально - одномерная схема для уравнения теплопроводности будет

$$\begin{gather*} \frac{{u^{{n} + 1/3} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_{xx}u^{{n} + 1/3} = 0, \\ \frac{{u^{n + 2/3} - u^{{n} + 1/3}}}{\tau} + {\mathbf{\Lambda}}_{yy}u^{{n} + 2/3} = 0, \\ \frac{{u^{{n} + 1} - u^{{n} + 2/3}}}{\tau} + {\mathbf{\Lambda}}_{zz}u^{{n} + 1} = 0. \end{gather*}$$

Для повышения порядка аппроксимации можно использовать схему с весами

$$\begin{gather*} \frac{{u^{{n} + 1/3} - u^{n}}}{\tau } + {\mathbf{\Lambda}}_{xx} \left[{(1 - {\sigma})u^n + {\sigma}u^{n + 1/3}}\right] = 0, \\ \frac{{u^{{n} + 2/3} - u^{{n} + 1/3}}}{\tau} + {\mathbf{\Lambda}}_{yy} \left[{(1 - {\sigma})u^{{n} + 1/3} + {\sigma}u^{{n} + 2/3}}\right] = 0, \\ \frac{{u^{{n} + 1} - u^{{n} + 2/3}}}{\tau} + {\mathbf{\Lambda}}_{zz}[(1 - {\sigma})u^{{n} + 2/3} + {\sigma}u^{{n} + 1} ] = 0. \end{gather*}$$

8.3. Методы двуциклического покомпонентного расщепления

Для этих методов отсутствует требование коммутативности операторов $$\Lambda _{i}.$$

Будем рассматривать численное решение (8.1) не на одном шаге по времени, отрезке [tn, tn + 1], а на двух последовательных шагах [tn - 1, tn + 1]. Пусть теперь разностные локально - одномерные операторы зависят явно от времени, тогда они определены в середине отрезка $${\mathbf{\Lambda}}_i = {\mathbf{\Lambda}}_i (t^{n} )$$. Запишем схему расщепления:

$$\begin{gather*} \frac{{u^{{n} - 1/2} - u^{{n} - 1}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{{n} - 1/2} + u^{{n} - 1}}}{2} = 0, \\ \frac{{u^{n} - u^{{n} - 1/2}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n} + u^{{n} - 1/2}}}{2} = 0, \\ \frac{{u^{{n} + 1/2} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{{n} + 1/2} - u^{n}}}{2} = 0, \\ {\frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{2} = 0.} \end{gather*}$$

В операторной форме этот метод записывается как $$u^{n + 1} = {\mathbf{T}}^{n}u^{n - 1}$$, где введено обозначение

$$\begin{gather*} {\mathbf{T}}^{n} = \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right) \left({{\mathbf{E}} + \frac{\tau}{2} {\mathbf{\Lambda}}_2}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right) \times \\ \times \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)^{- 1} \left({{\mathbf{E}}- \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right) \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right) = \mathbf{E}{- } 2 {\tau}{\Lambda}+ \frac{{(2 \tau)}^2}{2}{(\Lambda)}^2 + \ldots \end{gather*}$$

Если локально - одномерные операторы положительны $${\mathbf{A}}_i (t) > 0$$, то при достаточной гладкости решения и элементов матриц $${\mathbf{A}}_i (t)$$ схема (8.6) абсолютно устойчива и аппроксимирует (8.1) со вторым порядком.

Для неоднородного дифференциального уравнения $$$ \frac{{\partial u}}{{\partial}t} + {\mathbf{A}}u = f $$$ разностная аппроксимация метода расщепления может быть представлена в виде

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - 1/2} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - 1}, \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)(u^{n} - {\tau}f^n) = \left({{\mathbf{E}} - \frac{\tau }{2}{\mathbf{\Lambda}}_2}\right)u^{n - 1/2}, \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)u^{n + 1/2} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)(u^{n} + {\tau}f^{n} ), \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + 1} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + 1/2}, \end{gather*}$$

где fn = f(tn).

В операторной форме записи решение неоднородной задачи имеет вид: $$u^{n + 1} = {\mathbf{T}}^{n} u^{n - 1} + 2 {\tau}{\mathbf{T}}_1^{n}{\mathbf{T}}_2^{n} f^{n}$$, где введены обозначения $$$ {\mathbf{T}}^{n} = {\mathbf{T}}_1^{n}{\mathbf{T}}_2^{n}{\mathbf{T}}_2^{n}{\mathbf{T}}_1^{n}$$ $${\mathbf{T}}_i^{n} = \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_i^{n}}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau }{2}{\mathbf{\Lambda}}_i^{n}}\right). $$$ Представим разностную аппроксимацию неоднородного дифференциального уравнения с помощью последовательного применения операторов $$\Lambda _{1}, \dots , \Lambda _{n}$$, $$\Lambda _{n}, \dots , \Lambda _{1}:$$

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - \frac{N - 1}{N}} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - 1}, \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)(u^{n} - {\tau}f^{n} ) = \left({{\mathbf{E}} - \frac{\tau }{2}{\mathbf{\Lambda}}_n }\right)u^{n - 1/N}, \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)u^{n + 1/N} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)(u^{n} + {\tau}f^{n} ), \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + 1} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + \frac{{N - 1}}{N}} . \end{gather*}$$

Рассмотрим примеры использования метода двуциклического расщепления для некоторых задач математической физики.

Пример 1. Трехмерное нестационарное уравнение диффузии, область интегрирования — параллелепипед. Полагаем, что в вертикальном направлении (ось 0z ) коэффициент диффузии в вертикальной плоскости $$\gamma$$ зависит от координаты, что характерно для задач геофизики, $$\mu$$ — коэффициент диффузии в горизонтальной плоскости. Задача может быть представлена в виде

$$$ \frac{{\partial}u}{{\partial}t} = \frac{\partial }{{{\partial}z}} \gamma \frac{{\partial}u}{{{\partial}z}} + \mu \Delta u + f. $$$

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

$$$ \frac{{{\partial}u_1}}{{\partial}t} = \frac{{\partial}}{{{\partial}z}} \gamma \frac{{{\partial}u_1}}{{{\partial}z}} + f, $$$

она описывает диффузию в вертикальной плоскости. Вторую и третью задачи запишем

$$$ \frac{{{\partial}u_2}}{{\partial}t} = \mu \frac{{{\partial}^2 u_2}} {{{\partial}x^2}}, \frac{{{\partial}u_3 }}{{\partial}t} = \mu \frac{{{\partial}^2 u_3 }} {{{\partial}y^2}}. $$$

Теперь рассмотрим разностную аппроксимацию исходного дифференциального уравнения

$$\begin{gather*} \frac{{\partial}u}{{\partial}t} + ({\mathbf{\Lambda}}_1 + {\mathbf{\Lambda}}_2 + {\mathbf{\Lambda}}_3 )u = f, \mbox{ где} \\ {\mathbf{\Lambda}}_1 u = -\frac{\mu}{h_x^2}(u_{m, j, k + 1} - 2u_{mjk}+ u_{m, j, k - 1} ),\\ {\mathbf{\Lambda}}_2 u = -\frac{\mu}{h_y^2}(u_{m, j - 1, k} - 2u_{mjk} + u_{m, j + 1, k}),\\ {\mathbf{\Lambda}}_3 u = \frac{1}{{h_{z}}} \left[{- \frac{{\gamma_{m + 1/2, jk}}}{{h_{z}}}(u_{m + 1, jk} - u_{mjk} ) + \frac{{\gamma_{m - 1/2}}}{{h_{z}}} (u_{mjk} - u_{m - 1, jk} )}\right]. \end{gather*}$$

Разностная схема двуциклического покомпонентного расщепления приобретает вид

$$\begin{gather*} \frac{{u^{n + 1/6} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1/6} + u^{n}}}{2} = 0, \\ \frac{{u^{n + 2/6} - u^{n + 1/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 2/6} + u^{n + 1/6}}}{2} = 0, \\ \frac{{u^{n + 3/6} - u^{n + 2/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 3/6} + u^{n + 2/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 4/6} - u^{n + 3/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 4/6} + u^{n + 3/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 5/6} - u^{n + 4/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 5/6} + u^{n + 4/6}}}{2} = 0, \\ \frac{{u^{n + 1} - u^{n + 5/6}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1} + u^{n + 5/6}}}{2} = 0. \end{gather*}$$

Пример 2. Сопряженное нестационарное уравнение переноса и диффузии

$$$ \frac{{\partial}u}{{\partial}t} + ({\mathbf{A}}_1 + {\mathbf{A}}_2 + {\mathbf{A}}_3 )u = f, $$$

где операторы определены как

$$$ {\mathbf{A}}_1 = \frac{{{\partial}(v_{1} u)}}{{\partial}x} - \mu \frac{{{\partial}^2 u}}{{{\partial}x^2}}, {\mathbf{A}}_2 = \frac{{{\partial}(v_2u)}}{{\partial}y} - \mu \frac{{{\partial}^2 u}}{{{\partial}y^2}}, {\mathbf{A}}_3 = \frac{{{\partial}(v_3u)}}{{{\partial}z}} - \frac{{\partial}}{{{\partial}z}} \gamma \frac{{\partial}u}{{{\partial}z}} + {\sigma}u. $$$

Здесь v1, v2, v3 — компоненты вектора скорости, u — концентрация субстанции, $$\sigma$$ — коэффициент поглощения субстанции внешней средой, $$\sigma > 0$$. Соответствующую схему расщепления представим в виде

$$\begin{gather*} \frac{{u^{n + 1/6} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1/6} + u^{n}}}{2} = 0, \\ \frac{{u^{n + 2/6} - u^{n + 1/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 2/6} + u^{n + 1/6}}}{2} = 0, \\ \frac{{u^{n + 3/6} - u^{n + 2/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 3/6} + u^{n + 2/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 4/6} - u^{n + 3/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 4/6} + u^{n + 3/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 5/6} - u^{n + 4/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 5/6} + u^{n + 4/6}}}{2} = 0, \\ \frac{{u^{n + 1} - u^{n + 5/6}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1} + u^{n + 5/6}}}{2} = 0, \end{gather*}$$

где разностные операторы аппроксимируют соответствующие локально - одномерные дифференциальные операторы. Так для рассматриваемой задачи

$$\begin{gather*} {\mathbf{\Lambda}}_1 u = - \frac{\mu }{{h_x^2}}(u_{m + 1, jk} - 2u_{mjk} + u_{m - 1j, k} ) + \frac{{(v_1 u)_{m + 1, jk} - (v_1 u)_{m - 1, jk}}}{{2h_x }}, \\ {\mathbf{\Lambda}}_2 u = - \frac{\mu }{{h_y^2}}(u_{m, j + 1, k} - 2u_{mjk} + u_{m, j - 1, k} ) + \frac{{(v_2 u)_{m, j + 1, k} - (v_2 u)_{m, j - 1, k}}}{{2h_y}}, \\ {\mathbf{\Lambda}}_3 u = \frac{1}{{h_{z}}} \left[{- \frac{{\gamma_{m + 1/2, jk}}}{{h_{z}}}(u_{m + 1, jk} - u_{mjk} ) + \frac{{\gamma_{m - 1/2}}}{h}(u_{mjk} - u_{m - 1, jk} )}\right] + \\ + \frac{{(v_3 u)_{m, j, k + 1} - (v_3 u)_{m, j, k - 1}}}{{2h_{z}}} + {\sigma}u_{mjk} . \end{gather*}$$

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

$$$ \frac{{\partial}u}{{\partial}t} + \frac{{{\partial}(v_1 , u)}} {{\partial}x} + \frac{{{\partial}(v_2 , u)}}{{\partial}y} - \mu \Delta u - {\sigma}u = f, $$$

в котором компоненты скорости движения среды v1, v2 удовлетворяют уравнению неразрывности:

$$$ \frac{{{\partial}v_1}}{{\partial}x} + \frac{{{\partial}v_2}} {{\partial}y} = 0. $$$

Первый этап расщепления задачи по физическим процессам связан с переносом, на нем решается разностный аналог уравнения

$$$ \frac{{{\partial}u_1}}{{\partial}t} + \frac{{{\partial}(v_1 , u)}} {{\partial}x} + \frac{{{\partial}(v_2 , u)}}{{\partial}y} = 0. $$$

Второй этап расщепления описывает процессы диффузии и поглощения субстанций

$$$ \frac{{{\partial}u_2}}{{\partial}t} - \mu \Delta u + {\sigma}u_2 = f. $$$

Пример 3. Расщепление по физическим процессам системы уравнений газовой динамики (метод крупных частиц).

Система

$$\begin{gather*} \frac{{\partial}\rho}{{\partial}t} + div ({\rho}{\mathbf{v}}) = 0, \\ \frac{{\partial}({\rho}u_1)}{{\partial}t} + div ({\rho}{\mathbf{v}}u_1 ) + \frac{{\partial P}}{{\partial}x} = 0, \\ \frac{{{\partial}({\rho}u_2 )}}{{\partial}t} + div ({\rho}{\mathbf{v}}u_2 ) + \frac{{\partial P}}{{\partial}y} = 0, \\ \frac{{{\partial}({\rho}e)}}{{\partial}t} + div ({\rho}e{\mathbf{v}}) + \div (P{\mathbf{v}}) = 0, \\ P = P({\rho}, \varepsilon ), e = \varepsilon + \frac{{u_1^2 + u_{2}^{2}}}{2}, \end{gather*}$$

u1 , u2 — компоненты вектора скорости $$\mathbf{v}$$, P — давление газа, $$\rho$$ — плотность, $$\varepsilon$$ — внутренняя энергия.

Первый (Эйлеров) этап. Измеряются лишь величины, относящиеся к ячейке в целом, а жидкость внутри каждой ячейки сетки считается моментально замороженной. Расчет производится по формулам, аппроксимирующим уравнения

$$\begin{gather*} \frac{{{\partial}{\rho}}}{{\partial}t} = 0, \\ {\rho}\frac{{{\partial}u_1}}{{\partial}t} + \frac{{\partial}p}{{\partial}x} = 0, {\rho}\frac{{{\partial}u_2}}{{\partial}t} + \frac{{\partial}p}{{\partial}y} = 0, \\ {\rho}\frac{{\partial e}}{{\partial}t} + div (P{v}) = 0. \end{gather*}$$

На втором (Лагранжевом) этапе происходит движение газа массы через границы эйлеровых ячеек и перераспределение массы, импульса, энергии по пространству; определяются поля параметров течения газа. Аппроксимируется система уравнений

$$\begin{gather*} \frac{{{\partial}\rho }}{{\partial}t} + div ({\rho}{\mathbf{v}}) = 0, \\ \frac{{{\partial}({\rho}u_1 )}}{{\partial}t} + div ({\rho}u_1 {\mathbf{v}}) = 0, \frac{{{\partial}({\rho}u_2 )}}{{\partial}t} + div ({\rho}u_2 {\mathbf{v}}) = 0, \\ \frac{{{\partial}({\rho}e)}}{{\partial}t} + div ({\rho}e{\mathbf{v}}) = 0. \end{gather*}$$

8.4. Методы расщепления с факторизацией оператора

8.4.1. Факторизованная схема расщепления

Пусть для решения дифференциальной задачи

$$$ \mathbf{B} \frac{{\partial}u}{{\partial}t} + \mathbf{A}u = f, u(t_0 ) = u_0 $$$

используется разностная схема $${\mathbf{B}}u^{n + 1} = F^{n}$$, где $$F^{n} = ({\mathbf{B}} - {\tau}{\mathbf{A}})u^{n} + {\tau}f^{n}$$, n = 0, 1, ... Пусть для вычисления Fn затрачивается O(N) действий, число арифметических операций пропорционально числу узлов сетки N. Такие разностные операторы называются экономичными.

Пусть $$\mathbf{B}_i$$ (i = 1, 2, ..., N) - экономичные разностные операторы, такие, что $${\mathbf{B}}_i v = F$$. Назовем схему разностной схемой с факторизованным оператором $$\mathbf{B}$$, если возможно его представление в виде

$$\mathbf{B} = \mathbf{B}_1 \mathbf{B}_2 \ldots \mathbf{B}_n.$$

Эта схема будет также экономичной, так как для решения разностного уравнения по - прежнему потребуется O(N) действий. В самом деле, решение уравнения

$$\mathbf{B}_1 \mathbf{B}_2 \ldots \mathbf{B}_n u^{n + 1} = F^{n} $$

может быть найдено в результате последовательного решения p уравнений

$$\begin{gather*} {\mathbf{B}}_1 u_1 = F^{n}, \\ {\mathbf{B}}_2 u_2 = u_1 , \\ {\mathbf{B}}_3 u_3 = u_2, \\ \ldots \\ {\mathbf{B}}_i u_i = u_{i - 1}, \end{gather*} $$

здесь i = 2, 3, ..., N. Тогда un + 1 = un. В записи задачи введены обозначения $$$ u_1 = u^{n + 1/N}, \ldots , u_i = u^{n + i/N}, \ldots , u_{N - 1} = u^{n + \frac{{N - 1}}{N}} $$$ — промежуточные значения.

Схемы с факторизованным оператором иногда называются также факторизованными схемами. Устойчивая схема с факторизованным оператором $$\mathbf{B}$$, которая представляет собой произведение конечного числа операторов $$\mathbf{B}_1, \ldots , \mathbf{B}_n$$, является экономичной схемой.

Пример. Метод переменных направлений (продольно - поперечная схема). Приведем запись схемы для решения линейного двумерного уравнения теплопроводности. Расчетные формулы есть

$$\begin{gather*} \frac{{u^{n + 1/2} - u^{n}}}{{{{1 \over 2}} \tau }} - ({\mathbf{\Lambda}}_1 u^{n + 1/2} + {\mathbf{\Lambda}}_2 u^{n} ) = f^{n}, \\ \frac{{u^{n + 1} - u^{n + 1/2}}}{{{{1 \over 2}} {\tau}}} - ({\mathbf{\Lambda}}_1 u^{n + 1/2} + {\mathbf{\Lambda}}_2 u^{n + 1} ) = f^{n} . \end{gather*}$$

Тогда, исключая un + 1/2, получим в операторной форме записи

$$$ \left(\mathbf{E} - \frac{\tau}{2}\mathbf{\Lambda}_1\right) \left(\mathbf{E} - \frac{\tau}{2}\mathbf{\Lambda}_2\right) u^{n + 1} = \left(\mathbf{E} + \frac{\tau}{2}\mathbf{\Lambda}_1\right) \left(\mathbf{E} + \frac{\tau}{2} \mathbf{\Lambda}_2\right) u^{n}, $$$

или $$\mathbf{B}_1\mathbf{B}_2u^{n + 1} = F^{n}$$, где $$$ {\mathbf{B}}_1 = {\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1 $$$, $$$ {\mathbf{B}}_2 = {\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2 $$$, $$$ F^{n} = \left(\mathbf{E} + \frac{\tau}{2}\mathbf{\Lambda}_1\right) \left(\mathbf{E} + \frac{\tau}{2}\mathbf{\Lambda}_2\right)u^{n} $$$. Разностная схема может быть представлена в виде факторизованный схемы расщепления:

$${\mathbf{B}}_1 u_1 = F^{n}, {\mathbf{B}}_2 u^{n + 1} = u_1.$$

8.4.2. Неявная схема расщепления с приближенной факторизацией

Рассмотрим неявную разностную схему

$$$ {\frac{{u^{n + 1} - u^{n}}}{\tau} + {\Lambda}u^{n + 1} = 0, n = 0, 1, \ldots , {\mathbf{\Lambda}} = \sum\limits_{i = 1}^{N}{{\mathbf{\Lambda}}_i }, \mathbf{\Lambda}_i > 0.} $$$

Представим разностную схему (8.7) в виде

$${({\mathbf{E}} + {\tau}{\mathbf{\Lambda}})u^{n + 1} = u^{n} .}$$

Факторизуем разностную схему (8.8) приближенно с точностью до членов порядка $$O(\tau ^{2}).$$ Для этого заменим в (8.8) оператор $${\mathbf{E}} + {\tau}{\mathbf{\Lambda}}$$ на факторизованный

$$({\mathbf{E}} + \tau {\mathbf{\Lambda}}_1 )({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_2 ) \ldots ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_n ) = {\mathbf{E}} + {\tau}{\mathbf{\Lambda}} + {\tau}^2 {\mathbf{R}},$$

где введено обозначение

$${\mathbf{R}} = \sum\limits_{i < j}{{\mathbf{\Lambda}}_i {\mathbf{\Lambda}}_j + {\tau}\sum\limits_{i < j < k}{{\mathbf{\Lambda}}_i }}{\mathbf{\Lambda}}_j {\mathbf{\Lambda}}_k + \ldots + {\tau}^{n - 2}{\mathbf{\Lambda}}_1 \ldots {\mathbf{\Lambda}}_n.$$

Врезультате приходим к неявной схеме с приближенной факторизацией

$${\mathbf{B}}u^{n + 1} = u^{n}, \mathbf{B} = {\mathop \Pi\limits_{i = 1}^{n}{\mathbf{B}_i}}, {\mathbf{B}}_i = {\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_i,$$

или

$$\begin{gather*} ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_1 )u^{n + 1/N} = u^{n}, \\ ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_2 )u^{n + 2/N} = u^{n + 1/N}, \\ \ldots \\ ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_n )u^{n + 1} = u^{n + \frac{{N - 1}}{N}} . \end{gather*}$$

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

8.4.3. Метод "предиктор - корректор"

Основная идея методов типа "предиктор - корректор" заключается в следующем. На каждом отрезке [tn, tn + 1] задача решается в два приема: сначала по схеме первого порядка аппроксимации и со значительным запасом устойчивости находится решение в момент времени $$t^{n + 1/2} = t^{n} + \tau /2$$ — предиктор. После этого на втором этапе расписывается исходное уравнение по схеме более высокого порядка аппроксимации (чаще всего, второго) — корректор. Основная идея семейств таких методов близка к идее построения методов типа Рунге - Кутты для обыкновенных дифференциальных уравнений.

Представим эту схему как следующую схему расщепления:

$$\begin{gather*} \frac{{u^{n + 1/4} - u^{n}}}{{{\tau}/2}} + {\mathbf{\Lambda}}_1u^{n + 1/4} = 0, \\ \frac{{u^{n + 1/2} - u^{n + 1/4}}}{{{\tau}/2}} + {\mathbf{\Lambda}}_2 u^{n + 1/2} = 0, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = 0. \end{gather*} $$

Если в этой схеме расщепления исключить un + 1/4, то получим последовательность расчетных формул

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right) \left({{\mathbf{E}} + \frac{\tau}{2} {\mathbf{\Lambda}}_2}\right)u^{n + 1/2} = {\varphi}^{n}, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = 0, \end{gather*}$$

далее, исключив, $${\mathbf{\Lambda}}u^{n + 1/2}$$, получим

$$$ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)^{- 1} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)^{- 1}u^{n} = 0. $$$

Если для разностных операторов выполнены условия $$$ \frac{\tau}{2} \left\| {{\mathbf{\Lambda}}_i }\right\| < 1, $$$ $${\mathbf{\Lambda}}_1 \ge 0$$, $${\mathbf{\Lambda}}_2 \ge 0$$, а коэффициенты разностной схемы явно не зависят от времени, то при достаточной гладкости решения дифференциальной задачи разностная схема абсолютно устойчива и аппроксимирует исходную задачу со вторым порядком.

Далее рассмотрим случай, когда оператор $$\mathbf{\Lambda}$$ представляется в виде суммы операторов $$\Lambda _{i}:$$$${\mathbf{\Lambda}} = \sum\limits_i {{\mathbf{\Lambda}}_i } .$$ Пусть все эти разностные операторы положительны. Метод "предиктор - корректор" можно записать в виде последовательности расчетных формул

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2} {\mathbf{\Lambda}}_1}\right)u^{n + 1/2N} = u^{n} + \frac{\tau}{2}f^{n + 1/2}, \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)u^{n + 2/2N} = u^{n + 1/2N}, \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)u^{n + 1/2} = u^{n + \frac{N}{{2N}}}, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = f^{n + 1/2} . \end{gather*}$$

Эта последовательность после исключения промежуточных этапов сводится к одному разностному уравнению

$$$ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \prod\limits_{i = N}^1 {\left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_i } \right)} ^{- 1} (u^{n} + \frac{\tau}{2}f^{n + 1/2} ) = f^{n + 1/2} . $$$

Приведем пример построения такой схемы. Для нестационарного трехмерного уравнения теплопроводности

$$$ \frac{{\partial}u}{{\partial}t} - a^2 \sum\limits_{i = 1}^3 {\frac{{{\partial}^2 u}}{{{\partial}x_i^2}} = 0} $$$

получим следующую разностную схему типа "предиктор - корректор"

$$\begin{gather*} \frac{{u^{n + 1/6} - u^{n}}}{{{\tau}/2}} + {\mathbf{\Lambda}}_1 u^{n + 1/6} = 0, \\ \frac{u^{n + 1/3} - u^{n + 1/6}}{{\tau}/2} + {\mathbf{\Lambda}}_2 u^{n + 1/3} = 0, \\ \frac{u^{n + 1/2} - u^{n + 1/3}}{{\tau}/2} + {\mathbf{\Lambda}}_3 u^{n + 1/2} = 0, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = 0. \end{gather*}$$

Далее, исключая промежуточные слои, получим схему, записанную в каноническом виде

$$\begin{gather*} \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \frac{{u^{n + 1} + u^{n}}}{2} + \frac{{{\tau}^2}}{4}({\mathbf{\Lambda}}_1 {\mathbf{\Lambda}}_2 + {\mathbf{\Lambda}}_1 {\mathbf{\Lambda}}_3 + {\mathbf{\Lambda}}_2 {\mathbf{\Lambda}}_3 ) \frac{{u^{n + 1} - u^{n}}}{\tau} + \\ + \frac{{{\tau}^3 }}{8}{\mathbf{\Lambda}}_1 {\mathbf{\Lambda}}_2 {\mathbf{\Lambda}}_3 \frac{{u^{n + 1} - u^{n}}}{\tau} = 0. \end{gather*}$$

Схема абсолютно устойчива (для коммутирующих операторов), имеет второй порядок аппроксимации по $$\tau$$ и hi. Конечно, при практическом решении задач на компьютере используется именно последовательность разностных операторов. Канонический вид схемы удобен для ее теоретического исследования.

Страницы:

8.1. Понятие о методах расщепления

Рассмотрим дифференциальную задачу для уравнения в частных производных с постоянными коэффициентами:

$$$ {\frac{{\partial}u}{{\partial}t} + {\mathbf{A}}u = 0; x \in \Omega_x , t \in \Omega_t , \left. u\right|_\Gamma = u_\Gamma , u(t_0 ) = u_0 .} $$$

Здесь оператор $${\mathbf{A}} \ge 0$$ — положительный дифференциальный оператор с постоянными коэффициентами. В запись оператора $$\mathbf{A}$$ входят производные по пространственным переменным. Для любого ненулевого элемента выполнено $$(\mathbf{A} \varphi, \varphi) \ge 0$$. $$\Gamma$$ — граница области интегрирования $$\Omega _{x}$$ ; $$\Lambda$$ — разностный оператор, аппроксимирующий $$\mathbf{A}$$. Можно проверить, что разностное уравнение

$$$ {\frac{{u^{{n} + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \frac{{u^{{n} + 1} + u^{n}}}{2} = 0, u^0 = u_0} $$$

аппроксимирует (8.1) со вторым порядком по $$\tau$$ ( схема Кранка - Никольсон ). Заметим, что (8.2) можно трактовать как результат попеременного применения явной и неявной схем первого порядка аппроксимации, записанных на интервалах [tn, t + 1/2], [tn + 1/2, tn + 1]

$$\begin{gather*} \frac{{u^{{n} + 1/2} - u^{n}}}{{{\tau}/2}} + {\mathbf{\Lambda}}u^{n} = 0, \\ {\frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{{{\tau}/2}} + {\mathbf{\Lambda}}u^{{n} + 1} = 0.} \end{gather*}$$

Исключая из уравнений (8.3) значения функции на промежуточном слое по времени (с полуцелым индексом), получим (8.2). Если $$\mathbf{A} = \mathbf{A}(t)$$, то

$$$ {\frac{{u^{{n} + 1} - u^{{n}}}}{\tau} + {\mathbf{\Lambda}}^{n} \frac{{u^{{n} + 1} + u^{{n}}}}{{2}} = 0} $$$

при этом разностный оператор также является положительным:

$$(\Lambda_nu, u) \ge 0,$$

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

$$$ u^{{n} + 1} = ({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}^{n} )^{- 1}({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}^{n})u^{n}, $$$

или

$$u^{{n} + 1} = {\mathbf{T}}^{n} u^{n},$$

где $$$ {\mathbf{T}}^{n} = ({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}^{n} )^{- 1}({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}^{n} ). $$$

Для доказательства устойчивости полученного разностного уравнения умножим скалярно (8.4) на (un + un + 1)/2, получим

$$$ {\frac{(u^{{n} + 1}, u^{{n} + 1}) - (u^{{n}}, u^{{n}})} {2 \tau} + \left({\mathbf{\Lambda}}^{n} \frac{{u^{{n} + 1} + u^{{n}}}}{{2}}, \frac{{u^{{n} + 1} + u^{{n}}}}{{2}}\right) = 0} $$$

Так как в силу положительности разностного оператора $${(\Lambda^{n}u, u) \ge 0}$$, то из (8.5}) следует, что $$$ \left\| {u^{{n} + 1}\right\| \le \left\| {u^{n}}\right\|, $$$ чем и обеспечена устойчивость схемы. Если разностный оператор $$\Lambda$$ (пространственные разности) выбран в виде полусуммы разностных операторов на верхнем и нижнем слоях по времени $$$ \frac{1}{2}({\mathbf{\Lambda}}^{{n} + 1} + {\mathbf{\Lambda}}^{n} ) = {\mathbf{\Lambda}}^{{n} + 1/2}, $$$ то схема имеет второй порядок аппроксимации по $$\tau.$$

8.2. Метод расщепления первого и второго порядка точности по tau

8.2.1. Локально - одномерные схемы

Положим, что дифференциальный оператор $$\mathbf{A}$$ и соответствующий ему разностный оператор $$\Lambda$$ можно представить в виде суммы операторов, каждый из которых включает производные лишь по одной пространственной переменной и разности лишь вдоль одного направления соответственно. Всего пространственных направлений N. Такие дифференциальные и разностные операторы будем называть локально - одномерными. И дифференциальный, и разностный операторы записываются в виде суммы локально - одномерных:

$${\mathbf{A}} = \sum\limits_{{i} = 1}^{N}{{\mathbf{A}}_i }, {\mathbf{\Lambda}} = \sum\limits_i^{N}{{\mathbf{\Lambda}}_i }.$$

Для однородной задачи можно выписать схему расщепления по направлениям:

$$\begin{gather*} \frac{u^{n + 1/N} - u^n}{\tau} + {\mathbf{\Lambda}}_1 u^{n + 1/N} = 0, \\ \frac{u^{n + 2/N} - u^{n + 1/N}}{\tau} + {\mathbf{\Lambda}}_2 u^{n + 2/N} = 0, \\ \ldots \\ \frac{u^{n + 1} - u^{n + (N - 1)/N}}{\tau} + {\mathbf{\Lambda}}_n u^{{n} + 1} = 0. \end{gather*}$$

Получена система разностных уравнений, каждое из которых не аппроксимирует исходное дифференциальное, но может быть легко решено (методом прогонки вдоль соответствующего направления, если разностные операторы содержат лишь первые и вторые разности). Тем не менее, последовательно примененные друг за другом, они дают на следующем слое по времени решение с разумной точностью. Говорят, что имеет место суммарная аппроксимация — результирующий оператор послойного перехода получился аппроксимирующим. Описанный выше способ называется иногда методом дробных шагов, и уже встречался при решении многомерного уравнения теплопроводности.

Для неоднородной задачи один из возможных вариантов схемы расщепления имеет вид

$$\begin{gather*} \frac{u^{n + 1/N} - u^n}{\tau} + {\mathbf{\Lambda}}_1 u^{{n} + 1/{N}} = 0, \\ \frac{u^{n + 2/N} - u^{n + 1/N}}{\tau} + {\mathbf{\Lambda}}_2 u^{n + 2/N} = 0, \\ \ldots \\ \frac{u^{n + 1} - u^{n + (N - 1) /N}}{\tau} + {\mathbf{\Lambda}}_n u^{{n} + 1} = f^{n} . \end{gather*}$$

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

Приведенные выше схемы расщепления по направлениям абсолютно устойчивы.

8.2.2. Схемы Кранка - Никольсон

Рассмотрим обобщение схемы Кранка - Никольсон на случай многомерных уравнений с локально - одномерными операторами. Положим, как и ранее, $${\mathbf{\Lambda}} = \sum\limits_i^{N}{{\mathbf{\Lambda}}_i }$$. Если коэффициенты разностного оператора явно зависят от времени, они берутся на промежуточном временном слое $${\mathbf{\Lambda}} = {\mathbf{\Lambda}}(t^{{n} + 1/2} )$$. Для простоты изложения рассмотрим двумерный случай.

Схему расщепления по направлениям представим в виде

$$\begin{gather*} \frac{{u^{{n} + 1/2} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{{n} + 1/2} + u^{n}}}{2} = 0 \\ \frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{{n} + 1} + u^{{n} + 1/2}}}{2} = 0, \end{gather*}$$

а решение на следующем слое по времени в операторной форме выписывается как $$u^{{n} + 1} = {\mathbf{T}}^{n} u^{n}$$. Для оператора послойного перехода получается следующая формула:

$$$ {\mathbf{T}}^{n} = ({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2 )^{- 1} ({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2)({\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1 )^{- 1} ({\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1). $$$

При выполнении условия

$$$ \frac{\tau}{2} \left\|{{\mathbf{\Lambda}}_i }\right\| < 1 $$$

схема устойчива, обладает вторым порядком аппроксимации по времени, если операторы $$\Lambda _{1}, \Lambda _{2}$$ коммутативны, и первым — если нет.

8.2.3. Общая формулировка методов расщепления

Заменим локально - одномерные дифференциальные операторы $$\mathbf{A}_i$$ разностными операторами на каждом шаге по времени $$t_n \le t \le t_{{n} + 1}$$. Представим схему расщепления в следующем общем виде:

$$\begin{gather*} \frac{{u^{{n} + 1/{N}} - u^{n}}}{\tau} + {{{\Lambda}}}_{10}u^{n} + {{{\Lambda}}}_{11} u^{{n} + 1/{N}} = 0, \\ \frac{{u^{{n} + 2/{N}} - u^{{n} + 1/{N}}}}{\tau } + {{{\Lambda}}}_{20} u^{n} + {{{\Lambda}}}_{21} u^{{n} + 1/{N}} + {{{\Lambda}}}_{22} u^{{n} + 2/{N}} = 0, \\ \ldots \\ \frac{u^{n + 1} - u^{n + \frac{N - 1}{N}}}{\tau} + {{{\Lambda}}}_{{N}0} u^{n} + {{{\Lambda}}}_{{N}1} u^{{n} + 1/2} + \ldots + {{{\Lambda}}}_{{N}{N}} u^{{n} + 1} = 0. \end{gather*}$$

Условие устойчивости такой схемы расщепления будет

$$\left\|{C_i C_{{i} - 1} \ldots C_1}\right\| \le 1 + c {\tau}, \quad c = const,$$

где $$C_i = {(\mathbf{E} + {\tau}\Lambda_{ii})}^{- 1} (\mathbf{E} + \tau \Lambda_{i, i - 1}), i = 1, 2, \ldots , N$$ Двухслойная схема расщепления с весовыми коэффициентами представлена в виде

$$\begin{gather*} \frac{{u^{{n} + 1/{N}} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \left[{(1 - {\sigma})u^{n} + {\sigma}u^{{n} + 1/{N}}}\right] = 0, \\ \frac{{u^{{n} + 2/{N}} - u^{{n} + 1/{N}}}}{\tau } + {\mathbf{\Lambda}}_2 \left[{(1 - \sigma )u^{{n} + 1/{N}} + \sigma u^{{n} + 2/{N}}}\right] = 0, \\ \ldots \\ \frac{{u^{{n} + 1} - u^{{n} + \frac{{{N} - 1}}{{N}}}}}{\tau} + {\mathbf{\Lambda}}_n \left[ {(1 - {\sigma})u^{{n} + \frac{{{N} - 1}}{{N}}} + {\sigma}u^{{n} + 1}}\right] = 0. \end{gather*}$$

Если в этой схеме расщепления положить веса верхнего и нижнего слоев по времени равными, $$\gamma = 0, 5$$, то в случае коммутирующих операторов $$\Lambda _{i}$$ (каждый такой разностный оператор аппроксимирует со вторым порядком соответствующий локально - одномерный дифференциальный оператор) схема будет иметь второй порядок аппроксимации и по времени. Если при этом каждый оператор $${\mathbf{\Lambda}}_i \ge 0$$, то схема будет абсолютно устойчивой.

8.2.4. Схемы расщепления для уравнения теплопроводности

Рассматриваем нестационарное уравнение теплопроводности

$$$ \frac{{\partial}u}{{\partial}t} - \Delta u = 0, t \in \Omega_t , \left\{{x, y, z}\right\} \in \Omega . $$$

Здесь оператор Лапласа определен как $$\Delta = {\partial}^2 / {\partial}x^2 + {\partial}^2 / {\partial}y^2 + {\partial}^2 / {\partial}z^2$$. Его также можно записать в виде суммы трех локально - одномерных операторов $$\mathbf{A} = \mathbf{A}_x + \mathbf{A}_y + \mathbf{A}_{z}$$. Соответствующие разностные операторы будут $${\mathbf{\Lambda}} = {\mathbf{\Lambda}}_{xx} + {\mathbf{\Lambda}}_{yy} + {\mathbf{\Lambda}}_{zz}$$, где $$$ \Lambda_{xx} = \frac{u_{m - 1, jk} - 2u_{mjk} + u_{m + 1, jk}}{h_x^2} $$$ аналогично определяются операторы вычисления второй разностной производной и по остальным направлениям $$\Lambda _{yy}, \Lambda _{zz}$$ Локально - одномерная схема для уравнения теплопроводности будет

$$\begin{gather*} \frac{{u^{{n} + 1/3} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_{xx}u^{{n} + 1/3} = 0, \\ \frac{{u^{n + 2/3} - u^{{n} + 1/3}}}{\tau} + {\mathbf{\Lambda}}_{yy}u^{{n} + 2/3} = 0, \\ \frac{{u^{{n} + 1} - u^{{n} + 2/3}}}{\tau} + {\mathbf{\Lambda}}_{zz}u^{{n} + 1} = 0. \end{gather*}$$

Для повышения порядка аппроксимации можно использовать схему с весами

$$\begin{gather*} \frac{{u^{{n} + 1/3} - u^{n}}}{\tau } + {\mathbf{\Lambda}}_{xx} \left[{(1 - {\sigma})u^n + {\sigma}u^{n + 1/3}}\right] = 0, \\ \frac{{u^{{n} + 2/3} - u^{{n} + 1/3}}}{\tau} + {\mathbf{\Lambda}}_{yy} \left[{(1 - {\sigma})u^{{n} + 1/3} + {\sigma}u^{{n} + 2/3}}\right] = 0, \\ \frac{{u^{{n} + 1} - u^{{n} + 2/3}}}{\tau} + {\mathbf{\Lambda}}_{zz}[(1 - {\sigma})u^{{n} + 2/3} + {\sigma}u^{{n} + 1} ] = 0. \end{gather*}$$

8.3. Методы двуциклического покомпонентного расщепления

Для этих методов отсутствует требование коммутативности операторов $$\Lambda _{i}.$$

Будем рассматривать численное решение (8.1) не на одном шаге по времени, отрезке [tn, tn + 1], а на двух последовательных шагах [tn - 1, tn + 1]. Пусть теперь разностные локально - одномерные операторы зависят явно от времени, тогда они определены в середине отрезка $${\mathbf{\Lambda}}_i = {\mathbf{\Lambda}}_i (t^{n} )$$. Запишем схему расщепления:

$$\begin{gather*} \frac{{u^{{n} - 1/2} - u^{{n} - 1}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{{n} - 1/2} + u^{{n} - 1}}}{2} = 0, \\ \frac{{u^{n} - u^{{n} - 1/2}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n} + u^{{n} - 1/2}}}{2} = 0, \\ \frac{{u^{{n} + 1/2} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{{n} + 1/2} - u^{n}}}{2} = 0, \\ {\frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{{n} + 1} - u^{{n} + 1/2}}}{2} = 0.} \end{gather*}$$

В операторной форме этот метод записывается как $$u^{n + 1} = {\mathbf{T}}^{n}u^{n - 1}$$, где введено обозначение

$$\begin{gather*} {\mathbf{T}}^{n} = \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right) \left({{\mathbf{E}} + \frac{\tau}{2} {\mathbf{\Lambda}}_2}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right) \times \\ \times \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)^{- 1} \left({{\mathbf{E}}- \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right) \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right) = \mathbf{E}{- } 2 {\tau}{\Lambda}+ \frac{{(2 \tau)}^2}{2}{(\Lambda)}^2 + \ldots \end{gather*}$$

Если локально - одномерные операторы положительны $${\mathbf{A}}_i (t) > 0$$, то при достаточной гладкости решения и элементов матриц $${\mathbf{A}}_i (t)$$ схема (8.6) абсолютно устойчива и аппроксимирует (8.1) со вторым порядком.

Для неоднородного дифференциального уравнения $$$ \frac{{\partial u}}{{\partial}t} + {\mathbf{A}}u = f $$$ разностная аппроксимация метода расщепления может быть представлена в виде

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - 1/2} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - 1}, \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)(u^{n} - {\tau}f^n) = \left({{\mathbf{E}} - \frac{\tau }{2}{\mathbf{\Lambda}}_2}\right)u^{n - 1/2}, \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)u^{n + 1/2} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)(u^{n} + {\tau}f^{n} ), \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + 1} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + 1/2}, \end{gather*}$$

где fn = f(tn).

В операторной форме записи решение неоднородной задачи имеет вид: $$u^{n + 1} = {\mathbf{T}}^{n} u^{n - 1} + 2 {\tau}{\mathbf{T}}_1^{n}{\mathbf{T}}_2^{n} f^{n}$$, где введены обозначения $$$ {\mathbf{T}}^{n} = {\mathbf{T}}_1^{n}{\mathbf{T}}_2^{n}{\mathbf{T}}_2^{n}{\mathbf{T}}_1^{n}$$ $${\mathbf{T}}_i^{n} = \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_i^{n}}\right)^{- 1} \left({{\mathbf{E}} - \frac{\tau }{2}{\mathbf{\Lambda}}_i^{n}}\right). $$$ Представим разностную аппроксимацию неоднородного дифференциального уравнения с помощью последовательного применения операторов $$\Lambda _{1}, \dots , \Lambda _{n}$$, $$\Lambda _{n}, \dots , \Lambda _{1}:$$

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - \frac{N - 1}{N}} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n - 1}, \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)(u^{n} - {\tau}f^{n} ) = \left({{\mathbf{E}} - \frac{\tau }{2}{\mathbf{\Lambda}}_n }\right)u^{n - 1/N}, \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)u^{n + 1/N} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)(u^{n} + {\tau}f^{n} ), \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + 1} = \left({{\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)u^{n + \frac{{N - 1}}{N}} . \end{gather*}$$

Рассмотрим примеры использования метода двуциклического расщепления для некоторых задач математической физики.

Пример 1. Трехмерное нестационарное уравнение диффузии, область интегрирования — параллелепипед. Полагаем, что в вертикальном направлении (ось 0z ) коэффициент диффузии в вертикальной плоскости $$\gamma$$ зависит от координаты, что характерно для задач геофизики, $$\mu$$ — коэффициент диффузии в горизонтальной плоскости. Задача может быть представлена в виде

$$$ \frac{{\partial}u}{{\partial}t} = \frac{\partial }{{{\partial}z}} \gamma \frac{{\partial}u}{{{\partial}z}} + \mu \Delta u + f. $$$

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

$$$ \frac{{{\partial}u_1}}{{\partial}t} = \frac{{\partial}}{{{\partial}z}} \gamma \frac{{{\partial}u_1}}{{{\partial}z}} + f, $$$

она описывает диффузию в вертикальной плоскости. Вторую и третью задачи запишем

$$$ \frac{{{\partial}u_2}}{{\partial}t} = \mu \frac{{{\partial}^2 u_2}} {{{\partial}x^2}}, \frac{{{\partial}u_3 }}{{\partial}t} = \mu \frac{{{\partial}^2 u_3 }} {{{\partial}y^2}}. $$$

Теперь рассмотрим разностную аппроксимацию исходного дифференциального уравнения

$$\begin{gather*} \frac{{\partial}u}{{\partial}t} + ({\mathbf{\Lambda}}_1 + {\mathbf{\Lambda}}_2 + {\mathbf{\Lambda}}_3 )u = f, \mbox{ где} \\ {\mathbf{\Lambda}}_1 u = -\frac{\mu}{h_x^2}(u_{m, j, k + 1} - 2u_{mjk}+ u_{m, j, k - 1} ),\\ {\mathbf{\Lambda}}_2 u = -\frac{\mu}{h_y^2}(u_{m, j - 1, k} - 2u_{mjk} + u_{m, j + 1, k}),\\ {\mathbf{\Lambda}}_3 u = \frac{1}{{h_{z}}} \left[{- \frac{{\gamma_{m + 1/2, jk}}}{{h_{z}}}(u_{m + 1, jk} - u_{mjk} ) + \frac{{\gamma_{m - 1/2}}}{{h_{z}}} (u_{mjk} - u_{m - 1, jk} )}\right]. \end{gather*}$$

Разностная схема двуциклического покомпонентного расщепления приобретает вид

$$\begin{gather*} \frac{{u^{n + 1/6} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1/6} + u^{n}}}{2} = 0, \\ \frac{{u^{n + 2/6} - u^{n + 1/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 2/6} + u^{n + 1/6}}}{2} = 0, \\ \frac{{u^{n + 3/6} - u^{n + 2/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 3/6} + u^{n + 2/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 4/6} - u^{n + 3/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 4/6} + u^{n + 3/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 5/6} - u^{n + 4/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 5/6} + u^{n + 4/6}}}{2} = 0, \\ \frac{{u^{n + 1} - u^{n + 5/6}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1} + u^{n + 5/6}}}{2} = 0. \end{gather*}$$

Пример 2. Сопряженное нестационарное уравнение переноса и диффузии

$$$ \frac{{\partial}u}{{\partial}t} + ({\mathbf{A}}_1 + {\mathbf{A}}_2 + {\mathbf{A}}_3 )u = f, $$$

где операторы определены как

$$$ {\mathbf{A}}_1 = \frac{{{\partial}(v_{1} u)}}{{\partial}x} - \mu \frac{{{\partial}^2 u}}{{{\partial}x^2}}, {\mathbf{A}}_2 = \frac{{{\partial}(v_2u)}}{{\partial}y} - \mu \frac{{{\partial}^2 u}}{{{\partial}y^2}}, {\mathbf{A}}_3 = \frac{{{\partial}(v_3u)}}{{{\partial}z}} - \frac{{\partial}}{{{\partial}z}} \gamma \frac{{\partial}u}{{{\partial}z}} + {\sigma}u. $$$

Здесь v1, v2, v3 — компоненты вектора скорости, u — концентрация субстанции, $$\sigma$$ — коэффициент поглощения субстанции внешней средой, $$\sigma > 0$$. Соответствующую схему расщепления представим в виде

$$\begin{gather*} \frac{{u^{n + 1/6} - u^{n}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1/6} + u^{n}}}{2} = 0, \\ \frac{{u^{n + 2/6} - u^{n + 1/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 2/6} + u^{n + 1/6}}}{2} = 0, \\ \frac{{u^{n + 3/6} - u^{n + 2/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 3/6} + u^{n + 2/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 4/6} - u^{n + 3/6}}}{\tau} + {\mathbf{\Lambda}}_3 \frac{{u^{n + 4/6} + u^{n + 3/6}}}{2} = \frac{{f^{n + 1/2}}}{2}, \\ \frac{{u^{n + 5/6} - u^{n + 4/6}}}{\tau} + {\mathbf{\Lambda}}_2 \frac{{u^{n + 5/6} + u^{n + 4/6}}}{2} = 0, \\ \frac{{u^{n + 1} - u^{n + 5/6}}}{\tau} + {\mathbf{\Lambda}}_1 \frac{{u^{n + 1} + u^{n + 5/6}}}{2} = 0, \end{gather*}$$

где разностные операторы аппроксимируют соответствующие локально - одномерные дифференциальные операторы. Так для рассматриваемой задачи

$$\begin{gather*} {\mathbf{\Lambda}}_1 u = - \frac{\mu }{{h_x^2}}(u_{m + 1, jk} - 2u_{mjk} + u_{m - 1j, k} ) + \frac{{(v_1 u)_{m + 1, jk} - (v_1 u)_{m - 1, jk}}}{{2h_x }}, \\ {\mathbf{\Lambda}}_2 u = - \frac{\mu }{{h_y^2}}(u_{m, j + 1, k} - 2u_{mjk} + u_{m, j - 1, k} ) + \frac{{(v_2 u)_{m, j + 1, k} - (v_2 u)_{m, j - 1, k}}}{{2h_y}}, \\ {\mathbf{\Lambda}}_3 u = \frac{1}{{h_{z}}} \left[{- \frac{{\gamma_{m + 1/2, jk}}}{{h_{z}}}(u_{m + 1, jk} - u_{mjk} ) + \frac{{\gamma_{m - 1/2}}}{h}(u_{mjk} - u_{m - 1, jk} )}\right] + \\ + \frac{{(v_3 u)_{m, j, k + 1} - (v_3 u)_{m, j, k - 1}}}{{2h_{z}}} + {\sigma}u_{mjk} . \end{gather*}$$

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

$$$ \frac{{\partial}u}{{\partial}t} + \frac{{{\partial}(v_1 , u)}} {{\partial}x} + \frac{{{\partial}(v_2 , u)}}{{\partial}y} - \mu \Delta u - {\sigma}u = f, $$$

в котором компоненты скорости движения среды v1, v2 удовлетворяют уравнению неразрывности:

$$$ \frac{{{\partial}v_1}}{{\partial}x} + \frac{{{\partial}v_2}} {{\partial}y} = 0. $$$

Первый этап расщепления задачи по физическим процессам связан с переносом, на нем решается разностный аналог уравнения

$$$ \frac{{{\partial}u_1}}{{\partial}t} + \frac{{{\partial}(v_1 , u)}} {{\partial}x} + \frac{{{\partial}(v_2 , u)}}{{\partial}y} = 0. $$$

Второй этап расщепления описывает процессы диффузии и поглощения субстанций

$$$ \frac{{{\partial}u_2}}{{\partial}t} - \mu \Delta u + {\sigma}u_2 = f. $$$

Пример 3. Расщепление по физическим процессам системы уравнений газовой динамики (метод крупных частиц).

Система

$$\begin{gather*} \frac{{\partial}\rho}{{\partial}t} + div ({\rho}{\mathbf{v}}) = 0, \\ \frac{{\partial}({\rho}u_1)}{{\partial}t} + div ({\rho}{\mathbf{v}}u_1 ) + \frac{{\partial P}}{{\partial}x} = 0, \\ \frac{{{\partial}({\rho}u_2 )}}{{\partial}t} + div ({\rho}{\mathbf{v}}u_2 ) + \frac{{\partial P}}{{\partial}y} = 0, \\ \frac{{{\partial}({\rho}e)}}{{\partial}t} + div ({\rho}e{\mathbf{v}}) + \div (P{\mathbf{v}}) = 0, \\ P = P({\rho}, \varepsilon ), e = \varepsilon + \frac{{u_1^2 + u_{2}^{2}}}{2}, \end{gather*}$$

u1 , u2 — компоненты вектора скорости $$\mathbf{v}$$, P — давление газа, $$\rho$$ — плотность, $$\varepsilon$$ — внутренняя энергия.

Первый (Эйлеров) этап. Измеряются лишь величины, относящиеся к ячейке в целом, а жидкость внутри каждой ячейки сетки считается моментально замороженной. Расчет производится по формулам, аппроксимирующим уравнения

$$\begin{gather*} \frac{{{\partial}{\rho}}}{{\partial}t} = 0, \\ {\rho}\frac{{{\partial}u_1}}{{\partial}t} + \frac{{\partial}p}{{\partial}x} = 0, {\rho}\frac{{{\partial}u_2}}{{\partial}t} + \frac{{\partial}p}{{\partial}y} = 0, \\ {\rho}\frac{{\partial e}}{{\partial}t} + div (P{v}) = 0. \end{gather*}$$

На втором (Лагранжевом) этапе происходит движение газа массы через границы эйлеровых ячеек и перераспределение массы, импульса, энергии по пространству; определяются поля параметров течения газа. Аппроксимируется система уравнений

$$\begin{gather*} \frac{{{\partial}\rho }}{{\partial}t} + div ({\rho}{\mathbf{v}}) = 0, \\ \frac{{{\partial}({\rho}u_1 )}}{{\partial}t} + div ({\rho}u_1 {\mathbf{v}}) = 0, \frac{{{\partial}({\rho}u_2 )}}{{\partial}t} + div ({\rho}u_2 {\mathbf{v}}) = 0, \\ \frac{{{\partial}({\rho}e)}}{{\partial}t} + div ({\rho}e{\mathbf{v}}) = 0. \end{gather*}$$

8.4. Методы расщепления с факторизацией оператора

8.4.1. Факторизованная схема расщепления

Пусть для решения дифференциальной задачи

$$$ \mathbf{B} \frac{{\partial}u}{{\partial}t} + \mathbf{A}u = f, u(t_0 ) = u_0 $$$

используется разностная схема $${\mathbf{B}}u^{n + 1} = F^{n}$$, где $$F^{n} = ({\mathbf{B}} - {\tau}{\mathbf{A}})u^{n} + {\tau}f^{n}$$, n = 0, 1, ... Пусть для вычисления Fn затрачивается O(N) действий, число арифметических операций пропорционально числу узлов сетки N. Такие разностные операторы называются экономичными.

Пусть $$\mathbf{B}_i$$ (i = 1, 2, ..., N) - экономичные разностные операторы, такие, что $${\mathbf{B}}_i v = F$$. Назовем схему разностной схемой с факторизованным оператором $$\mathbf{B}$$, если возможно его представление в виде

$$\mathbf{B} = \mathbf{B}_1 \mathbf{B}_2 \ldots \mathbf{B}_n.$$

Эта схема будет также экономичной, так как для решения разностного уравнения по - прежнему потребуется O(N) действий. В самом деле, решение уравнения

$$\mathbf{B}_1 \mathbf{B}_2 \ldots \mathbf{B}_n u^{n + 1} = F^{n} $$

может быть найдено в результате последовательного решения p уравнений

$$\begin{gather*} {\mathbf{B}}_1 u_1 = F^{n}, \\ {\mathbf{B}}_2 u_2 = u_1 , \\ {\mathbf{B}}_3 u_3 = u_2, \\ \ldots \\ {\mathbf{B}}_i u_i = u_{i - 1}, \end{gather*} $$

здесь i = 2, 3, ..., N. Тогда un + 1 = un. В записи задачи введены обозначения $$$ u_1 = u^{n + 1/N}, \ldots , u_i = u^{n + i/N}, \ldots , u_{N - 1} = u^{n + \frac{{N - 1}}{N}} $$$ — промежуточные значения.

Схемы с факторизованным оператором иногда называются также факторизованными схемами. Устойчивая схема с факторизованным оператором $$\mathbf{B}$$, которая представляет собой произведение конечного числа операторов $$\mathbf{B}_1, \ldots , \mathbf{B}_n$$, является экономичной схемой.

Пример. Метод переменных направлений (продольно - поперечная схема). Приведем запись схемы для решения линейного двумерного уравнения теплопроводности. Расчетные формулы есть

$$\begin{gather*} \frac{{u^{n + 1/2} - u^{n}}}{{{{1 \over 2}} \tau }} - ({\mathbf{\Lambda}}_1 u^{n + 1/2} + {\mathbf{\Lambda}}_2 u^{n} ) = f^{n}, \\ \frac{{u^{n + 1} - u^{n + 1/2}}}{{{{1 \over 2}} {\tau}}} - ({\mathbf{\Lambda}}_1 u^{n + 1/2} + {\mathbf{\Lambda}}_2 u^{n + 1} ) = f^{n} . \end{gather*}$$

Тогда, исключая un + 1/2, получим в операторной форме записи

$$$ \left(\mathbf{E} - \frac{\tau}{2}\mathbf{\Lambda}_1\right) \left(\mathbf{E} - \frac{\tau}{2}\mathbf{\Lambda}_2\right) u^{n + 1} = \left(\mathbf{E} + \frac{\tau}{2}\mathbf{\Lambda}_1\right) \left(\mathbf{E} + \frac{\tau}{2} \mathbf{\Lambda}_2\right) u^{n}, $$$

или $$\mathbf{B}_1\mathbf{B}_2u^{n + 1} = F^{n}$$, где $$$ {\mathbf{B}}_1 = {\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_1 $$$, $$$ {\mathbf{B}}_2 = {\mathbf{E}} - \frac{\tau}{2}{\mathbf{\Lambda}}_2 $$$, $$$ F^{n} = \left(\mathbf{E} + \frac{\tau}{2}\mathbf{\Lambda}_1\right) \left(\mathbf{E} + \frac{\tau}{2}\mathbf{\Lambda}_2\right)u^{n} $$$. Разностная схема может быть представлена в виде факторизованный схемы расщепления:

$${\mathbf{B}}_1 u_1 = F^{n}, {\mathbf{B}}_2 u^{n + 1} = u_1.$$

8.4.2. Неявная схема расщепления с приближенной факторизацией

Рассмотрим неявную разностную схему

$$$ {\frac{{u^{n + 1} - u^{n}}}{\tau} + {\Lambda}u^{n + 1} = 0, n = 0, 1, \ldots , {\mathbf{\Lambda}} = \sum\limits_{i = 1}^{N}{{\mathbf{\Lambda}}_i }, \mathbf{\Lambda}_i > 0.} $$$

Представим разностную схему (8.7) в виде

$${({\mathbf{E}} + {\tau}{\mathbf{\Lambda}})u^{n + 1} = u^{n} .}$$

Факторизуем разностную схему (8.8) приближенно с точностью до членов порядка $$O(\tau ^{2}).$$ Для этого заменим в (8.8) оператор $${\mathbf{E}} + {\tau}{\mathbf{\Lambda}}$$ на факторизованный

$$({\mathbf{E}} + \tau {\mathbf{\Lambda}}_1 )({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_2 ) \ldots ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_n ) = {\mathbf{E}} + {\tau}{\mathbf{\Lambda}} + {\tau}^2 {\mathbf{R}},$$

где введено обозначение

$${\mathbf{R}} = \sum\limits_{i < j}{{\mathbf{\Lambda}}_i {\mathbf{\Lambda}}_j + {\tau}\sum\limits_{i < j < k}{{\mathbf{\Lambda}}_i }}{\mathbf{\Lambda}}_j {\mathbf{\Lambda}}_k + \ldots + {\tau}^{n - 2}{\mathbf{\Lambda}}_1 \ldots {\mathbf{\Lambda}}_n.$$

Врезультате приходим к неявной схеме с приближенной факторизацией

$${\mathbf{B}}u^{n + 1} = u^{n}, \mathbf{B} = {\mathop \Pi\limits_{i = 1}^{n}{\mathbf{B}_i}}, {\mathbf{B}}_i = {\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_i,$$

или

$$\begin{gather*} ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_1 )u^{n + 1/N} = u^{n}, \\ ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_2 )u^{n + 2/N} = u^{n + 1/N}, \\ \ldots \\ ({\mathbf{E}} + {\tau}{\mathbf{\Lambda}}_n )u^{n + 1} = u^{n + \frac{{N - 1}}{N}} . \end{gather*}$$

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

8.4.3. Метод "предиктор - корректор"

Основная идея методов типа "предиктор - корректор" заключается в следующем. На каждом отрезке [tn, tn + 1] задача решается в два приема: сначала по схеме первого порядка аппроксимации и со значительным запасом устойчивости находится решение в момент времени $$t^{n + 1/2} = t^{n} + \tau /2$$ — предиктор. После этого на втором этапе расписывается исходное уравнение по схеме более высокого порядка аппроксимации (чаще всего, второго) — корректор. Основная идея семейств таких методов близка к идее построения методов типа Рунге - Кутты для обыкновенных дифференциальных уравнений.

Представим эту схему как следующую схему расщепления:

$$\begin{gather*} \frac{{u^{n + 1/4} - u^{n}}}{{{\tau}/2}} + {\mathbf{\Lambda}}_1u^{n + 1/4} = 0, \\ \frac{{u^{n + 1/2} - u^{n + 1/4}}}{{{\tau}/2}} + {\mathbf{\Lambda}}_2 u^{n + 1/2} = 0, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = 0. \end{gather*} $$

Если в этой схеме расщепления исключить un + 1/4, то получим последовательность расчетных формул

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right) \left({{\mathbf{E}} + \frac{\tau}{2} {\mathbf{\Lambda}}_2}\right)u^{n + 1/2} = {\varphi}^{n}, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = 0, \end{gather*}$$

далее, исключив, $${\mathbf{\Lambda}}u^{n + 1/2}$$, получим

$$$ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)^{- 1} \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_1}\right)^{- 1}u^{n} = 0. $$$

Если для разностных операторов выполнены условия $$$ \frac{\tau}{2} \left\| {{\mathbf{\Lambda}}_i }\right\| < 1, $$$ $${\mathbf{\Lambda}}_1 \ge 0$$, $${\mathbf{\Lambda}}_2 \ge 0$$, а коэффициенты разностной схемы явно не зависят от времени, то при достаточной гладкости решения дифференциальной задачи разностная схема абсолютно устойчива и аппроксимирует исходную задачу со вторым порядком.

Далее рассмотрим случай, когда оператор $$\mathbf{\Lambda}$$ представляется в виде суммы операторов $$\Lambda _{i}:$$$${\mathbf{\Lambda}} = \sum\limits_i {{\mathbf{\Lambda}}_i } .$$ Пусть все эти разностные операторы положительны. Метод "предиктор - корректор" можно записать в виде последовательности расчетных формул

$$\begin{gather*} \left({{\mathbf{E}} + \frac{\tau}{2} {\mathbf{\Lambda}}_1}\right)u^{n + 1/2N} = u^{n} + \frac{\tau}{2}f^{n + 1/2}, \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_2}\right)u^{n + 2/2N} = u^{n + 1/2N}, \\ \ldots \\ \left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_n }\right)u^{n + 1/2} = u^{n + \frac{N}{{2N}}}, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = f^{n + 1/2} . \end{gather*}$$

Эта последовательность после исключения промежуточных этапов сводится к одному разностному уравнению

$$$ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \prod\limits_{i = N}^1 {\left({{\mathbf{E}} + \frac{\tau}{2}{\mathbf{\Lambda}}_i } \right)} ^{- 1} (u^{n} + \frac{\tau}{2}f^{n + 1/2} ) = f^{n + 1/2} . $$$

Приведем пример построения такой схемы. Для нестационарного трехмерного уравнения теплопроводности

$$$ \frac{{\partial}u}{{\partial}t} - a^2 \sum\limits_{i = 1}^3 {\frac{{{\partial}^2 u}}{{{\partial}x_i^2}} = 0} $$$

получим следующую разностную схему типа "предиктор - корректор"

$$\begin{gather*} \frac{{u^{n + 1/6} - u^{n}}}{{{\tau}/2}} + {\mathbf{\Lambda}}_1 u^{n + 1/6} = 0, \\ \frac{u^{n + 1/3} - u^{n + 1/6}}{{\tau}/2} + {\mathbf{\Lambda}}_2 u^{n + 1/3} = 0, \\ \frac{u^{n + 1/2} - u^{n + 1/3}}{{\tau}/2} + {\mathbf{\Lambda}}_3 u^{n + 1/2} = 0, \\ \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}}u^{n + 1/2} = 0. \end{gather*}$$

Далее, исключая промежуточные слои, получим схему, записанную в каноническом виде

$$\begin{gather*} \frac{{u^{n + 1} - u^{n}}}{\tau} + {\mathbf{\Lambda}} \frac{{u^{n + 1} + u^{n}}}{2} + \frac{{{\tau}^2}}{4}({\mathbf{\Lambda}}_1 {\mathbf{\Lambda}}_2 + {\mathbf{\Lambda}}_1 {\mathbf{\Lambda}}_3 + {\mathbf{\Lambda}}_2 {\mathbf{\Lambda}}_3 ) \frac{{u^{n + 1} - u^{n}}}{\tau} + \\ + \frac{{{\tau}^3 }}{8}{\mathbf{\Lambda}}_1 {\mathbf{\Lambda}}_2 {\mathbf{\Lambda}}_3 \frac{{u^{n + 1} - u^{n}}}{\tau} = 0. \end{gather*}$$

Схема абсолютно устойчива (для коммутирующих операторов), имеет второй порядок аппроксимации по $$\tau$$ и hi. Конечно, при практическом решении задач на компьютере используется именно последовательность разностных операторов. Канонический вид схемы удобен для ее теоретического исследования.

Вернуться к учебному плану