Контрольная работа № 2 — все варианты

Уравнение и условия

Рассматривается двумерное параболическое уравнение с реакционным членом:

$$\dfrac{\partial u}{\partial t}=0{,}2\,\dfrac{\partial^2 u}{\partial x^2}+0{,}5\,\dfrac{\partial^2 u}{\partial y^2}-5tu,$$

граничные условия

$$\dfrac{\partial u}{\partial x}(t,0,y)=ty,\qquad \dfrac{\partial u}{\partial x}(t,1,y)=5ty,$$ $$u(t,x,0)=tx,\qquad u(t,x,1)=2tx,$$

начальное условие $u(t{=}0,x,y)=0$.

Здесь коэффициенты диффузии $\sigma_x=0{,}2$, $\sigma_y=0{,}5$, реакция $q(t)=-5t$. По направлению $x$ заданы условия 2-го рода (Неймана) на обоих концах, по направлению $y$ — условия 1-го рода (Дирихле).

Сетка и разностные операторы

Вводим равномерную сетку: $x_j=(j-1)h_x,\ j=\overline{1,N_x}$ ($h_x=1/(N_x-1)$); $y_k=(k-1)h_y,\ k=\overline{1,N_y}$ ($h_y=1/(N_y-1)$); $t^n=n\,\Delta t$. Сеточная функция $u_{j,k}^{n}\approx u(t^n,x_j,y_k)$. Вторые разности:

$$\Lambda_{xx}u_{j,k}=\dfrac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^{2}},\qquad \Lambda_{yy}u_{j,k}=\dfrac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^{2}}.$$

Идея схемы переменных направлений (СПН)

Полностью неявная по обоим направлениям схема даёт пятидиагональную матрицу — прогонка неприменима. Метод переменных направлений (схема Писмена–Рэкфорда) дробит шаг $\Delta t$ промежуточным слоём $t^{n+1/2}=t^n+\Delta t/2$ на две подсхемы. В каждой подсхеме неявен только один диффузионный оператор — тогда система трёхдиагональна и решается прогонкой; оператор по другому направлению берётся явно. Каждый полушаг имеет «длину» $\Delta t/2$, поэтому расщеплённые диффузионные операторы входят с множителем $\tfrac{\sigma}{2}\,\Delta t$. Схема устойчива безусловно; диффузионная часть имеет порядок $O(\Delta t^2+h_x^2+h_y^2)$.

Подсхема ① ($n\to n+1/2$): неявно по $x$, явно по $y$

$$\dfrac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t/2}=0{,}2\,\Lambda_{xx}u_{j,k}^{n+1/2}+0{,}5\,\Lambda_{yy}u_{j,k}^{n}-5t^{n+1/2}u_{j,k}^{n+1/2}.$$

Неизвестные слоя $n+1/2$ стоят при операторе $\Lambda_{xx}$ и при реакционном члене; член $\Lambda_{yy}$ взят с известного слоя $n$.

Подсхема ② ($n+1/2\to n+1$): неявно по $y$, явно по $x$

$$\dfrac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t/2}=0{,}2\,\Lambda_{xx}u_{j,k}^{n+1/2}+0{,}5\,\Lambda_{yy}u_{j,k}^{n+1}-5t^{n+1}u_{j,k}^{n+1}.$$

Теперь неявен оператор $\Lambda_{yy}$ (слой $n+1$), а $\Lambda_{xx}$ берётся с уже найденного слоя $n+1/2$.

Приведение подсхемы ① к виду для прогонки (прогонка по $j$, при фиксированном $k$)

Умножаем уравнение на $\Delta t/2$ и переносим все неизвестные слоя $n+1/2$ влево. Обозначим $r_x=\dfrac{0{,}2\,\Delta t}{2h_x^{2}}=\dfrac{0{,}1\,\Delta t}{h_x^{2}}$. Получаем трёхдиагональную по $x$ систему

$$c_j\,u_{j-1,k}^{n+1/2}+b_j\,u_{j,k}^{n+1/2}+a_j\,u_{j+1,k}^{n+1/2}=\xi_{j,k},$$

с коэффициентами (проверено символьно):

$$a_j=c_j=-\dfrac{0{,}2\,\Delta t}{2h_x^{2}}=-\dfrac{0{,}1\,\Delta t}{h_x^{2}},\qquad b_j=1+\dfrac{0{,}2\,\Delta t}{h_x^{2}}+\dfrac{5t^{n+1/2}\Delta t}{2},$$ $$\xi_{j,k}=u_{j,k}^{n}+\dfrac{0{,}5\,\Delta t}{2}\,\Lambda_{yy}u_{j,k}^{n}.$$

Внедиагональные $a_j,c_j$ отрицательны ($\Lambda_{xx}$ при переносе влево даёт $-r_x(u_{j+1}+u_{j-1})$); на диагонали $+2r_x=\dfrac{0{,}2\Delta t}{h_x^2}$ плюс единица из $\partial/\partial t$ плюс неотрицательный вклад реакции.

Проверка сходимости (устойчивости) прогонки для подсхемы ①

Достаточное условие корректности прогонки — диагональное преобладание $|b_j|>|a_j|+|c_j|$:

$$|a_j|+|c_j|=\dfrac{0{,}2\,\Delta t}{h_x^{2}}\;<\;1+\dfrac{0{,}2\,\Delta t}{h_x^{2}}+\dfrac{5t^{n+1/2}\Delta t}{2}=|b_j|\qquad (t\ge 0).$$

Неравенство выполнено всегда: единица от $\partial/\partial t$ уже обеспечивает строгое преобладание, а реакционный член $\dfrac{5t^{n+1/2}\Delta t}{2}\ge 0$ лишь усиливает диагональ. Прогонка устойчива.

Рекуррентное прогоночное соотношение (подсхема ①)

Ищем решение в виде $u_{j,k}^{n+1/2}=\alpha_j\,u_{j+1,k}^{n+1/2}+\beta_j$. Прогоночные коэффициенты:

$$\alpha_j=\dfrac{-a_j}{b_j+c_j\,\alpha_{j-1}},\qquad \beta_j=\dfrac{\xi_{j,k}-c_j\,\beta_{j-1}}{b_j+c_j\,\alpha_{j-1}}.$$

Стартовые коэффициенты $\alpha_1,\beta_1$ из левого ГУ по $x$ (2-го рода)

Левая граница $x=0$ — условие Неймана $\dfrac{\partial u}{\partial x}(t,0,y)=ty$. Аппроксимируем правой (односторонней) разностью первого порядка в узле $j=1$:

$$\dfrac{u_{2,k}-u_{1,k}}{h_x}=t\,y_k\;\Longrightarrow\;u_{1,k}=u_{2,k}-h_x\,t\,y_k.$$

Сравнивая с $u_{1,k}=\alpha_1 u_{2,k}+\beta_1$, получаем стартовые значения (берётся $t=t^{n+1/2}$):

$$\boxed{\;\alpha_1=1,\qquad \beta_1=-h_x\,t^{n+1/2}\,y_k\;}$$

Замечание о точности: такая односторонняя аппроксимация даёт на границе порядок $O(h_x)$, что формально понижает общий порядок схемы. Для согласования с $O(h_x^2)$ можно использовать аппроксимацию 2-го порядка через фиктивный узел: $\dfrac{u_{2,k}-u_{0,k}}{2h_x}=t y_k$ с подстановкой $u_{0,k}$ в уравнение узла $j=1$; тогда $\alpha_1=\dfrac{-2a_1}{b_1}$ и т.п. Ниже для краткости оставлен вариант 1-го порядка.

Правая граница $x=1$ ($j=N_x$) — тоже Нейман $\dfrac{u_{N_x,k}-u_{N_x-1,k}}{h_x}=5t y_k$. Подставляя $u_{N_x-1,k}=\alpha_{N_x-1}u_{N_x,k}+\beta_{N_x-1}$, замыкаем прогонку:

$$u_{N_x,k}^{n+1/2}=\dfrac{\beta_{N_x-1}+h_x\,5t^{n+1/2}y_k}{1-\alpha_{N_x-1}},$$

после чего обратным ходом $u_{j,k}^{n+1/2}=\alpha_j u_{j+1,k}^{n+1/2}+\beta_j$ для $j=N_x-1,\dots,1$. (Делитель $1-\alpha_{N_x-1}\neq0$, так как при $\alpha_1=1$ все последующие $\alpha_j<1$.)

Подсхема ② — приведение к прогонке по $k$ (при фиксированном $j$)

Аналогично умножаем на $\Delta t/2$, неизвестные слоя $n+1$ влево, $r_y=\dfrac{0{,}5\,\Delta t}{2h_y^{2}}=\dfrac{0{,}25\,\Delta t}{h_y^{2}}$. Трёхдиагональная по $y$ система $\tilde c_k u_{j,k-1}^{n+1}+\tilde b_k u_{j,k}^{n+1}+\tilde a_k u_{j,k+1}^{n+1}=\tilde\xi_{j,k}$:

$$\tilde a_k=\tilde c_k=-\dfrac{0{,}5\,\Delta t}{2h_y^{2}}=-\dfrac{0{,}25\,\Delta t}{h_y^{2}},\qquad \tilde b_k=1+\dfrac{0{,}5\,\Delta t}{h_y^{2}}+\dfrac{5t^{n+1}\Delta t}{2},$$ $$\tilde\xi_{j,k}=u_{j,k}^{n+1/2}+\dfrac{0{,}2\,\Delta t}{2}\,\Lambda_{xx}u_{j,k}^{n+1/2}.$$

Диагональное преобладание: $|\tilde a_k|+|\tilde c_k|=\dfrac{0{,}5\,\Delta t}{h_y^{2}}<1+\dfrac{0{,}5\,\Delta t}{h_y^{2}}+\dfrac{5t^{n+1}\Delta t}{2}=|\tilde b_k|$ — выполнено всегда. Рекуррентное соотношение $u_{j,k}^{n+1}=\tilde\alpha_k u_{j,k+1}^{n+1}+\tilde\beta_k$ с теми же формулами для $\tilde\alpha_k,\tilde\beta_k$.

Стартовые коэффициенты по $y$ (ГУ 1-го рода)

По направлению $y$ заданы условия Дирихле, поэтому граничные узлы известны прямо: $u_{j,1}^{n+1}=t^{n+1}x_j$ (нижняя грань $y=0$) и $u_{j,N_y}^{n+1}=2t^{n+1}x_j$ (верхняя грань $y=1$). Прогонка по $k$ стартует с точного граничного значения, что эквивалентно

$$\tilde\alpha_1=0,\qquad \tilde\beta_1=t^{n+1}x_j,$$

а на правом конце сразу берётся $u_{j,N_y}^{n+1}=2t^{n+1}x_j$ как замыкающее значение для обратного хода.

Итог по алгоритму одного шага

1. Подсхема ①: для каждой строки $k$ прогонкой по $j$ найти $u_{\cdot,k}^{n+1/2}$ (ГУ по $x$ — 2-го рода, $\alpha_1=1,\ \beta_1=-h_x t^{n+1/2}y_k$).
2. Подсхема ②: для каждого столбца $j$ прогонкой по $k$ найти $u_{j,\cdot}^{n+1}$ (ГУ по $y$ — 1-го рода, граничные узлы заданы явно, $\tilde\alpha_1=0,\ \tilde\beta_1=t^{n+1}x_j$).
3. Обе подсхемы трёхдиагональны и устойчивы (диагональное преобладание выполняется всегда, реакция $-5tu$ при $t\ge0$ лишь усиливает диагональ).

Уравнение и условия

Двумерное уравнение первого порядка (перенос):

$$\frac{\partial u}{\partial t}=2\,\frac{\partial u}{\partial x}-0{,}05\,\frac{\partial u}{\partial y}.$$

ГУ и НУ: $u(t,0,y)=0$, $u(t,1,y)=t$, $u(t,x,0)=0$, $u(t,x,1)=e^{tx}$, $u(0,x,y)=e^{x}$. Требуется выбрать соответствующие ГУ и записать неявную схему методом дробных шагов.

1. Каноническая форма

Переносим конвективные члены влево: $\dfrac{\partial u}{\partial t}-2\dfrac{\partial u}{\partial x}+0{,}05\dfrac{\partial u}{\partial y}=0$. Сравнивая с $\dfrac{\partial u}{\partial t}+v_1\dfrac{\partial u}{\partial x}+v_2\dfrac{\partial u}{\partial y}=f$:

$$v_1=-2,\qquad v_2=+0{,}05,\qquad f\equiv0.$$

Случай 3 ($v_1<0,\,v_2>0$).

2. Выбор разностей и ГУ (правило по знаку $v$)

• $v_1=-2<0$ — по $x$ правая разность $\dfrac{u_{j+1,k}-u_{j,k}}{h_x}$, правое ГУ при $x=1$: $u(t,1,y)=t$.
• $v_2=+0{,}05>0$ — по $y$ левая разность $\dfrac{u_{j,k}-u_{j,k-1}}{h_y}$, левое ГУ при $y=0$: $u(t,x,0)=0$.

Используемый набор: $u(0,x,y)=e^{x}$, $u(t,1,y)=t$, $u(t,x,0)=0$. Условия $u(t,0,y)=0$ и $u(t,x,1)=e^{tx}$ не используются: задача переноса формально переопределена ГУ на всех четырёх сторонах, но на «выходных» границах ($x=0$ для $v_1<0$ и $y=1$ для $v_2>0$) граничные условия игнорируются — схема сама вычисляет там значения по характеристикам. Сетка $u_{j,k}^n=u(t^n,x_j,y_k)$.

3. Неявная схема и расщепление

Полная неявная схема: $\dfrac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}+v_1\dfrac{u_{j+1,k}^{n+1}-u_{j,k}^{n+1}}{h_x}+v_2\dfrac{u_{j,k}^{n+1}-u_{j,k-1}^{n+1}}{h_y}=0.$ Метод дробных шагов делит $\Delta t$ точкой $t^{n+1/2}$ на две одномерные подсхемы:

① ($n\to n+1/2$, неявная по $x$): $\dfrac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}+v_1\dfrac{u_{j+1,k}^{n+1/2}-u_{j,k}^{n+1/2}}{h_x}=0.$

② ($n+1/2\to n+1$, неявная по $y$): $\dfrac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}+v_2\dfrac{u_{j,k}^{n+1}-u_{j,k-1}^{n+1}}{h_y}=0.$

Сложение подсхем: $\dfrac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}+L_x u^{n+1/2}+L_y u^{n+1}=0$ — согласуется с полной неявной схемой с точностью $O(\Delta t)$ (оператор по $x$ берётся на полуслое). Порядок аппроксимации $O(\Delta t,h_x,h_y)$; обе подсхемы абсолютно (безусловно) устойчивы.

4. Рекуррентное соотношение ① (по $x$)

Уравнение первого порядка ⇒ на полуслое связаны только два узла, прогонка не нужна. $u_{j,k}^{n+1/2}\!\left(1-v_1\tfrac{\Delta t}{h_x}\right)=u_{j,k}^{n}-v_1\tfrac{\Delta t}{h_x}u_{j+1,k}^{n+1/2}$, откуда

$$u_{j,k}^{n+1/2}=\dfrac{u_{j,k}^{n}-v_1\tfrac{\Delta t}{h_x}u_{j+1,k}^{n+1/2}}{1-v_1\tfrac{\Delta t}{h_x}}\ \xrightarrow{v_1=-2}\ \boxed{\,u_{j,k}^{n+1/2}=\dfrac{u_{j,k}^{n}+2\tfrac{\Delta t}{h_x}u_{j+1,k}^{n+1/2}}{1+2\tfrac{\Delta t}{h_x}}\,}$$

Знаменатель $1+2\tfrac{\Delta t}{h_x}>0$ всегда (абсолютная устойчивость подтверждена анализом фон Неймана: при $v_1<0$ имеем $|G|\le1$ для любого $\Delta t$). Счёт по $x$ убывающий: $j=N_x-1,\dots,1$ (обновляется и ребро $x=0$); старт $u_{N_x,k}^{n+1/2}=t^{n+1/2}=(n+\tfrac12)\Delta t$.

5. Рекуррентное соотношение ② (по $y$)

$u_{j,k}^{n+1}\!\left(1+v_2\tfrac{\Delta t}{h_y}\right)=u_{j,k}^{n+1/2}+v_2\tfrac{\Delta t}{h_y}u_{j,k-1}^{n+1}$, откуда

$$u_{j,k}^{n+1}=\dfrac{u_{j,k}^{n+1/2}+v_2\tfrac{\Delta t}{h_y}u_{j,k-1}^{n+1}}{1+v_2\tfrac{\Delta t}{h_y}}\ \xrightarrow{v_2=0{,}05}\ \boxed{\,u_{j,k}^{n+1}=\dfrac{u_{j,k}^{n+1/2}+0{,}05\tfrac{\Delta t}{h_y}u_{j,k-1}^{n+1}}{1+0{,}05\tfrac{\Delta t}{h_y}}\,}$$

Знаменатель $1+0{,}05\tfrac{\Delta t}{h_y}>0$ всегда (при $v_2>0$ с левой разностью $|G|\le1$ безусловно). Счёт по $y$ возрастающий: $k=2,\dots,N_y$ (обновляется и ребро $y=1$, ГУ при $y=1$ не используется); старт $u_{j,1}^{n+1}=0$ (ГУ при $y=0$).

6. Алгоритм

1. $u_{j,k}^{0}=e^{x_j}$.
2. Цикл по $n$: первый полушаг по $x$ (ГУ $u_{N_x,k}^{n+1/2}=t^{n+1/2}$, $j=N_x-1,\dots,1$ по ①), затем второй полушаг по $y$ (ГУ $u_{j,1}^{n+1}=0$, $k=2,\dots,N_y$ по ②).

Итог: неявная схема дробных шагов, порядок $O(\Delta t,h_x,h_y)$, абсолютно устойчива; решается прямыми рекуррентными соотношениями (без прогонки, т.к. уравнение первого порядка).

Условие

Привести к виду для метода установления со схемой Кранка–Николсона уравнение

$$\frac{du}{dx}+0{,}3\,\frac{d^{2}u}{dx^{2}}=3x^{2},\qquad \frac{du}{dx}(x{=}0)=0,\qquad \frac{du}{dx}(x{=}1)=1.$$

1. Приведение к стандартному виду $v\,u'=\sigma\,u''-k\,u+f$

Канонический вид курса (гл.10): вторая производная — справа с положительным коэффициентом $\sigma>0$, первая производная — слева. В исходном уравнении при $u''$ стоит $+0{,}3$. Если просто перенести $0{,}3\,u''$ вправо, получим $u'=-0{,}3\,u''+3x^{2}$, то есть $\sigma=-0{,}3<0$ — недопустимо. Поэтому оставляем $u''$ справа с плюсом и переносим $u'$ влево с тем знаком, который для этого нужен:

$$0{,}3\,u''=3x^{2}-u'\ \Longrightarrow\ -\,u'=0{,}3\,u''-3x^{2}.$$

Проверка: $-u'=0{,}3u''-3x^{2}\Rightarrow u'+0{,}3u''=3x^{2}$ — совпадает с исходным. Сравнивая с $v\,u'=\sigma\,u''-k\,u+f$, получаем

$$\boxed{v=-1,\qquad \sigma=0{,}3,\qquad k=0,\qquad f(x)=-3x^{2}.}$$

Так как $k=0$, для стационарной задачи достаточное условие сходимости прогонки $|a_j|+|c_j|\le|b_j|$ обращается в равенство $|a_j|+|c_j|=|b_j|$ (строгого диагонального преобладания нет — проверено прямым подсчётом) — прогонка напрямую неприменима. Поэтому используем метод установления. Поскольку $v=-1<0$, конвективный член $v\,u'$ аппроксимируется правой (по потоку, upwind) конечной разностью $\dfrac{u_{j+1}-u_j}{h}$.

2. Введение фиктивной производной по времени

Добавляем в левую часть фиктивную производную по «времени» (итерациям) со знаком плюс; искомая функция становится функцией двух переменных $u(x)\to\widetilde u(x,t)$:

$$\frac{\partial\widetilde u}{\partial t}+v\,\frac{\partial\widetilde u}{\partial x}=\sigma\,\frac{\partial^{2}\widetilde u}{\partial x^{2}}-k\,\widetilde u+f(x),\qquad\text{т.е.}\qquad \frac{\partial\widetilde u}{\partial t}-\frac{\partial\widetilde u}{\partial x}=0{,}3\,\frac{\partial^{2}\widetilde u}{\partial x^{2}}-3x^{2}.$$

Это параболическое уравнение. Граничные условия берутся из исходной стационарной задачи и от времени не зависят, поэтому при $t\to\infty$ решение «устанавливается»:

$$t\to\infty:\quad \widetilde u(x,t)\to [u(x)],\qquad \frac{\partial\widetilde u}{\partial t}\to 0.$$

Пошаговый переход $n\to n+1$ — это итерация, шаг $\Delta t$ — шаг итерации.

3. Схема Кранка–Николсона

Сетка по координате $x_j=(j-1)h$, $j=1,\dots,N$, $h=\dfrac{1}{N-1}$; «время» — итерации с шагом $\Delta t$. Пространственные операторы берутся как полусумма аппроксимаций на слоях $n$ и $n{+}1$; производная по времени — центральная относительно $(n{+}1/2)$. Так как $v<0$ — правая разность для $\partial u/\partial x$:

$$\frac{u_j^{n+1}-u_j^{n}}{\Delta t}+\frac{v}{2}\frac{u_{j+1}^{n+1}-u_j^{n+1}}{h}+\frac{v}{2}\frac{u_{j+1}^{n}-u_j^{n}}{h}=\frac{\sigma}{2}\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^{2}}+\frac{\sigma}{2}\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^{2}}+f(x_j),$$

где $v=-1$, $\sigma=0{,}3$, $f(x_j)=-3x_j^{2}$. Схема абсолютно устойчива; порядок $O(\Delta t^{2},h)$: по времени и диффузии — 2-й, по сносу (upwind 1-го порядка) — 1-й, что и ограничивает суммарный пространственный порядок до $O(h)$.

4. Приведение к виду для прогонки

Умножаем на $\Delta t$ и собираем неизвестные слоя $(n{+}1)$ слева. Получаем трёхдиагональную систему

$$a_j\,u_{j+1}^{n+1}+b_j\,u_j^{n+1}+c_j\,u_{j-1}^{n+1}=\xi_j^{n}.$$

Вклад конвективного (правая разность, слой $n{+}1$): $\dfrac{v}{2}\dfrac{\Delta t}{h}\bigl(u_{j+1}^{n+1}-u_j^{n+1}\bigr)$ — в $a_j$ даёт $+\dfrac{v}{2}\dfrac{\Delta t}{h}$, в $b_j$ даёт $-\dfrac{v}{2}\dfrac{\Delta t}{h}$. Вклад диффузии (слой $n{+}1$): $-\dfrac{\sigma}{2}\dfrac{\Delta t}{h^{2}}\bigl(u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}\bigr)$. Плюс единица от производной по времени в $b_j$. Итого:

$$a_j=\frac{v}{2}\frac{\Delta t}{h}-\frac{\sigma}{2}\frac{\Delta t}{h^{2}},\qquad b_j=1-\frac{v}{2}\frac{\Delta t}{h}+\sigma\frac{\Delta t}{h^{2}},\qquad c_j=-\frac{\sigma}{2}\frac{\Delta t}{h^{2}},$$ $$\xi_j^{n}=u_j^{n}+\frac{\sigma}{2}\frac{\Delta t}{h^{2}}\bigl(u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}\bigr)-\frac{v}{2}\frac{\Delta t}{h}\bigl(u_{j+1}^{n}-u_j^{n}\bigr)+\Delta t\,f(x_j).$$

Подставляя $v=-1$, $\sigma=0{,}3$, $f(x_j)=-3x_j^{2}$:

$$a_j=-\frac{\Delta t}{2h}-\frac{0{,}15\,\Delta t}{h^{2}},\qquad b_j=1+\frac{\Delta t}{2h}+\frac{0{,}3\,\Delta t}{h^{2}},\qquad c_j=-\frac{0{,}15\,\Delta t}{h^{2}},$$ $$\xi_j^{n}=u_j^{n}+\frac{0{,}15\,\Delta t}{h^{2}}\bigl(u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}\bigr)+\frac{\Delta t}{2h}\bigl(u_{j+1}^{n}-u_j^{n}\bigr)-3\,\Delta t\,x_j^{2}.$$

5. Проверка сходимости прогонки

Достаточное условие $|b_j|\ge|a_j|+|c_j|$. При $v=-1<0$ все три коэффициента имеют однозначные знаки: $a_j<0$, $c_j<0$, $b_j>0$; $|a_j|=\dfrac{\Delta t}{2h}+\dfrac{0{,}15\,\Delta t}{h^{2}}$, $|c_j|=\dfrac{0{,}15\,\Delta t}{h^{2}}$:

$$|a_j|+|c_j|=\frac{\Delta t}{2h}+\frac{0{,}3\,\Delta t}{h^{2}}<1+\frac{\Delta t}{2h}+\frac{0{,}3\,\Delta t}{h^{2}}=|b_j|.$$

То есть $|b_j|-\bigl(|a_j|+|c_j|\bigr)=1>0$ — строгое преобладание при любых $\Delta t,h>0$ благодаря единице в $b_j$ (она появилась именно из-за фиктивной производной по времени; в стационарном случае без неё было бы равенство). Прогонка сходится и устойчива.

6. Итерационное (прогоночное) соотношение

Полагаем $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$ и подставляем $u_{j-1}^{n+1}=\alpha_{j-1}u_j^{n+1}+\beta_{j-1}$. Прямым ходом слева направо:

$$\alpha_j=-\frac{a_j}{b_j+c_j\alpha_{j-1}},\qquad \beta_j=\frac{\xi_j^{n}-c_j\beta_{j-1}}{b_j+c_j\alpha_{j-1}}.$$

Обратным ходом справа налево находим $u_j^{n+1}$.

7. Коэффициенты $\alpha_1,\beta_1$ (левое ГУ, 2-й род)

Левое граничное условие $\dfrac{du}{dx}(x{=}0)=0$ аппроксимируем правой разностью на границе:

$$\frac{u_2^{n+1}-u_1^{n+1}}{h}=0\ \Longrightarrow\ u_1^{n+1}=u_2^{n+1}.$$

Сравнивая с $u_1^{n+1}=\alpha_1 u_2^{n+1}+\beta_1$:

$$\boxed{\alpha_1=1,\qquad \beta_1=0.}$$

8. Решение на правой границе $u_N$ (правое ГУ, 2-й род)

Правое условие $\dfrac{du}{dx}(x{=}1)=1$ аппроксимируем левой разностью на границе: $\dfrac{u_N^{n+1}-u_{N-1}^{n+1}}{h}=1$. Подставляя $u_{N-1}^{n+1}=\alpha_{N-1}u_N^{n+1}+\beta_{N-1}$:

$$\frac{u_N^{n+1}-(\alpha_{N-1}u_N^{n+1}+\beta_{N-1})}{h}=1\ \Longrightarrow\ u_N^{n+1}=\frac{h+\beta_{N-1}}{1-\alpha_{N-1}}.$$

9. Начальное приближение и условие окончания итераций

Нулевая итерация (начальное условие фиктивной нестационарной задачи) — берём в качестве удобного приближения свободный член; в методе установления годится любое разумное начальное приближение, так как решение сходится к стационару независимо от $u^0$:

$$u_j^{0}=f(x_j)=-3x_j^{2}.$$

Итерации продолжают, пока норма разности соседних приближений не станет меньше заданной точности $\varepsilon$:

$$\bigl\|u^{n+1}-u^{n}\bigr\|=\sqrt{h\sum_{j=1}^{N}\bigl(u_j^{n+1}-u_j^{n}\bigr)^{2}}\le\varepsilon.$$

10. Алгоритм

1. Задать начальное приближение $u_j^{0}=-3x_j^{2}$, $j=1,\dots,N$.
2. На каждой итерации $n$: из левого ГУ взять $\alpha_1=1$, $\beta_1=0$; прямым ходом $j=2,\dots,N-1$ вычислить $a_j,b_j,c_j,\xi_j^{n}$ и прогоночные коэффициенты $\alpha_j,\beta_j$.
3. Из правого ГУ найти $u_N^{n+1}=\dfrac{h+\beta_{N-1}}{1-\alpha_{N-1}}$.
4. Обратным ходом $j=N-1,\dots,1$ найти $u_j^{n+1}=\alpha_j u_{j+1}^{n+1}+\beta_j$.
5. Проверить условие окончания; если не выполнено — следующая итерация $n{+}1$.

1. Канонический вид

Приводим уравнение к виду $\dfrac{\partial u}{\partial t}+v_2\dfrac{\partial u}{\partial y}=\sigma_x\dfrac{\partial^2 u}{\partial x^2}+\sigma_y\dfrac{\partial^2 u}{\partial y^2}+\sigma_z\dfrac{\partial^2 u}{\partial z^2}+f$, замораживая переменные коэффициенты в узле $(j,k,m)$:

$$\sigma_x=7t\,y_k,\quad \sigma_y=5t\,z_m,\quad \sigma_z=6t\,x_j,\qquad v_2=8>0,\qquad f_{j,k,m}=-3\big(u_{j,k,m}\big)^2.$$

В источник $f$ входит только подлинная нелинейная реакция $-3u^2$. Конвективный член $8\,\partial u/\partial y$ (первая производная по $y$) НЕ переносится в $f$ и не берётся явно: согласно лекции 9, п. 10, первые производные учитываются неявно внутри соответствующей направленной подсхемы предиктора (здесь — подсхемы по $y$) с использованием противопотоковой (upwind) разности. Так как $v_2=8>0$, берётся левая (назад) разность $\dfrac{u_{j,k,m}-u_{j,k-1,m}}{h_y}$.

Операторы вторых производных: $\Lambda_{xx}u_{j,k,m}=\dfrac{u_{j+1,k,m}-2u_{j,k,m}+u_{j-1,k,m}}{h_x^2}$ и аналогично $\Lambda_{yy},\Lambda_{zz}$.

2. Расщепление времени (лекция 9, п. 8.1)

Интервал $[t^n,t^{n+1}]$ делится пополам точкой $t^{n+1/2}$; первая половина $\Delta t/2$ делится на три равные части точками $t^{n+1/6},t^{n+1/3}$. Предиктор (три неявные подсхемы) доводит решение до слоя $n+1/2$; корректор одной явной подсхемой продвигает от $t^n$ до $t^{n+1}$ на полном $\Delta t$. Итого четыре подсхемы; за шаг решение продвигается на $\Delta t$.

3. Предиктор — три неявные подсхемы

(1) неявно по $x$ (чистая диффузия): $\dfrac{u_{j,k,m}^{n+1/6}-u_{j,k,m}^{n}}{\Delta t/2}=\sigma_x\Lambda_{xx}u_{j,k,m}^{n+1/6}$.

(2) неявно по $y$ (диффузия + конвекция, upwind, $v_2>0$ — левая разность): $$\dfrac{u_{j,k,m}^{n+1/3}-u_{j,k,m}^{n+1/6}}{\Delta t/2}+8\,\dfrac{u_{j,k,m}^{n+1/3}-u_{j,k-1,m}^{n+1/3}}{h_y}=\sigma_y\Lambda_{yy}u_{j,k,m}^{n+1/3}.$$

(3) неявно по $z$ (чистая диффузия): $\dfrac{u_{j,k,m}^{n+1/2}-u_{j,k,m}^{n+1/3}}{\Delta t/2}=\sigma_z\Lambda_{zz}u_{j,k,m}^{n+1/2}$.

В каждой подсхеме неявны только производные по одной координате, поэтому каждая решается прогонкой. Источник $f=-3u^2$ в предикторе не участвует (он включён в корректор).

4. Прогоночный вид и сходимость (лекция 9, п. 8.2, 10.1)

Подсхемы (1) и (3) — диффузионные. Умножив (1) на $\Delta t/2$: $a_j=c_j=-\dfrac{\sigma_x}{2}\dfrac{\Delta t}{h_x^2}$, $b_j=1+\sigma_x\dfrac{\Delta t}{h_x^2}$; для (3) аналогично с $\sigma_z,h_z$. Здесь $a=c$ (симметрия). $|a|+|c|=\sigma\dfrac{\Delta t}{h^2}<1+\sigma\dfrac{\Delta t}{h^2}=|b|$.

Подсхема (2) — с конвекцией, коэффициенты НЕсимметричны. Умножив (2) на $\Delta t/2$ и собрав по узлам $y$:

$$\widetilde a_k=-\,8\,\frac{\Delta t}{2h_y}-\frac{\sigma_y}{2}\frac{\Delta t}{h_y^2}\ \ (u_{k-1}),\qquad \widetilde c_k=-\frac{\sigma_y}{2}\frac{\Delta t}{h_y^2}\ \ (u_{k+1}),\qquad \widetilde b_k=1+8\,\frac{\Delta t}{2h_y}+\sigma_y\frac{\Delta t}{h_y^2}.$$

Все коэффициенты $\widetilde a_k,\widetilde c_k<0$, $\widetilde b_k>0$, причём $|\widetilde a_k|+|\widetilde c_k|=8\dfrac{\Delta t}{2h_y}+\sigma_y\dfrac{\Delta t}{h_y^2}<\widetilde b_k$ (разница ровно $1$). Диагональное преобладание выполнено, прогонка сходится. Так как конвекция учтена неявно с правильным upwind-знаком, все три подсхемы предиктора абсолютно устойчивы (никаких ограничений на $\Delta t$ не возникает).

5. Прогоночные коэффициенты, $\alpha_1,\beta_1$

$u_{j,k,m}=\alpha_{j+1}u_{j+1,k,m}+\beta_{j+1}$, $\alpha_{j+1}=\dfrac{-c_j}{b_j+a_j\alpha_j}$, $\beta_{j+1}=\dfrac{\xi_j-a_j\beta_j}{b_j+a_j\alpha_j}$, где для (1) $\xi_j=u_{j,k,m}^{n}$. ГУ по $x$ 1-го рода, поэтому $\alpha_1=0$, $\beta_1=u_{0,k,m}=t^{\,n+1/6}\,y_k z_m$ (значение левого ГУ $u(t,x=0,y,z)=tyz$ на слое предиктора (1)). Для подсхем по $y$ и $z$ — аналогично со своими ГУ.

6. Корректор (лекция 9, формулы 9.20–9.21)

Аппроксимирует всё уравнение, центрируясь по времени относительно $t^{n+1/2}$; все операторы и источник берутся на известном слое $n+1/2$:

$$\frac{u_{j,k,m}^{n+1}-u_{j,k,m}^{n}}{\Delta t}=\sigma_x\Lambda_{xx}u^{n+1/2}+\sigma_y\Lambda_{yy}u^{n+1/2}+\sigma_z\Lambda_{zz}u^{n+1/2}-8\,\frac{u_{j,k,m}^{n+1/2}-u_{j,k-1,m}^{n+1/2}}{h_y}+f_{j,k,m}^{n+1/2}.$$

Все величины справа известны (слой $n+1/2$), поэтому корректор решается явно, рекуррентным соотношением, без прогонки; три направления, конвекция (upwind) и нелинейность учтены одновременно и согласованно.

7. Рекуррентное соотношение корректора

$$u_{j,k,m}^{n+1}=u_{j,k,m}^{n}+\Delta t\Big(\sigma_x\Lambda_{xx}u^{n+1/2}+\sigma_y\Lambda_{yy}u^{n+1/2}+\sigma_z\Lambda_{zz}u^{n+1/2}-8\,\frac{u_{j,k,m}^{n+1/2}-u_{j,k-1,m}^{n+1/2}}{h_y}+f_{j,k,m}^{n+1/2}\Big),$$

где $f_{j,k,m}^{n+1/2}=-3\big(u_{j,k,m}^{n+1/2}\big)^2$. База — слой $u^n$, правая часть отнесена к $t^{n+1/2}$: разность по времени центральная, что даёт 2-й порядок по времени.

8. Порядок и устойчивость

Предиктор сам даёт $O(\Delta t)$ по времени; центрированный корректор повышает временной порядок до $\Delta t^2$. По пространству: по $x$ и $z$ — вторые центральные разности $h_x^2,h_z^2$; по $y$ присутствует первая производная, аппроксимированная upwind (1-й порядок), поэтому по $y$ — только $h_y$. Итог:

$$O\big(\Delta t^2,\ h_x^2,\ h_y,\ h_z^2\big).$$

Конвекция учтена неявно с верным знаком, реакция вынесена в явный корректор, поэтому схема абсолютно устойчива (условия типа $\Delta t\le h_y/8$ не требуется).

9. Итог

• Предиктор = три неявные подсхемы (по $x,y,z$), каждая прогонка, даёт $u^{n+1/2}$; конвекция $8\,\partial u/\partial y$ входит НЕЯВНО в y-подсхему противопотоковой (левой) разностью.
• Корректор = одна явная подсхема, центрирована относительно $t^{n+1/2}$, учитывает все операторы, upwind-конвекцию и нелинейность $-3u^2$, даёт $u^{n+1}$.
• Диффузионные подсхемы: $a=c=-\dfrac{\sigma}{2}\dfrac{\Delta t}{h^2}$, $b=1+\sigma\dfrac{\Delta t}{h^2}$. Подсхема по $y$ — несимметрична: $\widetilde a_k=-8\dfrac{\Delta t}{2h_y}-\dfrac{\sigma_y}{2}\dfrac{\Delta t}{h_y^2}$, $\widetilde c_k=-\dfrac{\sigma_y}{2}\dfrac{\Delta t}{h_y^2}$, $\widetilde b_k=1+8\dfrac{\Delta t}{2h_y}+\sigma_y\dfrac{\Delta t}{h_y^2}$; во всех $|a|+|c|<|b|$. $\alpha_1=0$, $\beta_1$ — граничное значение.
• Порядок $O(\Delta t^2,\ h_x^2,\ h_y,\ h_z^2)$; схема абсолютно устойчива.

Уравнение и условия

Двумерное параболическое уравнение с нелинейной реакцией (Вариант 2, по условию требуется именно схема расщепления):

$$\frac{\partial u}{\partial t}=8t\,\frac{\partial^2 u}{\partial x^2}+5t\,\frac{\partial^2 u}{\partial y^2}-9u^2,\qquad u=u(t,x,y),\ \ x\in[0,1],\ y\in[0,1].$$ $$u(t,0,y)=ty,\quad u(t,1,y)=2ty,\qquad \frac{\partial u}{\partial y}(t,x,0)=0,\quad u(t,x,1)=1,\qquad u(0,x,y)=\sin(xy).$$

Коэффициенты диффузии $\sigma_x=8t$, $\sigma_y=5t$, реакционный член $g(u)=-9u^2$. По $x$ — оба граничных условия 1-го рода (Дирихле); по $y$ — слева условие 2-го рода (Нейман, $u_y=0$), справа 1-го рода (Дирихле). Сетка $u_{j,k}^{n}=u(t^n,x_j,y_k)$, $x_j=(j-1)h_x$, $y_k=(k-1)h_y$, $t^{n+1/2}=t^n+\Delta t/2$.

Разностные операторы:

$$\Lambda_{xx}u_{j,k}=\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^2},\qquad \Lambda_{yy}u_{j,k}=\frac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^2}.$$

Идея схемы расщепления (дробных шагов)

В отличие от схемы переменных направлений (где в каждой подсхеме присутствуют оба оператора, один неявно, другой явно, с множителем $\tfrac{\sigma}{2}$), схема расщепления по физическим процессам на каждом дробном шаге оставляет только один пространственный оператор — целиком, с полным коэффициентом $\sigma$. Шаг $\Delta t$ делится точкой $t^{n+1/2}$:

• Полушаг ① ($n\to n+1/2$): «работает» только диффузия по $x$ — неявно по $x$ (прогонка по индексу $j$); производных по $y$ нет.
• Полушаг ② ($n+1/2\to n+1$): «работает» только диффузия по $y$ — неявно по $y$ (прогонка по индексу $k$); производных по $x$ нет.

Нелинейный реакционный член $-9u^2$ линеаризуем (метод запаздывающего коэффициента): $-9u^2\approx -9\,u^{(\text{стар})}u^{(\text{нов})}$, и распределяем поровну между полушагами (по $\Delta t/2$). На полушаге ① берём $-9\,u_{j,k}^{n}\,u_{j,k}^{n+1/2}$, на полушаге ② — $-9\,u_{j,k}^{n+1/2}\,u_{j,k}^{n+1}$. Так трёхдиагональность сохраняется, а коэффициент при неизвестном (запаздывающий множитель) неотрицателен при $u\ge0$.

Подсхема ① ($n\to n+1/2$): неявно по $x$

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t/2}=8t^{n+1/2}\,\Lambda_{xx}u_{j,k}^{n+1/2}-9\,u_{j,k}^{n}\,u_{j,k}^{n+1/2}.$$

$y$-оператора нет — это суть расщепления. Умножив на $\Delta t/2$ и собрав неизвестные слоя $n+1/2$ слева, получаем трёхдиагональную по $j$ систему

$$a_j\,u_{j-1,k}^{n+1/2}+b_j\,u_{j,k}^{n+1/2}+c_j\,u_{j+1,k}^{n+1/2}=\xi_{j,k},$$ $$a_j=c_j=-\frac{8t^{n+1/2}\Delta t}{2h_x^2}=-\frac{4t^{n+1/2}\Delta t}{h_x^2},\qquad b_j=1+\frac{8t^{n+1/2}\Delta t}{h_x^2}+\frac{9\,u_{j,k}^{n}\,\Delta t}{2},\qquad \xi_{j,k}=u_{j,k}^{n}.$$

Сходимость прогонки (диагональное преобладание ①)

$$|a_j|+|c_j|=\frac{8t^{n+1/2}\Delta t}{h_x^2},\qquad |b_j|-\bigl(|a_j|+|c_j|\bigr)=1+\frac{9\,u_{j,k}^{n}\,\Delta t}{2}.$$

При $t\ge 0$ уже диффузионная часть даёт нестрогое равенство, а свободный «$+1$» обеспечивает строгое преобладание; член $\tfrac{9}{2}u^{n}\Delta t\ge0$ (решение неотрицательно: $u^0=\sin(xy)\ge0$ на $[0,1]^2$, ГУ $ty,2ty,1\ge0$ при $t\ge0$) лишь дополнительно усиливает диагональ. Условие $|b_j|>|a_j|+|c_j|$ выполнено безусловно $\Rightarrow$ прогонка устойчива.

Рекуррентное прогоночное соотношение (по $x$)

$$u_{j,k}^{n+1/2}=\alpha_{j+1}\,u_{j+1,k}^{n+1/2}+\beta_{j+1},\qquad \alpha_{j+1}=\frac{-c_j}{b_j+a_j\alpha_j},\quad \beta_{j+1}=\frac{\xi_{j,k}-a_j\beta_j}{b_j+a_j\alpha_j}.$$

Старт прогонки $\alpha_1,\beta_1$ из ГУ по $x$ (оба 1-го рода)

Левая граница $x=0$ (узел $j=1$) — условие Дирихле $u_{1,k}^{n+1/2}=t^{n+1/2}y_k$. Записав его в форме $u_{1,k}=\alpha_1 u_{2,k}+\beta_1$, получаем

$$\boxed{\;\alpha_1=0,\qquad \beta_1=t^{n+1/2}\,y_k.\;}$$

Правая граница $x=1$ (узел $j=N_x$) — тоже Дирихле, берётся напрямую: $u_{N_x,k}^{n+1/2}=2\,t^{n+1/2}\,y_k$. Дальше обратный ход прогонки от $j=N_x-1$ к $j=2$.

Подсхема ② ($n+1/2\to n+1$): неявно по $y$

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t/2}=5t^{n+1}\,\Lambda_{yy}u_{j,k}^{n+1}-9\,u_{j,k}^{n+1/2}\,u_{j,k}^{n+1}.$$

Теперь $x$-оператора нет. Умножив на $\Delta t/2$, неизвестные слоя $n+1$ — слева. Трёхдиагональная по $k$ система

$$\tilde a_k\,u_{j,k-1}^{n+1}+\tilde b_k\,u_{j,k}^{n+1}+\tilde c_k\,u_{j,k+1}^{n+1}=\tilde\xi_{j,k},$$ $$\tilde a_k=\tilde c_k=-\frac{5t^{n+1}\Delta t}{2h_y^2},\qquad \tilde b_k=1+\frac{5t^{n+1}\Delta t}{h_y^2}+\frac{9\,u_{j,k}^{n+1/2}\,\Delta t}{2},\qquad \tilde\xi_{j,k}=u_{j,k}^{n+1/2}.$$

Сходимость прогонки (диагональное преобладание ②)

$$|\tilde a_k|+|\tilde c_k|=\frac{5t^{n+1}\Delta t}{h_y^2},\qquad |\tilde b_k|-\bigl(|\tilde a_k|+|\tilde c_k|\bigr)=1+\frac{9\,u_{j,k}^{n+1/2}\Delta t}{2}>0.$$

Выполнено безусловно при $t\ge0$, $u^{n+1/2}\ge0$.

Рекуррентное прогоночное соотношение (по $y$)

$$u_{j,k}^{n+1}=\hat\alpha_{k+1}\,u_{j,k+1}^{n+1}+\hat\beta_{k+1},\qquad \hat\alpha_{k+1}=\frac{-\tilde c_k}{\tilde b_k+\tilde a_k\hat\alpha_k},\quad \hat\beta_{k+1}=\frac{\tilde\xi_{j,k}-\tilde a_k\hat\beta_k}{\tilde b_k+\tilde a_k\hat\alpha_k}.$$

Старт прогонки $\hat\alpha_1,\hat\beta_1$ из ГУ по $y$ (слева 2-й род, справа 1-й род)

Левая граница $y=0$ (узел $k=1$) — условие Неймана $\dfrac{\partial u}{\partial y}(t,x,0)=0$. Аппроксимируем правой разностью первого порядка:

$$\frac{u_{j,2}^{n+1}-u_{j,1}^{n+1}}{h_y}=0\ \Longrightarrow\ u_{j,1}^{n+1}=u_{j,2}^{n+1}.$$

Сравнивая с $u_{j,1}^{n+1}=\hat\alpha_1 u_{j,2}^{n+1}+\hat\beta_1$, получаем

$$\boxed{\;\hat\alpha_1=1,\qquad \hat\beta_1=0.\;}$$

Правая граница $y=1$ (узел $k=N_y$) — условие Дирихле $u_{j,N_y}^{n+1}=1$, берётся напрямую; затем обратный ход прогонки от $k=N_y-1$ к $k=2$.

Порядок аппроксимации и алгоритм

Схема расщепления (дробных шагов) имеет первый порядок по времени и второй по пространству во внутренних узлах, $O(\Delta t,\,h_x^2,\,h_y^2)$, и абсолютно устойчива (каждый дробный шаг — неявная одномерная схема). Алгоритм перехода $n\to n+1$:

1. Начальное поле $u_{j,k}^{0}=\sin(x_j y_k)$.
2. Полушаг ①: для каждой строки $k$ — прогонка по $x$ с коэффициентами $a_j=c_j,\,b_j,\,\xi_{j,k}$; старт $\alpha_1=0,\ \beta_1=t^{n+1/2}y_k$, справа $u_{N_x,k}^{n+1/2}=2t^{n+1/2}y_k$ $\Rightarrow$ получаем $u^{n+1/2}$.
3. Полушаг ②: для каждого столбца $j$ — прогонка по $y$ с $\tilde a_k=\tilde c_k,\,\tilde b_k,\,\tilde\xi_{j,k}$; старт $\hat\alpha_1=1,\ \hat\beta_1=0$ (Нейман), справа $u_{j,N_y}^{n+1}=1$ (Дирихле) $\Rightarrow$ получаем $u^{n+1}$.
4. Повторять по $n$ до $t=t_k$.

Шаг 1. Тип уравнения и канонический вид

Перенесём конвективный член по \(x\) и источник так, чтобы слева остался только \(u_t\):

$$\dfrac{\partial u}{\partial t} + 3\,\dfrac{\partial u}{\partial x} = t - 2\,\dfrac{\partial u}{\partial y}\ \Longrightarrow\ \dfrac{\partial u}{\partial t} = -3\,\dfrac{\partial u}{\partial x} - 2\,\dfrac{\partial u}{\partial y} + t.$$

Вторых производных по координатам нет — это уравнение переноса (конвекции) первого порядка с двумя пространственными переменными. Канонические параметры:

$$v_1 = 3\ (\text{скорость по }x),\qquad v_2 = 2\ (\text{скорость по }y),\qquad f(t,x,y) = t.$$

Шаг 2. Сетка, выбор разностей и граничных условий

Узлы: \(x_j=(j-1)h_x,\ y_k=(k-1)h_y,\ t^n=n\,\Delta t\). Обе скорости положительны (\(v_1>0,\ v_2>0\)), характеристики идут в сторону роста \(x\) и \(y\), информация приходит «слева». Поэтому для устойчивости берём левые (против потока) разности:

$$\Lambda_x u_{jk} = \dfrac{u_{jk}-u_{j-1,k}}{h_x},\qquad \Lambda_y u_{jk} = \dfrac{u_{jk}-u_{j,k-1}}{h_y}.$$

Входными (используемыми) являются границы со стороны притока — \(x=0\) и \(y=0\):

$$u(t,\,x=0,\,y)=1,\qquad u(t,\,x,\,y=0)=1.$$

Граничные условия на выходных границах \(u(t,\,x=1,\,y)=t\) и \(u(t,\,x,\,y=1)=t\) для уравнения переноса 1-го порядка в рекуррентные соотношения не входят (поток их «выносит»). Начальное условие:

$$u^{0}_{jk}=x_j\,y_k=(j-1)h_x\,(k-1)h_y.$$

Шаг 3. Схема переменных направлений (Писмена–Рэкфорда)

По канону курса (семинар С9) каждая подсхема продвигает решение через промежуточный слой \(t^{n+1/2}\) со знаменателем \(\Delta t\), при этом коэффициент при каждом операторе берётся в половинном размере (по \(v_i/2\)), а источник учитывается один раз в подсхеме ② на слое \(t^{n+1/2}=(n+\tfrac12)\Delta t\). В подсхеме ① оператор по \(x\) — неявный (на \(n+\tfrac12\)), по \(y\) — явный (на \(n\)); в подсхеме ② наоборот.

Подсхема ① — неявная по \(x\), явная по \(y\):

$$\dfrac{u^{\,n+1/2}_{jk}-u^{\,n}_{jk}}{\Delta t} = -\dfrac{3}{2}\,\dfrac{u^{\,n+1/2}_{jk}-u^{\,n+1/2}_{j-1,k}}{h_x} \;-\; \dfrac{u^{\,n}_{jk}-u^{\,n}_{j,k-1}}{h_y}.$$

(множитель при операторе по \(y\) равен \(v_2/2=1\)).

Подсхема ② — неявная по \(y\), явная по \(x\):

$$\dfrac{u^{\,n+1}_{jk}-u^{\,n+1/2}_{jk}}{\Delta t} = -\dfrac{3}{2}\,\dfrac{u^{\,n+1/2}_{jk}-u^{\,n+1/2}_{j-1,k}}{h_x} \;-\; \dfrac{u^{\,n+1}_{jk}-u^{\,n+1}_{j,k-1}}{h_y} \;+\; t^{\,n+1/2}.$$

Шаг 4. Проверка согласованности (аппроксимации)

Сложим обе подсхемы (промежуточный слой \(u^{\,n+1/2}\) сокращается):

$$\dfrac{u^{\,n+1}_{jk}-u^{\,n}_{jk}}{\Delta t} = -3\,\Lambda_x u^{\,n+1/2}_{jk} \;-\; 2\cdot\dfrac{\Lambda_y u^{\,n}_{jk}+\Lambda_y u^{\,n+1}_{jk}}{2} \;+\; t^{\,n+1/2}.$$

Коэффициент при операторе по \(x\) восстановился как \(\tfrac32+\tfrac32=3\); оператор по \(y\) вошёл как полусумма слоёв \(n\) и \(n{+}1\) с полным коэффициентом \(2\); источник вошёл один раз с полным шагом \(\Delta t\). Схема аппроксимирует исходное уравнение \(u_t=-3u_x-2u_y+t\). Правая часть центрирована относительно \(t^{n+1/2}\), поэтому по времени достигается второй порядок:

$$O\!\left(\Delta t^{2},\ h_x,\ h_y\right).$$

Шаг 5. Рекуррентные соотношения для подсхем

Подсхема ①. Введём \(r_x=\dfrac{3}{2}\dfrac{\Delta t}{h_x}\). Соберём неявные члены (слой \(n+\tfrac12\)) слева:

$$(1+r_x)\,u^{\,n+1/2}_{jk} - r_x\,u^{\,n+1/2}_{j-1,k} = u^{\,n}_{jk} - \dfrac{\Delta t}{h_y}\big(u^{\,n}_{jk}-u^{\,n}_{j,k-1}\big).$$

Из-за односторонней (левой) разности связь по \(x\) двухдиагональная (нижнетреугольная), поэтому прогонка не требуется — значения находятся прямой подстановкой слева направо от входной границы \(x=0\):

$$u^{\,n+1/2}_{1,k}=1,\qquad u^{\,n+1/2}_{jk}=\dfrac{r_x\,u^{\,n+1/2}_{j-1,k} + u^{\,n}_{jk} - \dfrac{\Delta t}{h_y}\big(u^{\,n}_{jk}-u^{\,n}_{j,k-1}\big)}{1+r_x},\quad j=2,\dots,N_x.$$

Подсхема ②. Введём \(r_y=\dfrac{\Delta t}{h_y}\) (коэффициент по \(y\) равен \(v_2/2=1\)) и \(t^{\,n+1/2}=(n+\tfrac12)\Delta t\). Соберём неявные члены (слой \(n+1\)) слева:

$$(1+r_y)\,u^{\,n+1}_{jk} - r_y\,u^{\,n+1}_{j,k-1} = u^{\,n+1/2}_{jk} - \dfrac{3}{2}\dfrac{\Delta t}{h_x}\big(u^{\,n+1/2}_{jk}-u^{\,n+1/2}_{j-1,k}\big) + \Delta t\,t^{\,n+1/2}.$$

Связь по \(y\) также двухдиагональная — прямая подстановка снизу вверх от входной границы \(y=0\), где по условию \(u(t,x,y=0)=1\):

$$u^{\,n+1}_{j,1}=1,\qquad u^{\,n+1}_{jk}=\dfrac{r_y\,u^{\,n+1}_{j,k-1} + u^{\,n+1/2}_{jk} - \dfrac{3}{2}\dfrac{\Delta t}{h_x}\big(u^{\,n+1/2}_{jk}-u^{\,n+1/2}_{j-1,k}\big) + \Delta t\,t^{\,n+1/2}}{1+r_y},\quad k=2,\dots,N_y.$$

Шаг 6. Устойчивость

В обеих подсхемах диагональный коэффициент по модулю превосходит сумму внедиагональных: \(|1+r_x|>|r_x|\) и \(|1+r_y|>|r_y|\) (диагональное преобладание), что гарантирует корректность прямой подстановки. Неявная по направлениям запись (схема переменных направлений) безусловно устойчива; ограничение на шаг типа КФУ (например \(r_x+r_y\le 1\)) потребовалось бы лишь для явной схемы.

Итог

• Уравнение — перенос 1-го порядка; канонические \(v_1=3,\ v_2=2,\ f=t\).
• Используются входные ГУ \(u(t,0,y)=1,\ u(t,x,0)=1\); выходные ГУ не используются; НУ \(u^0_{jk}=x_jy_k\).
• Схема переменных направлений с половинными коэффициентами при операторах, источник один раз на \(t^{n+1/2}\); порядок \(O(\Delta t^2,h_x,h_y)\).
• Подсхема ① — неявная по \(x\), подсхема ② — неявная по \(y\); связь двухдиагональная, решение прямой подстановкой от входной границы (\(u^{n+1/2}_{1,k}=1\), \(u^{n+1}_{j,1}=1\)).

Условие

$$0{,}2\left(\dfrac{\partial^2 u}{\partial x^2}+\dfrac{\partial^2 u}{\partial y^2}\right)+3e^{xy}=0,\qquad u=u(x,y),\ \ 0\le x\le1,\ 0\le y\le1,$$

$$\begin{cases}u(0,y)=y\\ u(1,y)=2y\end{cases}\qquad\begin{cases}u(x,0)=0\\ u(x,1)=x\end{cases}$$

Вторые производные по обеим координатам входят с одинаковым (положительным) знаком — это эллиптическое уравнение (уравнение Пуассона). Граничные условия — 1-го рода (значения функции на всех четырёх сторонах квадрата). По формулировке варианта требуется применить метод простой итерации и записать: шаг итерации, итерационное соотношение, условие окончания и начальное приближение.

Замечание о данных. Угловые значения, задаваемые двумя соседними сторонами, несовместны: в углу $(0,1)$ из $u(0,y)=y$ следует $1$, а из $u(x,1)=x$ — $0$; в углу $(1,1)$ из $u(1,y)=2y$ следует $2$, а из $u(x,1)=x$ — $1$. Это особенность исходных данных; на построение и сходимость итераций во внутренних узлах она не влияет (углы в пятиточечный шаблон не входят).

1. Приведение к каноническому виду

Делим уравнение на $0{,}2$ (коэффициент $3/0{,}2=15$), вторые производные оставляем слева, источник — справа:

$$\dfrac{\partial^2 u}{\partial x^2}+\dfrac{\partial^2 u}{\partial y^2}=-15\,e^{xy},\qquad F(x,y)=15\,e^{xy}>0.$$

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

2. Введение фиктивного времени (итераций)

Добавим производную по фиктивному времени:

$$\dfrac{\partial\widetilde u}{\partial t}=\dfrac{\partial^2 \widetilde u}{\partial x^2}+\dfrac{\partial^2 \widetilde u}{\partial y^2}+15\,e^{xy}.$$

Граничные условия и источник от времени не зависят, поэтому при $t\to\infty$ имеем $\partial\widetilde u/\partial t\to0$ и $\widetilde u(x,y,t)\to u(x,y)$ — решение исходной стационарной задачи. Номер слоя по времени $s$ есть номер итерации.

3. Сетка и аппроксимация

Равномерная сетка $x_i=i\,h$, $y_k=k\,h$, $i,k=0,\dots,N$, $h=1/N$ ($h_x=h_y=h$). Внутренние узлы $i,k=1,\dots,N-1$. Вторые производные — центральными разностями:

$$\Lambda_x u=\dfrac{u_{i+1,k}-2u_{i,k}+u_{i-1,k}}{h^2},\qquad \Lambda_y u=\dfrac{u_{i,k+1}-2u_{i,k}+u_{i,k-1}}{h^2}.$$

4. Шаг итерации (явная схема)

Производную по времени берём вперёд, пространственные операторы — на известном $s$-м слое (шаг $\tau$):

$$\dfrac{u_{i,k}^{\,s+1}-u_{i,k}^{\,s}}{\tau}=\Lambda_x u^{\,s}+\Lambda_y u^{\,s}+15\,e^{x_iy_k}.$$

Шаг $\tau$ выбирается из условия устойчивости явной схемы для двумерного параболического уравнения:

$$\tau\le\dfrac{h^2}{4}\qquad\Big(\text{в общем случае }\tau\le\dfrac{1}{2}\dfrac{1}{1/h_x^2+1/h_y^2}\Big),$$

причём $\tau=h^2/4$ — оптимальный (наибольшая скорость сходимости простой итерации). Число итераций до установления $s\sim 1/h^2$.

5. Итерационное соотношение

Выражая новое приближение явно, получаем рабочую формулу метода простой итерации:

$$\boxed{\,u_{i,k}^{\,s+1}=u_{i,k}^{\,s}+\tau\left[\dfrac{u_{i+1,k}^{\,s}+u_{i-1,k}^{\,s}+u_{i,k+1}^{\,s}+u_{i,k-1}^{\,s}-4u_{i,k}^{\,s}}{h^2}+15\,e^{x_iy_k}\right],\qquad i,k=1,\dots,N-1.}$$

При оптимальном $\tau=h^2/4$ член $u_{i,k}^{\,s}$ сокращается, и формула переходит в усреднение по четырём соседям (классическая схема Якоби):

$$u_{i,k}^{\,s+1}=\dfrac{1}{4}\left(u_{i+1,k}^{\,s}+u_{i-1,k}^{\,s}+u_{i,k+1}^{\,s}+u_{i,k-1}^{\,s}\right)+\dfrac{h^2}{4}\cdot15\,e^{x_iy_k}.$$

6. Граничные узлы

Метод простой итерации (схема Якоби) не использует прогонку: значения во внутренних узлах пересчитываются по явной формуле п.5. Поэтому прогоночные коэффициенты $\alpha_1,\beta_1$ для этого варианта не требуются (они появляются лишь в вариантах с установлением по неявным схемам / схеме переменных направлений). Граничные значения 1-го рода задаются напрямую и фиксированы на всех итерациях:

$$u_{0,k}=y_k,\quad u_{N,k}=2y_k,\quad u_{i,0}=0,\quad u_{i,N}=x_i.$$

7. Начальное приближение (нулевая итерация)

Метод установления сходится независимо от старта, поэтому во внутренних узлах берут произвольную (проще — нулевую) функцию:

$$u_{i,k}^{\,0}=0,\qquad i,k=1,\dots,N-1,$$

а в граничных узлах сразу заносят точные значения из ГУ: $u_{0,k}^{0}=y_k$, $u_{N,k}^{0}=2y_k$, $u_{i,0}^{0}=0$, $u_{i,N}^{0}=x_i$.

8. Условие окончания итерационного процесса

Итерации продолжают до установления — пока норма разности двух соседних приближений не станет меньше точности $\varepsilon$:

$$\max_{i,k}\big|u_{i,k}^{\,s+1}-u_{i,k}^{\,s}\big|<\varepsilon\qquad\Big(\text{или }\big\|u^{s+1}-u^{s}\big\|=\sqrt{h^2\!\!\sum_{i,k}\big(u_{i,k}^{\,s+1}-u_{i,k}^{\,s}\big)^2}<\varepsilon\Big).$$

9. Алгоритм

1. Задать сетку $h=1/N$, узлы $x_i=ih$, $y_k=kh$; выбрать $\tau\le h^2/4$ и точность $\varepsilon$.
2. Занести начальное приближение: $u_{i,k}^{0}=0$ внутри; на границе — значения из ГУ.
3. Для всех внутренних узлов $i,k=1,\dots,N-1$ вычислить $u_{i,k}^{\,s+1}$ по соотношению п.5.
4. Граничные значения оставить неизменными.
5. Если $\max_{i,k}|u_{i,k}^{\,s+1}-u_{i,k}^{\,s}|<\varepsilon$ — решение установилось; иначе $s:=s+1$ и вернуться к шагу 3.

Уравнение и условия (Вариант 2)

Дано трёхмерное по пространству ($x,y,z$) уравнение переноса с диффузией:

$$\frac{\partial u}{\partial t}+5\,\frac{\partial u}{\partial x}+7\,\frac{\partial u}{\partial y}=0{,}2\,\frac{\partial^2 u}{\partial y^2}+0{,}3\,\frac{\partial^2 u}{\partial z^2}-4t,$$

с граничными и начальным условиями

$$u(t,0,y,z)=yz,\quad \begin{cases}u(t,x,0,z)=xz,\\ u(t,x,1,z)=4xz,\end{cases}\quad \begin{cases}u(t,x,y,0)=xy,\\ u(t,x,y,1)=5xy,\end{cases}\quad u(0,x,y,z)=1.$$

Сеточная функция $u_{j,k,m}^{\,n}$ задана в узлах $x_j,\ y_k,\ z_m$ с шагами $h_x,h_y,h_z$; верхний индекс $n$ — временной слой, $t^{\,n}=n\,\Delta t$.

1. Анализ операторов уравнения

• Вторые производные (диффузия): по $y$ с коэффициентом $\sigma_y=0{,}2$ и по $z$ с коэффициентом $\sigma_z=0{,}3$. По $x$ второй производной нет.
• Первые производные (конвекция): $5\,\partial u/\partial x$ и $7\,\partial u/\partial y$ с положительными скоростями $v_x=5>0$, $v_y=7>0$.
• Свободный член: $f=-4t$ (зависит только от времени).

Таким образом, по $x$ — только конвекция, по $y$ — и конвекция, и диффузия, по $z$ — только диффузия.

Аппроксимация производных. Вторые производные — центральные разности второго порядка:

$$\Lambda_{yy}u_{j,k,m}=\frac{u_{j,k+1,m}-2u_{j,k,m}+u_{j,k-1,m}}{h_y^2},\qquad \Lambda_{zz}u_{j,k,m}=\frac{u_{j,k,m+1}-2u_{j,k,m}+u_{j,k,m-1}}{h_z^2}.$$

Первые производные стоят в левой части, поэтому при $v>0$ берётся левая (противопоточная) разность с левым граничным условием. Так как $v_x=5>0$, $v_y=7>0$:

$$\Lambda_x u_{j,k,m}=\frac{u_{j,k,m}-u_{j-1,k,m}}{h_x},\qquad \Lambda_y^{(1)}u_{j,k,m}=\frac{u_{j,k,m}-u_{j,k-1,m}}{h_y}.$$

Эти разности имеют первый порядок по соответствующей координате.

2. Идея схемы предиктор-корректор (по учебнику, гл.9 §8)

Интервал $\delta t=[t^{\,n},t^{\,n+1}]$ расщепляется пополам точкой $t^{\,n+1/2}=t^{\,n}+\Delta t/2$. Первый полуинтервал $\Delta t/2$ (предиктор) дополнительно делится на части (узлы $t^{\,n+1/6}, t^{\,n+1/3}$), на каждой из которых записывается неявная подсхема со знаменателем $\Delta t/2$, неявная только по одной координате (трёхдиагональная, решается прогонкой). Это обеспечивает абсолютную устойчивость. Корректор — одно соотношение на полный шаг $\Delta t$ от слоя $u^{\,n}$, правая часть которого вычислена в средней точке $t^{\,n+1/2}$ (через предиктор $\bar u$); это придаёт центральный характер разности по времени и поднимает порядок до $O(\Delta t^2)$. Всего схема состоит из четырёх подсхем: трёх подсхем предиктора и корректора.

Конвекция по координате внедряется в ту же неявную подсхему, что и диффузия по этой координате (как в семинаре 8, пример 2), а не в правую часть.

Предиктор (на интервале $\Delta t/2$)

① Подсхема по $x$ ($t^{\,n}\to t^{\,n+1/6}$, знаменатель $\Delta t/2$, неявно — конвекция по $x$):

$$\frac{u_{j,k,m}^{\,n+1/6}-u_{j,k,m}^{\,n}}{\Delta t/2}=-5\,\Lambda_x u_{j,k,m}^{\,n+1/6}.$$

② Подсхема по $y$ ($t^{\,n+1/6}\to t^{\,n+1/3}$, знаменатель $\Delta t/2$, неявно — диффузия и конвекция по $y$):

$$\frac{u_{j,k,m}^{\,n+1/3}-u_{j,k,m}^{\,n+1/6}}{\Delta t/2}=0{,}2\,\Lambda_{yy}u_{j,k,m}^{\,n+1/3}-7\,\Lambda_y^{(1)}u_{j,k,m}^{\,n+1/3}.$$

③ Подсхема по $z$ ($t^{\,n+1/3}\to t^{\,n+1/2}$, знаменатель $\Delta t/2$, неявно — диффузия по $z$, сюда же отнесём свободный член):

$$\frac{u_{j,k,m}^{\,n+1/2}-u_{j,k,m}^{\,n+1/3}}{\Delta t/2}=0{,}3\,\Lambda_{zz}u_{j,k,m}^{\,n+1/2}-4t^{\,n}.$$

Результат предиктора — промежуточные значения $\bar u_{j,k,m}\equiv u_{j,k,m}^{\,n+1/2}$ в средней точке $t^{\,n+1/2}$.

Корректор (полный шаг $\Delta t$ от $u^{\,n}$)

Все операторы берутся в средней точке $t^{\,n+1/2}$, то есть по предиктору $\bar u$:

$$\boxed{\ \frac{u_{j,k,m}^{\,n+1}-u_{j,k,m}^{\,n}}{\Delta t} =0{,}2\,\Lambda_{yy}\bar u_{j,k,m}+0{,}3\,\Lambda_{zz}\bar u_{j,k,m} -5\,\Lambda_x \bar u_{j,k,m}-7\,\Lambda_y^{(1)}\bar u_{j,k,m}-4t^{\,n+1/2}.\ }$$

Именно полный шаг $\Delta t$ от $u^{\,n}$ с правой частью в точке $t^{\,n+1/2}$ (а не «полушаг от $\bar u$») делает разность по времени центральной. Тейлор-проверка даёт погрешность аппроксимации (невязку разностного уравнения) корректора $$\psi_t=\frac{u^{\,n+1}-u^{\,n}}{\Delta t}-\big(Lu+f\big)^{\,n+1/2}=\frac{\Delta t^2}{24}\,u'''+\frac{\Delta t^3}{48}\,u^{(4)}=O(\Delta t^{2}),$$ то есть локальная (за шаг) погрешность $O(\Delta t^{3})$, а глобально — второй порядок по времени $O(\Delta t^2)$. (Для сравнения, вариант «$(u^{\,n+1}-\bar u)/(\Delta t/2)=\dots$» даёт невязку $O(\Delta t)$ и лишь первый порядок по времени — поэтому он отвергнут.)

3. Рекуррентные соотношения подсхем (прогонка)

Каждая неявная подсхема — трёхдиагональная система вида $a_i u_{i-1}+b_i u_i+c_i u_{i+1}=\xi_i$, решаемая прогонкой

$$\alpha_{i+1}=\frac{-c_i}{b_i+a_i\alpha_i},\qquad \beta_{i+1}=\frac{\xi_i-a_i\beta_i}{b_i+a_i\alpha_i},\qquad u_i=\alpha_{i+1}u_{i+1}+\beta_{i+1}.$$

Подсхема ① — прогонка по $j$ (направление $x$). Перенос конвективных членов налево даёт двухточечное (нижнетреугольное) рекуррентное соотношение «по потоку»:

$$a_j=-\frac{5\,\Delta t}{2h_x},\qquad b_j=1+\frac{5\,\Delta t}{2h_x},\qquad c_j=0,$$ $$\xi_{j,k,m}=u_{j,k,m}^{\,n},\qquad u_{j,k,m}^{\,n+1/6}=\frac{u_{j,k,m}^{\,n}+\dfrac{5\Delta t}{2h_x}\,u_{j-1,k,m}^{\,n+1/6}}{\,1+\dfrac{5\Delta t}{2h_x}\,}.$$

Стартовое значение даёт левое ГУ $u(t,0,y,z)=y_k z_m$ (так как $v_x>0$, считаем слева направо).

Подсхема ② — прогонка по $k$ (направление $y$). Неявны и диффузия $0{,}2\,\Lambda_{yy}$, и конвекция $7\,\Lambda_y^{(1)}$. С обозначением $r_y=\dfrac{0{,}2\,(\Delta t/2)}{h_y^2}=\dfrac{0{,}1\,\Delta t}{h_y^2}$:

$$a_k=-\frac{0{,}1\,\Delta t}{h_y^2}-\frac{7\,\Delta t}{2h_y},\qquad b_k=1+\frac{0{,}2\,\Delta t}{h_y^2}+\frac{7\,\Delta t}{2h_y},\qquad c_k=-\frac{0{,}1\,\Delta t}{h_y^2},$$ $$\xi_{j,k,m}=u_{j,k,m}^{\,n+1/6},\qquad u_{j,k,m}^{\,n+1/3}=\alpha_{k+1}\,u_{j,k+1,m}^{\,n+1/3}+\beta_{k+1}.$$

Подсхема ③ — прогонка по $m$ (направление $z$). Неявна диффузия $0{,}3\,\Lambda_{zz}$. С $r_z=\dfrac{0{,}3\,(\Delta t/2)}{h_z^2}=\dfrac{0{,}15\,\Delta t}{h_z^2}$:

$$a_m=c_m=-r_z=-\frac{0{,}15\,\Delta t}{h_z^2},\qquad b_m=1+2r_z=1+\frac{0{,}3\,\Delta t}{h_z^2},$$ $$\xi_{j,k,m}=u_{j,k,m}^{\,n+1/3}-\frac{\Delta t}{2}\cdot 4t^{\,n} =u_{j,k,m}^{\,n+1/3}-2\,\Delta t\,t^{\,n},\qquad u_{j,k,m}^{\,n+1/2}=\alpha_{m+1}\,u_{j,k,m+1}^{\,n+1/2}+\beta_{m+1}.$$

Корректор — рекуррентное (явное) соотношение. Корректор разрешён относительно $u^{\,n+1}$ непосредственно (правая часть полностью известна через $\bar u=u^{\,n+1/2}$):

$$u_{j,k,m}^{\,n+1}=u_{j,k,m}^{\,n}+\Delta t\Big(0{,}2\,\Lambda_{yy}\bar u_{j,k,m}+0{,}3\,\Lambda_{zz}\bar u_{j,k,m} -5\,\Lambda_x \bar u_{j,k,m}-7\,\Lambda_y^{(1)}\bar u_{j,k,m}-4t^{\,n+1/2}\Big).$$

4. Сходимость прогонок (диагональное преобладание)

• ① ($x$): $|a_j|+|c_j|=\dfrac{5\Delta t}{2h_x}<1+\dfrac{5\Delta t}{2h_x}=|b_j|$.
• ② ($y$): $|a_k|+|c_k|=\dfrac{0{,}2\,\Delta t}{h_y^2}+\dfrac{7\Delta t}{2h_y}<1+\dfrac{0{,}2\,\Delta t}{h_y^2}+\dfrac{7\Delta t}{2h_y}=|b_k|$ (разность $=1>0$).
• ③ ($z$): $|a_m|+|c_m|=2r_z<1+2r_z=|b_m|$.

Во всех подсхемах преобладание строгое — прогонка устойчива безусловно.

5. Граничные и начальное условия, старт прогонок

• По $x$ ($v_x=5>0$, левое ГУ): $u(t,0,y,z)=y_k z_m$ — стартовое значение прохода по $x$; правая граница $u(t,1,y,z)$ не требуется (схема односторонняя).
• По $y$ (прогонка подсхемы ②, ГУ 1-го рода): $u(t,x,0,z)=x_j z_m\Rightarrow\alpha_1=0,\ \beta_1=x_j z_m$; правая граница $u(t,x,1,z)=4x_j z_m$ замыкает прогонку.
• По $z$ (прогонка подсхемы ③, ГУ 1-го рода): $u(t,x,y,0)=x_j y_k\Rightarrow\alpha_1=0,\ \beta_1=x_j y_k$; правая граница $u(t,x,y,1)=5x_j y_k$.
• Начальное условие: $u_{j,k,m}^{\,0}=1$ во всех внутренних узлах.

6. Порядок аппроксимации

По времени корректор делает полный шаг $\Delta t$ от $u^{\,n}$ с операторами в средней точке $t^{\,n+1/2}$ — центральная разность, второй порядок $O(\Delta t^2)$. Диффузия по $y,z$ — центральные разности второго порядка $O(h_y^2),\,O(h_z^2)$. Конвекция аппроксимирована односторонними разностями первого порядка: по $x$ — $O(h_x)$, по $y$ — $O(h_y)$. По $y$ присутствуют и конвекция (1-й порядок), и диффузия (2-й порядок), поэтому суммарный порядок по $y$ — первый. Итог:

$$\psi=O\big(\Delta t^{2},\,h_x,\,h_y,\,h_z^{2}\big).$$

(Если бы конвекцию аппроксимировали центральными разностями, порядок по $x,y$ поднялся бы до второго: $O(\Delta t^2,h_x^2,h_y^2,h_z^2)$.) Схема абсолютно устойчива по всем неявно расщеплённым направлениям.

Уравнение и условия

$$\frac{\partial u}{\partial t}+7\,\frac{\partial u}{\partial x}-8\,\frac{\partial u}{\partial y}=0{,}2\,\frac{\partial^2 u}{\partial x^2}+0{,}3\,\frac{\partial^2 u}{\partial y^2}+t\sin(xy).$$

Перенесём конвективные члены вправо (приведём к «эволюционному» виду $u_t=\dots$):

$$\frac{\partial u}{\partial t}=0{,}2\,\frac{\partial^2 u}{\partial x^2}-7\,\frac{\partial u}{\partial x}+0{,}3\,\frac{\partial^2 u}{\partial y^2}+8\,\frac{\partial u}{\partial y}+t\sin(xy).$$

Итак, коэффициенты диффузии $\sigma_x=0{,}2$, $\sigma_y=0{,}3$; коэффициенты конвекции (со своими знаками после переноса) $C_x=-7$, $C_y=+8$; свободный член $f(t,x,y)=t\sin(xy)$. Все четыре граничных условия — 2-го рода (Неймана):

$$u_x(t,0,y)=t,\quad u_x(t,1,y)=2t,\qquad u_y(t,x,0)=t,\quad u_y(t,x,1)=2t;\qquad u(0,x,y)=y\,e^{x}.$$

Сетка $u_{j,k}^n=u(t^n,x_j,y_k)$, $x_j=(j-1)h_x$, $y_k=(k-1)h_y$, $t^n=n\Delta t$. Разностные операторы:

$$\Lambda_{xx}u_{j,k}=\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^2},\qquad \Lambda_{yy}u_{j,k}=\frac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^2},$$

первые производные — центральными разностями $\dfrac{u_{j+1,k}-u_{j-1,k}}{2h_x}$, $\dfrac{u_{j,k+1}-u_{j,k-1}}{2h_y}$.

Идея схемы расщепления (метода дробных шагов)

В отличие от СПН, схема расщепления делит шаг $\Delta t$ точкой $t^{n+1/2}$ так, что весь оператор по $x$ (диффузия $\sigma_x$ с полным коэффициентом + конвекция $C_x$) относится к первой подсхеме (неявной по $x$), а весь оператор по $y$ ($\sigma_y$ + $C_y$) — ко второй (неявной по $y$). Оператор «чужого» направления в данной подсхеме отсутствует (не берётся явно, как в СПН) — направления полностью расщеплены. Свободный член $f$ относим в первую подсхему. Каждая подсхема одномерна, трёхдиагональна и решается прогонкой.

$$\textcircled{1}\ (n\to n{+}\tfrac12,\ \text{неявно по }x):\quad \frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}=0{,}2\,\Lambda_{xx}u_{j,k}^{n+1/2}-7\,\frac{u_{j+1,k}^{n+1/2}-u_{j-1,k}^{n+1/2}}{2h_x}+t^{n+1/2}\sin(x_jy_k),$$

$$\textcircled{2}\ (n{+}\tfrac12\to n{+}1,\ \text{неявно по }y):\quad \frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}=0{,}3\,\Lambda_{yy}u_{j,k}^{n+1}+8\,\frac{u_{j,k+1}^{n+1}-u_{j,k-1}^{n+1}}{2h_y}.$$

Подсхема ① — приведение к прогонке по $j$ (направление $x$)

Умножаем на $\Delta t$ и переносим все неизвестные слоя $n{+}\tfrac12$ влево. Получаем трёхдиагональную систему вида $a_j u_{j+1,k}^{n+1/2}+b_j u_{j,k}^{n+1/2}+c_j u_{j-1,k}^{n+1/2}=\xi_{j,k}$. Вклад члена $-7\,\partial u/\partial x$ при центральной разности даёт добавки $-C_x\dfrac{\Delta t}{2h_x}$ к коэффициенту при $u_{j+1}$ и $+C_x\dfrac{\Delta t}{2h_x}$ при $u_{j-1}$ (при $C_x=-7$ это $+\tfrac{7\Delta t}{2h_x}$ и $-\tfrac{7\Delta t}{2h_x}$ соответственно):

$$a_j=-\frac{0{,}2\,\Delta t}{h_x^2}+\frac{7\Delta t}{2h_x}=-\frac{0{,}2\,\Delta t}{h_x^2}+\frac{3{,}5\,\Delta t}{h_x}\quad(\text{при }u_{j+1,k}),$$

$$c_j=-\frac{0{,}2\,\Delta t}{h_x^2}-\frac{7\Delta t}{2h_x}=-\frac{0{,}2\,\Delta t}{h_x^2}-\frac{3{,}5\,\Delta t}{h_x}\quad(\text{при }u_{j-1,k}),$$

$$b_j=1+\frac{0{,}4\,\Delta t}{h_x^2}\quad(\text{при }u_{j,k}),\qquad \xi_{j,k}=u_{j,k}^{n}+\Delta t\,t^{n+1/2}\sin(x_jy_k).$$

Подсхема ② — приведение к прогонке по $k$ (направление $y$)

Аналогично, система $\tilde a_k u_{j,k+1}^{n+1}+\tilde b_k u_{j,k}^{n+1}+\tilde c_k u_{j,k-1}^{n+1}=\tilde\xi_{j,k}$. Конвекция $C_y=+8$ даёт добавки $-8\dfrac{\Delta t}{2h_y}=-\dfrac{4\Delta t}{h_y}$ при $u_{j,k+1}$ и $+\dfrac{4\Delta t}{h_y}$ при $u_{j,k-1}$:

$$\tilde a_k=-\frac{0{,}3\,\Delta t}{h_y^2}-\frac{4\Delta t}{h_y}\quad(\text{при }u_{j,k+1}),\qquad \tilde c_k=-\frac{0{,}3\,\Delta t}{h_y^2}+\frac{4\Delta t}{h_y}\quad(\text{при }u_{j,k-1}),$$

$$\tilde b_k=1+\frac{0{,}6\,\Delta t}{h_y^2}\quad(\text{при }u_{j,k}),\qquad \tilde\xi_{j,k}=u_{j,k}^{n+1/2}.$$

Проверка сходимости прогонки (диагональное преобладание)

Подсхема ①. Если конвекция не превосходит диффузию ($\tfrac{3{,}5\Delta t}{h_x}\le\tfrac{0{,}2\Delta t}{h_x^2}$, т.е. $h_x\le\dfrac{2\sigma_x}{|C_x|}=\dfrac{0{,}4}{7}\approx0{,}057$), оба коэффициента $a_j,c_j$ остаются отрицательными и

$$|a_j|+|c_j|=\frac{0{,}4\,\Delta t}{h_x^2}\;<\;1+\frac{0{,}4\,\Delta t}{h_x^2}=|b_j|\quad\Rightarrow\quad\text{преобладание выполнено (запас }=1).$$

При нарушении этого условия (сеточное число Рейнольдса $\mathrm{Re}_h=\tfrac{|C_x|h_x}{2\sigma_x}>1$) один из коэффициентов меняет знак и для гарантии устойчивости центральной аппроксимации конвекции требуется измельчение сетки $h_x\le 0{,}4/7$ (либо переход к разностям против потока).

Подсхема ②. Аналогично при $h_y\le\dfrac{2\sigma_y}{|C_y|}=\dfrac{0{,}6}{8}=0{,}075$ имеем $|\tilde a_k|+|\tilde c_k|=\dfrac{0{,}6\Delta t}{h_y^2}<|\tilde b_k|$ — преобладание выполнено.

Рекуррентное прогоночное соотношение

Прогонка по $x$ (подсхема ①): $\;u_{j,k}^{n+1/2}=\alpha_j\,u_{j+1,k}^{n+1/2}+\beta_j$, где

$$\alpha_j=\frac{-a_j}{b_j+c_j\alpha_{j-1}},\qquad \beta_j=\frac{\xi_{j,k}-c_j\beta_{j-1}}{b_j+c_j\alpha_{j-1}}.$$

Прогонка по $y$ (подсхема ②): $\;u_{j,k}^{n+1}=\tilde\alpha_k\,u_{j,k+1}^{n+1}+\tilde\beta_k$, $\;\tilde\alpha_k=\dfrac{-\tilde a_k}{\tilde b_k+\tilde c_k\tilde\alpha_{k-1}}$, $\;\tilde\beta_k=\dfrac{\tilde\xi_{j,k}-\tilde c_k\tilde\beta_{k-1}}{\tilde b_k+\tilde c_k\tilde\alpha_{k-1}}.$

Благодаря диагональному преобладанию $|\alpha_j|<1$ при $j\ge2$, поэтому знаменатели и правое замыкание корректны.

Стартовые $\alpha_1,\beta_1$ из краевых условий (оба направления — 2-й род)

Подсхема ① (левое ГУ по $x$, $y=$const): $u_x(t,0,y)=t$ аппроксимируем односторонней разностью 1-го порядка:

$$\frac{u_{2,k}^{n+1/2}-u_{1,k}^{n+1/2}}{h_x}=t^{n+1/2}\;\Rightarrow\; u_{1,k}^{n+1/2}=u_{2,k}^{n+1/2}-h_x\,t^{n+1/2}.$$

Сравнивая с $u_{1,k}^{n+1/2}=\alpha_1 u_{2,k}^{n+1/2}+\beta_1$, получаем

$$\boxed{\alpha_1=1,\qquad \beta_1=-h_x\,t^{n+1/2}.}$$

Правое ГУ по $x$ ($u_x(t,1,y)=2t$): из $\dfrac{u_{N_x,k}^{n+1/2}-u_{N_x-1,k}^{n+1/2}}{h_x}=2t^{n+1/2}$ и $u_{N_x-1,k}^{n+1/2}=\alpha_{N_x-1}u_{N_x,k}^{n+1/2}+\beta_{N_x-1}$ замыкаем правый узел:

$$u_{N_x,k}^{n+1/2}=\frac{\beta_{N_x-1}+2t^{n+1/2}h_x}{1-\alpha_{N_x-1}}.$$

Подсхема ② (левое ГУ по $y$, $x=$const): $u_y(t,x,0)=t$:

$$\frac{u_{j,2}^{n+1}-u_{j,1}^{n+1}}{h_y}=t^{n+1}\;\Rightarrow\; u_{j,1}^{n+1}=u_{j,2}^{n+1}-h_y\,t^{n+1},$$

откуда

$$\boxed{\tilde\alpha_1=1,\qquad \tilde\beta_1=-h_y\,t^{n+1}.}$$

Правое ГУ по $y$ ($u_y(t,x,1)=2t$): $\;u_{j,N_y}^{n+1}=\dfrac{\tilde\beta_{N_y-1}+2t^{n+1}h_y}{1-\tilde\alpha_{N_y-1}}.$

Алгоритм и порядок аппроксимации

1. Начальный слой $u_{j,k}^0=y_k e^{x_j}$.
2. Подсхема ①: для каждого $k$ прогонка по $x$ (старт $\alpha_1=1,\ \beta_1=-h_x t^{n+1/2}$, замыкание справа) $\to u^{n+1/2}$.
3. Подсхема ②: для каждого $j$ прогонка по $y$ (старт $\tilde\alpha_1=1,\ \tilde\beta_1=-h_y t^{n+1}$, замыкание справа) $\to u^{n+1}$.
4. Повторять по $n$.

Порядок $O(\Delta t,\,h_x^2,\,h_y^2)$ внутри области; у границ односторонняя аппроксимация ГУ 2-го рода локально даёт $O(h)$. Обе подсхемы (по диффузии) абсолютно устойчивы, центральная конвекция накладывает условие $h\le 2\sigma/|C|$.

Условие

Для уравнения

$$\frac{\partial u}{\partial t}=0{,}2\,\frac{\partial u}{\partial x}-0{,}1\,\frac{\partial u}{\partial y}+\sin x+\sin y$$

с граничными и начальным условиями

$$\begin{cases}u(t,0,y)=\sin y,\\ u(t,1,y)=\cos y,\end{cases}\qquad\begin{cases}u(t,x,0)=\sin x,\\ u(t,x,1)=\cos x,\end{cases}\qquad u(0,x,y)=0$$

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

1. Тип уравнения и канонический вид

В уравнении присутствуют только первые производные по пространственным переменным — это двумерное уравнение 1-го порядка (уравнение переноса). Приведём его к стандартному виду $\dfrac{\partial u}{\partial t}+v_1\dfrac{\partial u}{\partial x}+v_2\dfrac{\partial u}{\partial y}=f$, перенося конвективные члены налево:

$$\frac{\partial u}{\partial t}-0{,}2\,\frac{\partial u}{\partial x}+0{,}1\,\frac{\partial u}{\partial y}=\sin x+\sin y.$$

Отсюда скорости переноса и свободный член:

$$v_1=-0{,}2<0,\qquad v_2=+0{,}1>0,\qquad f(x,y)=\sin x+\sin y.$$

Замечание. Это уравнение переноса, поэтому метод дробных шагов с прогонкой (трёхдиагональные системы, коэффициенты $\alpha_1,\beta_1$) здесь не возникает — по условию Варианта 3 требуется явная схема, в которой новый слой находится прямой подстановкой.

2. Выбор граничных условий и направления разностей

Для уравнения 1-го порядка по каждому направлению нужно одно граничное условие — на входной (по потоку) границе; разность по этому направлению берётся «против потока» (upwind), что и даёт устойчивую явную схему. Характеристики: $\dot x=v_1$, $\dot y=v_2$.

• По $x$: $v_1=-0{,}2<0$ — поток идёт справа налево, информация входит со стороны $x=1$. Входная граница $x=1$: берём $u(t,1,y)=\cos y$ и правую разность $\dfrac{u_{j+1,k}-u_{j,k}}{h_x}$.
• По $y$: $v_2=+0{,}1>0$ — поток идёт снизу вверх, информация входит со стороны $y=0$. Входная граница $y=0$: берём $u(t,x,0)=\sin x$ и левую разность $\dfrac{u_{j,k}-u_{j,k-1}}{h_y}$.

Условия на выходных границах — $u(t,0,y)=\sin y$ и $u(t,x,1)=\cos x$ — для явной схемы переноса не используются (иначе задача оказалась бы переопределена; значения в этих узлах вычисляются по самой схеме).

Сетка: $x_j=(j-1)h_x$, $j=1,\dots,N_x$, $h_x=\dfrac{1}{N_x-1}$; $\;y_k=(k-1)h_y$, $k=1,\dots,N_y$, $h_y=\dfrac{1}{N_y-1}$; $\;t^{n}=n\,\Delta t$.

3. Явная разностная схема

Производную по времени аппроксимируем правой разностью (значения с известного слоя $n$), пространственные производные — выбранными upwind-разностями на слое $n$. Свободный член берём в узле $(x_j,y_k)$:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=0{,}2\,\frac{u_{j+1,k}^{n}-u_{j,k}^{n}}{h_x}-0{,}1\,\frac{u_{j,k}^{n}-u_{j,k-1}^{n}}{h_y}+\sin x_j+\sin y_k.$$

Начальное условие: $u_{j,k}^{0}=0$.

Граничные условия (входные, 1-го рода): $u_{N_x,k}^{\,n+1}=\cos y_k$ (граница $x=1$), $\;u_{j,1}^{\,n+1}=\sin x_j$ (граница $y=0$).

4. Рекуррентное соотношение

Все пространственные члены явной схемы взяты со слоя $n$, поэтому единственное неизвестное $u_{j,k}^{n+1}$ выражается прямой подстановкой (без прогонки):

$$\boxed{\,u_{j,k}^{n+1}=u_{j,k}^{n}+0{,}2\,\frac{\Delta t}{h_x}\bigl(u_{j+1,k}^{n}-u_{j,k}^{n}\bigr)-0{,}1\,\frac{\Delta t}{h_y}\bigl(u_{j,k}^{n}-u_{j,k-1}^{n}\bigr)+\Delta t\,\bigl(\sin x_j+\sin y_k\bigr)\,}$$

Формула даёт значения нового слоя для узлов $j=1,\dots,N_x-1$, $k=2,\dots,N_y$ (входная граница $x=1$, $j=N_x$, и входная граница $y=0$, $k=1$, берутся из граничных условий).

Удобно ввести числа Куранта $r_x=\dfrac{|v_1|\,\Delta t}{h_x}=\dfrac{0{,}2\,\Delta t}{h_x}$ и $r_y=\dfrac{|v_2|\,\Delta t}{h_y}=\dfrac{0{,}1\,\Delta t}{h_y}$. Тогда

$$u_{j,k}^{n+1}=(1-r_x-r_y)\,u_{j,k}^{n}+r_x\,u_{j+1,k}^{n}+r_y\,u_{j,k-1}^{n}+\Delta t\,(\sin x_j+\sin y_k).$$

5. Условие устойчивости на шаг

Запись схемы через $r_x,r_y$ показывает, что новое значение есть линейная комбинация трёх соседних. Схема устойчива (по принципу максимума / Куранта), когда все коэффициенты неотрицательны, то есть $1-r_x-r_y\ge 0$:

$$\boxed{\,r_x+r_y\le 1\quad\Longleftrightarrow\quad 0{,}2\,\frac{\Delta t}{h_x}+0{,}1\,\frac{\Delta t}{h_y}\le 1\,}$$

(суммарное число Куранта по двум направлениям не превосходит 1; тот же результат даёт фон-неймановский анализ 2D-схемы «против потока»). Схема условно устойчива: шаг $\Delta t$ ограничен сверху

$$\Delta t\le\frac{1}{\dfrac{0{,}2}{h_x}+\dfrac{0{,}1}{h_y}}.$$

6. Алгоритм

1. Задать $h_x,h_y$ и $\Delta t$, удовлетворяющий $0{,}2\frac{\Delta t}{h_x}+0{,}1\frac{\Delta t}{h_y}\le 1$; задать $N_x,N_y,N_t$.
2. Начальный слой: $u_{j,k}^{0}=0$ для всех $j,k$.
3. Цикл по времени $n=0,\dots,N_t-1$:
   • Для узлов $j=1,\dots,N_x-1$, $k=2,\dots,N_y$ — по рекуррентной формуле (используются $u_{j,k}^{n}$, $u_{j+1,k}^{n}$, $u_{j,k-1}^{n}$).
   • Входные границы: $u_{N_x,k}^{\,n+1}=\cos y_k$, $\;u_{j,1}^{\,n+1}=\sin x_j$.
4. Результат — массив $u_{j,k}^{n}$. Порядок аппроксимации $O(\Delta t,h_x,h_y)$ (первый порядок по времени и пространству).

Уравнение и граничные условия

$$2\,\frac{\partial^2 u}{\partial x^2}+8\,\frac{\partial^2 u}{\partial y^2}=20,\qquad x\in[0,1],\;y\in[0,1];$$

$$\begin{cases} u(0,y)=y^2 \\ u(1,y)=y^2+1 \end{cases}\qquad \begin{cases} u(x,0)=x^2 \\ u(x,1)=x^2+1 \end{cases}.$$

Это эллиптическая (стационарная) краевая задача с граничными условиями первого рода (значения $u$ заданы на всех четырёх сторонах квадрата). По условию варианта 3 — метод установления со схемой переменных направлений (СПН, схема Писмена–Рэкфорда).

Шаг 1. Идея метода установления

Стационарную задачу превращаем в нестационарную, вводя фиктивную производную по «времени» $\partial u/\partial t$: при $t\to\infty$ решение «устанавливается» к искомому ($\partial u/\partial t\to0$). Вторые производные оставляем справа с положительными коэффициентами, а источник переносим вправо со знаком минус:

$$\frac{\partial u}{\partial t}=2\,\frac{\partial^2 u}{\partial x^2}+8\,\frac{\partial^2 u}{\partial y^2}-20.$$

Проверка: при $\partial u/\partial t=0$ получаем $2u_{xx}+8u_{yy}-20=0$, т.е. исходное $2u_{xx}+8u_{yy}=20$. ✓ Здесь $\sigma_x=2>0$, $\sigma_y=8>0$ — коэффициенты диффузии положительны, уравнение параболического типа, процесс установления устойчив. Роль «времени» играют итерации $s\to s+1$ с промежуточным полуслоем $s+1/2$, шаг по «времени» — $\Delta t$.

Сетка $x_j=(j-1)h_x$, $y_k=(k-1)h_y$; разностные операторы вторых производных (вторая разность $=(u_{j+1}-2u_j+u_{j-1})/h^2$):

$$\Lambda_{xx}u_{j,k}=\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^2},\qquad \Lambda_{yy}u_{j,k}=\frac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^2}.$$

Шаг 2. Схема переменных направлений (каноническая форма)

Записываем СПН в форме семинара С9. Шаг $\Delta t$ расщепляется промежуточным полуслоем $s+1/2$ на две подсхемы; каждый диффузионный оператор берётся с половинным коэффициентом ($\sigma_x/2=1$, $\sigma_y/2=4$) и присутствует в обеих подсхемах — «своё» направление неявно, «чужое» берётся на уже посчитанном слое. В обеих подсхемах стоит полный шаг $\Delta t$.

Ключевой момент (отличие от шаблона С9). В семинаре С9 свободный член $f$ для нестационарной задачи ставится только во вторую подсхему (на полуслое $s+1/2$) — это даёт 2-й порядок по времени при счёте переходного процесса. Но здесь задача стационарная, и важна лишь неподвижная точка итераций. Чтобы она не зависела от $\Delta t$ и совпала с точным стационарным решением, свободный член $-20$ нужно расщепить симметрично — по $-10$ в каждую подсхему ($-f/2$, $f=20$). Иначе, поместив весь $-20$ только во вторую подсхему, в неподвижной точке ($u^s=u^{s+1}=U,\ u^{s+1/2}=V$) вычитание ② $-$ ① даёт $2(U-V)/\Delta t=-20$, т.е. $V=U+10\,\Delta t$ во внутренних узлах при точных границах — полуслой смещён относительно слоя на $\propto\Delta t$, и установление сходится к решению с ошибкой $\propto\Delta t$ (см. численную проверку).

$$\boxed{\text{①}\quad \frac{u_{j,k}^{\,s+1/2}-u_{j,k}^{\,s}}{\Delta t}=1\cdot\Lambda_{xx}u_{j,k}^{\,s+1/2}+4\cdot\Lambda_{yy}u_{j,k}^{\,s}-10}$$

$$\boxed{\text{②}\quad \frac{u_{j,k}^{\,s+1}-u_{j,k}^{\,s+1/2}}{\Delta t}=1\cdot\Lambda_{xx}u_{j,k}^{\,s+1/2}+4\cdot\Lambda_{yy}u_{j,k}^{\,s+1}-10}$$

Подсхема ① — неявна по $x$ (уровень $s+1/2$), явна по $y$ (уровень $s$). Подсхема ② — неявна по $y$ (уровень $s+1$), а член $\Lambda_{xx}$ берётся на уже найденном полуслое $s+1/2$.

Согласованность подсхем. Сложим ① и ②:

$$\frac{u_{j,k}^{\,s+1}-u_{j,k}^{\,s}}{\Delta t}=2\,\Lambda_{xx}u_{j,k}^{\,s+1/2}+4\,\Lambda_{yy}u_{j,k}^{\,s}+4\,\Lambda_{yy}u_{j,k}^{\,s+1}-20.$$

Правая часть симметрична относительно полуслоя $s+1/2$, поэтому разность по «времени» — центральная, и схема имеет 2-й порядок по «времени»: $O(\Delta t^2,\,h_x^2,\,h_y^2)$. В стационаре ($u^s=u^{s+1/2}=u^{s+1}=U$) получаем $2\Lambda_{xx}U+8\Lambda_{yy}U-20=0$ — точную разностную аппроксимацию исходного уравнения, не содержащую $\Delta t$. Именно симметрия операторов и источника гарантирует корректность неподвижной точки.

Шаг 3. Приведение подсхем к виду для прогонки

Подсхема ① (неизвестные — вдоль строки $k$, индекс $j$). Раскрываем $\Lambda_{xx}u^{s+1/2}$, переносим неявные члены влево:

$$-\frac{\Delta t}{h_x^2}\,u_{j-1,k}^{\,s+1/2}+\Bigl(1+\frac{2\,\Delta t}{h_x^2}\Bigr)u_{j,k}^{\,s+1/2}-\frac{\Delta t}{h_x^2}\,u_{j+1,k}^{\,s+1/2}=u_{j,k}^{\,s}+4\,\Delta t\,\Lambda_{yy}u_{j,k}^{\,s}-10\,\Delta t.$$

Это трёхдиагональная система $a_j u_{j-1,k}^{\,s+1/2}+b_j u_{j,k}^{\,s+1/2}+c_j u_{j+1,k}^{\,s+1/2}=\xi_{j,k}$ с коэффициентами

$$a_j=c_j=-\frac{\Delta t}{h_x^2},\qquad b_j=1+\frac{2\,\Delta t}{h_x^2},\qquad \xi_{j,k}=u_{j,k}^{\,s}+4\,\Delta t\,\Lambda_{yy}u_{j,k}^{\,s}-10\,\Delta t.$$

Подсхема ② (неизвестные — вдоль столбца $j$, индекс $k$). Аналогично с оператором $\Lambda_{yy}$ и коэффициентом $4$:

$$-\frac{4\,\Delta t}{h_y^2}\,u_{j,k-1}^{\,s+1}+\Bigl(1+\frac{8\,\Delta t}{h_y^2}\Bigr)u_{j,k}^{\,s+1}-\frac{4\,\Delta t}{h_y^2}\,u_{j,k+1}^{\,s+1}=u_{j,k}^{\,s+1/2}+\Delta t\,\Lambda_{xx}u_{j,k}^{\,s+1/2}-10\,\Delta t,$$

$$\tilde a_k=\tilde c_k=-\frac{4\,\Delta t}{h_y^2},\qquad \tilde b_k=1+\frac{8\,\Delta t}{h_y^2},\qquad \tilde\xi_{j,k}=u_{j,k}^{\,s+1/2}+\Delta t\,\Lambda_{xx}u_{j,k}^{\,s+1/2}-10\,\Delta t.$$

Шаг 4. Проверка сходимости прогонки (диагональное преобладание)

Подсхема ①:

$$|a_j|+|c_j|=\frac{2\,\Delta t}{h_x^2}<1+\frac{2\,\Delta t}{h_x^2}=|b_j|.$$

Подсхема ②:

$$|\tilde a_k|+|\tilde c_k|=\frac{8\,\Delta t}{h_y^2}<1+\frac{8\,\Delta t}{h_y^2}=|\tilde b_k|.$$

В обоих случаях $|b|-(|a|+|c|)=1>0$ — строгое диагональное преобладание. Абсолютный запас, равный $1$, обеспечивается членом $+1$, который вносит производная по «времени» (он же делает схему безусловно устойчивой). Поэтому достаточное условие сходимости прогонки выполнено при любых $\Delta t,h_x,h_y>0$ — прогонка сходится безусловно.

Шаг 5. Итерационные (прогоночные) соотношения, $\alpha_1,\beta_1$

Подсхема ① по $x$. Прямой ход: $u_{j,k}^{\,s+1/2}=\alpha_{j+1}u_{j+1,k}^{\,s+1/2}+\beta_{j+1}$ с

$$\alpha_{j+1}=\frac{-c_j}{b_j+a_j\alpha_j},\qquad \beta_{j+1}=\frac{\xi_{j,k}-a_j\beta_j}{b_j+a_j\alpha_j}.$$

Левая граница $x=0$ (ГУ 1-го рода): $u_{1,k}^{\,s+1/2}=u(0,y_k)=y_k^2$. Записывая это как $u_{1,k}=\alpha_1 u_{2,k}+\beta_1$, получаем

$$\alpha_1=0,\qquad \beta_1=y_k^2.$$

Правая граница $x=1$ задаёт старт обратного хода: $u_{N_x,k}^{\,s+1/2}=u(1,y_k)=y_k^2+1$.

Подсхема ② по $y$. Аналогично $u_{j,k}^{\,s+1}=\tilde\alpha_{k+1}u_{j,k+1}^{\,s+1}+\tilde\beta_{k+1}$. Нижняя граница $y=0$: $u_{j,1}^{\,s+1}=u(x_j,0)=x_j^2$, поэтому

$$\tilde\alpha_1=0,\qquad \tilde\beta_1=x_j^2;$$

верхняя граница $y=1$: $u_{j,N_y}^{\,s+1}=u(x_j,1)=x_j^2+1$ (старт обратного хода).

Шаг 6. Начальное приближение и условие окончания итераций

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

$$u_{j,k}^{\,0}=f=20$$

во всех внутренних узлах, с обязательной расстановкой точных граничных значений на всех четырёх сторонах ($u(0,y)=y^2$ и т.д.). Поскольку граничные условия не зависят от «времени», установление сходится к стационарному решению независимо от выбора старта (подойдёт и $u^0_{j,k}=0$ или $x_jy_k$); канонический выбор $u^0=f$ лишь удобен и единообразен.

Условие окончания итераций. Канон гл. 10 — норма среднеквадратичного типа разности соседних приближений:

$$\big\|u^{\,s+1}-u^{\,s}\big\|=\sqrt{h_xh_y\sum_{j,k}\big(u_{j,k}^{\,s+1}-u_{j,k}^{\,s}\big)^2}\le\varepsilon$$

(равносильно можно использовать максимум-норму $\max_{j,k}\big|u_{j,k}^{\,s+1}-u_{j,k}^{\,s}\big|<\varepsilon$).

Шаг 7. Алгоритм

1. Задать $u_{j,k}^{0}=f=20$ (свободный член), выставить точные граничные значения.
2. ① для каждой строки $k$: прогонка по $x$ ($\alpha_1=0,\ \beta_1=y_k^2$; справа $u_{N_x,k}=y_k^2+1$) $\to u^{\,s+1/2}$.
3. ② для каждого столбца $j$: прогонка по $y$ ($\tilde\alpha_1=0,\ \tilde\beta_1=x_j^2$; сверху $u_{j,N_y}=x_j^2+1$) $\to u^{\,s+1}$.
4. Проверить $\sqrt{h_xh_y\sum_{j,k}(u^{\,s+1}_{j,k}-u^{\,s}_{j,k})^2}\le\varepsilon$: если нет — $s:=s+1$ и к п.2; если да — стоп.

Контроль точности (численная проверка)

Аналитическое решение задачи — $u(x,y)=x^2+y^2$: действительно $2\cdot2+8\cdot2=20$, и все четыре ГУ выполняются. Для $u=x^2+y^2$ трёхточечная вторая разность точна, поэтому корректная схема обязана сходиться к $x^2+y^2$ с точностью до округления при любом $\Delta t$. Прямой прогон (сетка $21\times21$, порог по изменению $10^{-13}$) даёт:

при симметричном расщеплении источника ($-10$ в каждую подсхему) и коэффициентах $\sigma/2$ в обеих подсхемах ошибка относительно $x^2+y^2$ составляет $\approx5{,}8\cdot10^{-12}$ при $\Delta t=0{,}1$; $2{,}8\cdot10^{-12}$ при $\Delta t=0{,}05$; $1{,}1\cdot10^{-12}$ при $\Delta t=0{,}02$; $5{,}4\cdot10^{-13}$ при $\Delta t=0{,}01$ — не зависит от $\Delta t$ и определяется лишь машинной точностью;

при несимметричном расщеплении (весь $-20$ во ②) ошибка $\propto\Delta t$: $0{,}80$; $0{,}40$; $0{,}16$; $0{,}08$ соответственно — к точному решению не сходится. Это подтверждает необходимость симметричного расщепления источника.

Уравнение и условия

\n

Трёхмерное (по координатам $x,y,z$) нестационарное уравнение конвекции-диффузии с линейной реакцией:

\n$$\\frac{\\partial u}{\\partial t}+5\\left(\\frac{\\partial u}{\\partial x}+\\frac{\\partial u}{\\partial y}-\\frac{\\partial u}{\\partial z}\\right)=7\\,\\frac{\\partial^2 u}{\\partial z^2}-u,\\qquad u=u(t,x,y,z).$$\n

Граничные и начальное условия:

\n$$u(t,0,y,z)=0,\\quad u(t,x,0,z)=0,\\quad u(t,x,y,0)=0,\\quad u(t,x,y,1)=x,\\quad u(0,x,y,z)=0.$$\n

Сетка: $j\\leftrightarrow x$, $k\\leftrightarrow y$, $m\\leftrightarrow z$; шаги $h_x,h_y,h_z$ и $\\Delta t$. Вводим разностный оператор второй производной по $z$ (вторая разность):

\n$$\\Lambda_{zz}u_{j,k,m}=\\frac{u_{j,k,m+1}-2u_{j,k,m}+u_{j,k,m-1}}{h_z^{2}}.$$\n\n

1. Анализ членов и выбор конечных разностей

\n

Запишем уравнение в форме $u_t+v_x u_x+v_y u_y+v_z u_z=\\sigma_z u_{zz}-k u$:

\n$$v_x=+5,\\qquad v_y=+5,\\qquad v_z=-5,\\qquad \\sigma_z=7,\\qquad k=+1.$$\n

Вторая производная присутствует только по $z$; по $x$ и $y$ уравнение первого порядка (чистая конвекция). Конвективные члены оставляем в левой части и аппроксимируем односторонними разностями по знаку скорости (как в семинаре C7: при $v>0$ — левая разность и левое ГУ, при $v<0$ — правая разность и правое ГУ):

\n

\n
• $v_x=+5>0$ ⇒ по $x$ левая разность $\\dfrac{u_{j,k,m}-u_{j-1,k,m}}{h_x}$, рабочее ГУ — левое $u(t,0,y,z)=0$;
\n
• $v_y=+5>0$ ⇒ по $y$ левая разность $\\dfrac{u_{j,k,m}-u_{j,k-1,m}}{h_y}$, рабочее ГУ — левое $u(t,x,0,z)=0$;
\n
• $v_z=-5<0$ ⇒ для члена $-5\\,u_z$ — правая (forward) разность $-5\\,\\dfrac{u_{j,k,m+1}-u_{j,k,m}}{h_z}$; вторая производная по $z$ — центральная $\\Lambda_{zz}$; по $z$ задача второго порядка, поэтому используются оба ГУ: $u(t,x,y,0)=0$ и $u(t,x,y,1)=x$.
\n

\n\n

2. Структура схемы предиктор-корректор (каноническая форма)

\n

Согласно учебнику (гл.9, п.8, ур. (9.17)–(9.21)) для трёхмерной задачи интервал $\\delta t$ между $t^{n}$ и $t^{n+1}$ расщепляется пополам точкой $t^{n+1/2}$, а полушаг $\\Delta t/2$ дополнительно делится на три равные части ($t^{n+1/6}$, $t^{n+1/3}$). Схема состоит из четырёх подсхем: трёх неявных подсхем предиктора (по одному пространственному направлению в каждой, каждая на полушаге $\\Delta t/2$) и одного явного корректора на целом шаге $\\Delta t$.

\n

Роль предиктора (9.17)–(9.19) — обеспечить абсолютную устойчивость; роль корректора (9.20) — повысить порядок аппроксимации по времени.

\n\n

Предиктор. Подсхема ① — неявно по $x$ (полушаг $\\Delta t/2$, $n\\to n+1/6$):

\n$$\\frac{u_{j,k,m}^{\\,n+1/6}-u_{j,k,m}^{\\,n}}{\\Delta t/2}+5\\,\\frac{u_{j,k,m}^{\\,n+1/6}-u_{j-1,k,m}^{\\,n+1/6}}{h_x}=0.$$\n\n

Предиктор. Подсхема ② — неявно по $y$ (полушаг $\\Delta t/2$, $n+1/6\\to n+1/3$):

\n$$\\frac{u_{j,k,m}^{\\,n+1/3}-u_{j,k,m}^{\\,n+1/6}}{\\Delta t/2}+5\\,\\frac{u_{j,k,m}^{\\,n+1/3}-u_{j,k-1,m}^{\\,n+1/3}}{h_y}=0.$$\n\n

Предиктор. Подсхема ③ — неявно по $z$ (полушаг $\\Delta t/2$, $n+1/3\\to n+1/2$; диффузия, $z$-конвекция и реакция):

\n$$\\frac{u_{j,k,m}^{\\,n+1/2}-u_{j,k,m}^{\\,n+1/3}}{\\Delta t/2}\n=7\\,\\Lambda_{zz}u_{j,k,m}^{\\,n+1/2}\n-5\\,\\frac{u_{j,k,m+1}^{\\,n+1/2}-u_{j,k,m}^{\\,n+1/2}}{h_z}\n-u_{j,k,m}^{\\,n+1/2}.$$\n

Последовательное решение ①→②→③ даёт промежуточную оценку $u_{j,k,m}^{\\,n+1/2}$.

\n\n

Корректор ④ — ЯВНОЕ соотношение на целом шаге $\\Delta t$ ($n\\to n+1$):

\n$$\\frac{u_{j,k,m}^{\\,n+1}-u_{j,k,m}^{\\,n}}{\\Delta t}\n=7\\,\\Lambda_{zz}u_{j,k,m}^{\\,n+1/2}\n-5\\,\\frac{u_{j,k,m}^{\\,n+1/2}-u_{j-1,k,m}^{\\,n+1/2}}{h_x}\n-5\\,\\frac{u_{j,k,m}^{\\,n+1/2}-u_{j,k-1,m}^{\\,n+1/2}}{h_y}\n+5\\,\\frac{u_{j,k,m+1}^{\\,n+1/2}-u_{j,k,m}^{\\,n+1/2}}{h_z}\n-u_{j,k,m}^{\\,n+1/2}.$$\n

Вся правая часть взята по уже найденному промежуточному слою $u^{n+1/2}$, поэтому корректор явный и реализуется прямым пересчётом (ур. (9.21)):

\n$$\\boxed{\\,u_{j,k,m}^{\\,n+1}=u_{j,k,m}^{\\,n}+\\Delta t\\!\\left[7\\,\\Lambda_{zz}u^{\\,n+1/2}-5\\,\\frac{u_{j,k,m}^{\\,n+1/2}-u_{j-1,k,m}^{\\,n+1/2}}{h_x}-5\\,\\frac{u_{j,k,m}^{\\,n+1/2}-u_{j,k-1,m}^{\\,n+1/2}}{h_y}+5\\,\\frac{u_{j,k,m+1}^{\\,n+1/2}-u_{j,k,m}^{\\,n+1/2}}{h_z}-u_{j,k,m}^{\\,n+1/2}\\right].}$$\n

Левая разность по времени $\\dfrac{u^{n+1}-u^{n}}{\\Delta t}$ центрирована относительно $t^{n+1/2}$ — это центральная конечная разность 2-го порядка по времени, что и повышает порядок схемы.

\n\n

3. Рекуррентные соотношения по подсхемам

\n\n

Подсхема ① (по $x$) — прямая двухточечная рекуррента. Член по $x$ первого порядка, левая разность ⇒ узел выражается через уже найденного «слева» соседа $u_{j-1}$; расчёт идёт слева направо по $j$. Умножая на $\\Delta t/2$:

\n$$u_{j,k,m}^{\\,n+1/6}\\Big(1+\\frac{5\\Delta t}{2h_x}\\Big)-\\frac{5\\Delta t}{2h_x}\\,u_{j-1,k,m}^{\\,n+1/6}=u_{j,k,m}^{\\,n}\n\\;\\Longrightarrow\\;\nu_{j,k,m}^{\\,n+1/6}=\\frac{5\\Delta t\\,u_{j-1,k,m}^{\\,n+1/6}+2h_x\\,u_{j,k,m}^{\\,n}}{5\\Delta t+2h_x}.$$\n

Запуск: левое ГУ $u_{0,k,m}^{\\,n+1/6}=0$.

\n\n

Подсхема ② (по $y$) — прямая двухточечная рекуррента (слева направо по $k$, запуск $u_{j,0,m}^{\\,n+1/3}=0$):

\n$$u_{j,k,m}^{\\,n+1/3}=\\frac{5\\Delta t\\,u_{j,k-1,m}^{\\,n+1/3}+2h_y\\,u_{j,k,m}^{\\,n+1/6}}{5\\Delta t+2h_y}.$$\n\n

Подсхема ③ (по $z$) — трёхдиагональная система, решается прогонкой. Фиксируем $(j,k)$, переносим неявные члены налево и умножаем на $\\Delta t/2$:

\n$$a_m\\,u_{j,k,m-1}^{\\,n+1/2}+b_m\\,u_{j,k,m}^{\\,n+1/2}+c_m\\,u_{j,k,m+1}^{\\,n+1/2}=\\xi_{j,k,m},$$\n$$a_m=-\\frac{7\\,\\Delta t}{2h_z^{2}},\\qquad\nb_m=1+\\frac{7\\,\\Delta t}{h_z^{2}}+\\frac{5\\,\\Delta t}{2h_z}+\\frac{\\Delta t}{2},\\qquad\nc_m=-\\frac{7\\,\\Delta t}{2h_z^{2}}-\\frac{5\\,\\Delta t}{2h_z},$$\n$$\\xi_{j,k,m}=u_{j,k,m}^{\\,n+1/3}.$$\n

(Член $-5\\,u_z$ при правой разности $\\big(u_{m+1}-u_m\\big)/h_z$ даёт $+\\dfrac{5\\Delta t}{2h_z}$ в $b_m$ и $-\\dfrac{5\\Delta t}{2h_z}$ в $c_m$; реакция $-u$ даёт $+\\dfrac{\\Delta t}{2}$ в $b_m$ — оба усиливают диагональ.) Рекуррента прогонки по $m$:

\n$$u_{j,k,m}^{\\,n+1/2}=\\alpha_{m+1}\\,u_{j,k,m+1}^{\\,n+1/2}+\\beta_{m+1},\\qquad\n\\alpha_{m+1}=\\frac{-c_m}{b_m+a_m\\alpha_m},\\quad\n\\beta_{m+1}=\\frac{\\xi_{j,k,m}-a_m\\beta_m}{b_m+a_m\\alpha_m}.$$\n

ГУ по $z$ (1-го рода): при $z=0$ — $u_{j,k,0}^{\\,n+1/2}=0\\Rightarrow\\alpha_1=0,\\ \\beta_1=0$; при $z=1$ — $u_{j,k,N_z}^{\\,n+1/2}=x_j$. Достаточное условие устойчивости прогонки выполнено: $c_m<0$ и $b_m-\\big(|a_m|+|c_m|\\big)=1+\\dfrac{\\Delta t}{2}>0$ (проверено символически), т.е. имеет место диагональное преобладание.

\n\n

Корректор ④ — прямой явный пересчёт. Отдельная система не решается: для каждого узла $(j,k,m)$ значение $u_{j,k,m}^{\\,n+1}$ вычисляется по выписанному выше боксу подстановкой $u^{n+1/2}$.

\n\n

4. Порядок аппроксимации

\n

Правая часть корректора аппроксимирована относительно $t^{n+1/2}$, левая временная разность центральная ⇒ второй порядок по времени. Вторая производная по $z$ — центральная разность $O(h_z^2)$; конвективные члены по $x,y,z$ — односторонние разности первого порядка $O(h_x),O(h_y),O(h_z)$. Итог:

\n$$\\psi=O\\big(\\Delta t^{2},\\;h_x,\\;h_y,\\;h_z\\big).$$\n

Если конвективные члены аппроксимировать центральными разностями, порядок по пространству повышается до $O(h_x^2,h_y^2,h_z^2)$, но возможна потеря устойчивости при больших скоростях; в курсе для конвекции используется устойчивая односторонняя разность по знаку скорости.

\n\n

5. Алгоритм на шаге $n\\to n+1$

\n

\n
1. $u^{0}_{j,k,m}=0$ (начальное условие).
\n
2. Предиктор ① (x): для каждой $x$-линии — прямая рекуррента слева направо (запуск $u_{0,k,m}=0$) ⇒ $u^{n+1/6}$.
\n
3. Предиктор ② (y): для каждой $y$-линии — прямая рекуррента слева направо (запуск $u_{j,0,m}=0$) ⇒ $u^{n+1/3}$.
\n
4. Предиктор ③ (z): для каждой $z$-линии $(j,k)$ — прогонка по $m$ (ГУ $z=0:0$, $z=1:x_j$) ⇒ $u^{n+1/2}$.
\n
5. Корректор ④: явный пересчёт $u^{n+1}$ по формуле (9.21) во всех внутренних узлах; на гранях — ГУ.
\n
6. Переход к слою $n+1$.
\n

\n

Уравнение и условия

$$\frac{\partial u}{\partial t}+4\,\frac{\partial u}{\partial x}=6t\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+txy,$$

$$u(t,0,y)=t,\quad u(t,1,y)=2t,\qquad u(t,x,0)=t,\quad u(t,x,1)=2t,\qquad u(t{=}0,x,y)=0.$$

Приведём к канонической форме $u_t=\sigma_x u_{xx}+\sigma_y u_{yy}+C_x u_x+C_y u_y+f$:

$$\frac{\partial u}{\partial t}=6t\,\frac{\partial^2 u}{\partial x^2}+6t\,\frac{\partial^2 u}{\partial y^2}-4\,\frac{\partial u}{\partial x}+txy.$$

Коэффициенты: диффузия $\sigma_x=\sigma_y=6t$, конвекция $C_x=-4$, $C_y=0$, источник $f=txy$. Все четыре граничных условия — 1-го рода (Дирихле). Сетка $u_{j,k}^n=u(t^n,x_j,y_k)$, $j=1..N_x$, $k=1..N_y$; шаги $h_x,h_y,\Delta t$; $x_j=(j-1)h_x$, $y_k=(k-1)h_y$.

Идея схемы переменных направлений (СПН)

Следуем канонической записи СПН (семинар 9). Шаг $\Delta t$ делится точкой $t^{n+1/2}=t^n+\tfrac{\Delta t}{2}$ на две подсхемы; знаменатель обеих подсхем равен $\Delta t$. В каждой подсхеме неявен только один разностный оператор (вдоль одного направления) — система трёхдиагональна и решается прогонкой; второй оператор берётся явно с соседнего слоя.

Диффузия расщепляется поровну: $\tfrac{\sigma_x}{2}=\tfrac{\sigma_y}{2}=3t^{n+1/2}$. Оператор $\Lambda_{xx}^{n+1/2}$ входит в обе подсхемы (в ② — явно, с уже найденного полуслоя), благодаря чему при сложении полушагов восстанавливается полный $6t^{n+1/2}\Lambda_{xx}^{n+1/2}$, а вся правая часть центрируется на $t^{n+1/2}$. Коэффициент $\sigma=6t$ берётся в обеих подсхемах на одном временно́м центре $t^{n+1/2}$.

Корректное расщепление конвекции. Конвективный член $C_x u_x$ — оператор по $x$, поэтому он относится только к подсхеме ① (неявной по $x$) и не повторяется в ②: иначе при суммировании полушагов конвекция задвоится ($-8u_x$ вместо $-4u_x$). Свободный член $f=txy$ берётся один раз (отнесём его к ②, на уровне $t^{n+1/2}$).

Конечная разность для конвекции — против потока (upwind), как в курсе (семинар 7). Конвекция в этом курсе всегда аппроксимируется односторонней разностью, выбранной так, чтобы усиливать диагональ (обеспечить безусловное диагональное преобладание). Для члена $+4\,u_x$ (он же $-4u_x$ в правой части) устойчивое направление — разность назад:

$$+4\,\frac{\partial u}{\partial x}\;\to\;+4\,\frac{u_{j,k}-u_{j-1,k}}{h_x}\qquad(\text{порядок }O(h_x)).$$

Вторые производные — центральными разностями второго порядка:

$$\Lambda_{xx}u_{j,k}=\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^2},\qquad \Lambda_{yy}u_{j,k}=\frac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^2}.$$

Подсхема ① ($n\to n+1/2$): неявно по $x$, явно по $y$

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}+4\,\frac{u_{j,k}^{n+1/2}-u_{j-1,k}^{n+1/2}}{h_x}=3t^{n+1/2}\,\Lambda_{xx}u_{j,k}^{n+1/2}+3t^{n+1/2}\,\Lambda_{yy}u_{j,k}^{n}.$$

Неявны по $x$ оператор $\Lambda_{xx}$ и вся конвекция $C_x u_x$ (на слое $n{+}1/2$); диффузия по $y$ — явно на слое $n$.

Подсхема ② ($n+1/2\to n+1$): неявно по $y$, явно по $x$

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}=3t^{n+1/2}\,\Lambda_{xx}u_{j,k}^{n+1/2}+3t^{n+1/2}\,\Lambda_{yy}u_{j,k}^{n+1}+t^{n+1/2}x_jy_k.$$

Оператор, неявный по $y$, на слое $n{+}1$; диффузия по $x$ — явно с полуслоя $n{+}1/2$; конвекции в ② нет (учтена целиком в ①); источник входит один раз на уровне $t^{n+1/2}$.

Контроль суммирования полушагов

Сложив ① и ② (оба со знаменателем $\Delta t$):

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=6t^{n+1/2}\Lambda_{xx}u^{n+1/2}-4\frac{u_{j,k}^{n+1/2}-u_{j-1,k}^{n+1/2}}{h_x}+3t^{n+1/2}\big(\Lambda_{yy}u^{n}+\Lambda_{yy}u^{n+1}\big)+t^{n+1/2}x_jy_k.$$

Диффузия по $x$ собралась в полный $6t$; диффузия по $y$ — полусумма уровней $n$ и $n{+}1$ (центр $t^{n+1/2}$); конвекция вошла ровно один раз с коэффициентом $-4$; источник учтён один раз. Это аппроксимация исходного уравнения относительно $t^{n+1/2}$.

Приведение к прогонке. Подсхема ① (прогонка по $j$)

Умножаем ① на $\Delta t$ и собираем слой $n{+}1/2$ слева: $a_j u_{j-1,k}^{n+1/2}+b_j u_{j,k}^{n+1/2}+c_j u_{j+1,k}^{n+1/2}=\xi_{j,k}$.

Знаки конвекции (upwind назад). Член $+4(u_j-u_{j-1})/h_x$ после переноса налево добавляет $+\dfrac{4\Delta t}{h_x}$ к диагонали $b_j$ и $-\dfrac{4\Delta t}{h_x}$ к $a_j$ (коэффициенту при $u_{j-1}$); на $c_j$ конвекция не влияет. Диффузия даёт $-\dfrac{3t^{n+1/2}\Delta t}{h_x^2}$ к обоим внедиагональным и $+\dfrac{6t^{n+1/2}\Delta t}{h_x^2}$ к диагонали. Итого (проверено символьно):

$$a_j=-\frac{3t^{n+1/2}\Delta t}{h_x^2}-\frac{4\Delta t}{h_x},\qquad b_j=1+\frac{6t^{n+1/2}\Delta t}{h_x^2}+\frac{4\Delta t}{h_x},\qquad c_j=-\frac{3t^{n+1/2}\Delta t}{h_x^2},$$

$$\xi_{j,k}=u_{j,k}^{n}+3t^{n+1/2}\Delta t\,\Lambda_{yy}u_{j,k}^{n}.$$

Сходимость прогонки (подсхема ①)

Введём $D=\dfrac{3t^{n+1/2}\Delta t}{h_x^2}\ge0$, $V=\dfrac{4\Delta t}{h_x}\ge0$. Тогда $a_j=-(D+V)\le0$, $c_j=-D\le0$, $b_j=1+2D+V$. Оба внедиагональных коэффициента неположительны при любых $h_x,\Delta t,t\ge0$ (M-матрица), и

$$|a_j|+|c_j|=(D+V)+D=2D+V<1+2D+V=|b_j|.$$

Таким образом, благодаря разности против потока строгое диагональное преобладание выполняется безусловно (никакого ограничения вида $h_x\le\tfrac32 t$ не требуется). В этом и состоит преимущество upwind перед центральной разностью для конвективного члена.

Рекуррентное прогоночное соотношение и $\alpha_1,\beta_1$ (①)

Прямой ход: $u_{j,k}^{n+1/2}=\alpha_{j+1}\,u_{j+1,k}^{n+1/2}+\beta_{j+1}$, где

$$\alpha_{j+1}=\frac{-c_j}{b_j+a_j\alpha_j},\qquad \beta_{j+1}=\frac{\xi_{j,k}-a_j\beta_j}{b_j+a_j\alpha_j}.$$

Левое ГУ по $x$ — 1-го рода: $u(t,0,y)=t\Rightarrow u_{1,k}^{n+1/2}=t^{n+1/2}$, откуда (чтобы $u_{1,k}^{n+1/2}=\alpha_1 u_{2,k}^{n+1/2}+\beta_1$ давало фиксированное значение):

$$\boxed{\alpha_1=0,\qquad \beta_1=t^{n+1/2}.}$$

Правое ГУ по $x$ — 1-го рода: $u_{N_x,k}^{n+1/2}=u(t,1,y)=2t^{n+1/2}$ — стартовое значение обратного хода.

Приведение к прогонке. Подсхема ② (прогонка по $k$)

В ② неявен только $\Lambda_{yy}$, конвекции нет ($C_y=0$, а $C_x u_x$ учтён в ①). Умножая на $\Delta t$ и собирая слой $n{+}1$ слева, $\tilde a_k u_{j,k-1}^{n+1}+\tilde b_k u_{j,k}^{n+1}+\tilde c_k u_{j,k+1}^{n+1}=\tilde\xi_{j,k}$:

$$\tilde a_k=\tilde c_k=-\frac{3t^{n+1/2}\Delta t}{h_y^2},\qquad \tilde b_k=1+\frac{6t^{n+1/2}\Delta t}{h_y^2},$$

$$\tilde\xi_{j,k}=u_{j,k}^{n+1/2}+3t^{n+1/2}\Delta t\,\Lambda_{xx}u_{j,k}^{n+1/2}+t^{n+1/2}x_jy_k\,\Delta t.$$

Конвективного слагаемого в $\tilde\xi$ нет; все слагаемые — с известного полуслоя $n{+}1/2$.

Сходимость прогонки (②) и рекуррентное соотношение

Внедиагональные коэффициенты строго отрицательны, конвекции нет, поэтому

$$|\tilde a_k|+|\tilde c_k|=\frac{6t^{n+1/2}\Delta t}{h_y^2}<1+\frac{6t^{n+1/2}\Delta t}{h_y^2}=|\tilde b_k|$$

выполнено безусловно (при $t\ge0$). Прямой ход: $u_{j,k}^{n+1}=\hat\alpha_{k+1}u_{j,k+1}^{n+1}+\hat\beta_{k+1}$, где

$$\hat\alpha_{k+1}=\frac{-\tilde c_k}{\tilde b_k+\tilde a_k\hat\alpha_k},\qquad \hat\beta_{k+1}=\frac{\tilde\xi_{j,k}-\tilde a_k\hat\beta_k}{\tilde b_k+\tilde a_k\hat\alpha_k}.$$

Левое ГУ по $y$ — 1-го рода: $u(t,x,0)=t\Rightarrow u_{j,1}^{n+1}=t^{n+1}$, откуда

$$\boxed{\hat\alpha_1=0,\qquad \hat\beta_1=t^{n+1}.}$$

Правое ГУ по $y$ — 1-го рода: $u_{j,N_y}^{n+1}=u(t,x,1)=2t^{n+1}$ — стартовое значение обратного хода.

Аппроксимация и алгоритм

Диффузионные члены — центральными разностями $O(h_x^2,h_y^2)$; конвективный член — разностью против потока $O(h_x)$, поэтому по пространству схема имеет порядок $O(h_x,h_y^2)$ (первый порядок по $x$ из-за upwind-конвекции — ср. формулировку семинара 8: $O(\Delta t,h_x^2,h_y)$ для уравнения с конвекцией по $y$). По времени диффузионная (симметрично расщеплённая) часть центрирована на $t^{n+1/2}$ и даёт $O(\Delta t^2)$, но конвекция входит лишь в один полушаг (несимметрично по структуре расщепления), поэтому в общем случае порядок по времени $O(\Delta t)$. Итого $O(\Delta t,\,h_x,\,h_y^2)$. По времени схема абсолютно устойчива (неявность по одному направлению в каждом полушаге); прогонки сходятся безусловно. Алгоритм одного шага:

1. $u^0=0$ (начальное условие $u(t{=}0,x,y)=0$).
2. Подсхема ①: для каждого $k$ — прогонка по $x$ ($\alpha_1=0,\ \beta_1=t^{n+1/2}$ слева; $u_{N_x,k}^{n+1/2}=2t^{n+1/2}$ справа) $\to u^{n+1/2}$.
3. Подсхема ②: для каждого $j$ — прогонка по $y$ ($\hat\alpha_1=0,\ \hat\beta_1=t^{n+1}$ слева; $u_{j,N_y}^{n+1}=2t^{n+1}$ справа) $\to u^{n+1}$.
4. Повторять по $n$ до $t_k$.

Уравнение

$$\frac{\partial u}{\partial t}+0{,}2\,\frac{\partial u}{\partial x}=0{,}1\,\frac{\partial u}{\partial y}+\cos(xy),\qquad u(t{=}0,x,y)=0.$$

1. Тип и канонический вид

В уравнении присутствуют только первые производные по пространству — это двумерное уравнение переноса (1-го порядка). Приводим к стандартному виду $\dfrac{\partial u}{\partial t}+v_1\dfrac{\partial u}{\partial x}+v_2\dfrac{\partial u}{\partial y}=f$, перенося $y$-слагаемое влево:

$$\frac{\partial u}{\partial t}+0{,}2\,\frac{\partial u}{\partial x}-0{,}1\,\frac{\partial u}{\partial y}=\cos(xy)\quad\Rightarrow\quad v_1=0{,}2,\;\; v_2=-0{,}1,\;\; f=\cos(xy).$$

2. Выбор граничных условий

Для уравнения 1-го порядка по каждому пространственному направлению нужно ровно одно граничное условие — на входной (по потоку) границе. Знак скорости определяет, откуда втекает информация, и сторону разности («против потока», upwind):

• по $x$: $v_1=0{,}2>0$ → поток идёт слева направо → входная граница $x{=}0$, берём ГУ $u(t,0,y)=\cos y$; разность левая (назад) $\dfrac{u_{j,k}-u_{j-1,k}}{h_x}$;
• по $y$: $v_2=-0{,}1<0$ → поток идёт сверху вниз → входная граница $y{=}1$, берём ГУ $u(t,x,1)=\sin x$; разность правая (вперёд) $\dfrac{u_{j,k+1}-u_{j,k}}{h_y}$.

Условия на выходных границах $u(t,1,y)=\sin y$ (при $x{=}1$) и $u(t,x,0)=\cos x$ (при $y{=}0$) для построения схемы не используются — иначе задача оказалась бы переопределённой.

3. Сетка

$x_j=j\,h_x,\ j=0,\dots,N_x$; $y_k=k\,h_y,\ k=0,\dots,N_y$; $t^n=n\,\Delta t$. Введём числа Куранта

$$r_x=\frac{0{,}2\,\Delta t}{h_x},\qquad r_y=\frac{0{,}1\,\Delta t}{h_y}.$$

4. Неявная схема методом дробных шагов

Шаг $\Delta t$ дробится на два полушага; на каждом неявен ровно один оператор переноса, второй «выключен». Источник $f=\cos(xy)$ относим целиком к первому полушагу. Сумма приращений двух полушагов воспроизводит полный оператор $\partial u/\partial t=-0{,}2\,\partial u/\partial x+0{,}1\,\partial u/\partial y+f$.

Подсхема ① ($n\to n+1/2$): неявно по $x$, явно (нулевой) по $y$, с источником. Аппроксимируем $\partial u/\partial x$ левой разностью (по потоку $v_1>0$):

$$\frac{u_{j,k}^{\,n+1/2}-u_{j,k}^{\,n}}{\Delta t}=-0{,}2\,\frac{u_{j,k}^{\,n+1/2}-u_{j-1,k}^{\,n+1/2}}{h_x}+\cos(x_j y_k).$$

Подсхема ② ($n+1/2\to n+1$): неявно по $y$. Аппроксимируем $\partial u/\partial y$ правой разностью (по потоку $v_2<0$); знак коэффициента $+0{,}1$ из исходной записи $\partial u/\partial t=+0{,}1\,\partial u/\partial y$:

$$\frac{u_{j,k}^{\,n+1}-u_{j,k}^{\,n+1/2}}{\Delta t}=0{,}1\,\frac{u_{j,k+1}^{\,n+1}-u_{j,k}^{\,n+1}}{h_y}.$$

5. Рекуррентные соотношения

Каждая подсхема двухдиагональна (входят лишь два соседних узла), поэтому решается прямой подстановкой (бегущим счётом) без полной прогонки. Обе матрицы строго диагонально преобладающие: $|1+r_x|>|r_x|$ и $|1+r_y|>|r_y|$.

Подсхема ①. Умножая на $\Delta t$ и группируя неизвестные слоя $n+1/2$:

$$(1+r_x)\,u_{j,k}^{\,n+1/2}-r_x\,u_{j-1,k}^{\,n+1/2}=u_{j,k}^{\,n}+\Delta t\,\cos(x_j y_k).$$

Идём от входной границы $x{=}0$ ($j=0$, где $u_{0,k}^{\,n+1/2}=\cos y_k$ известно) в сторону возрастания $j$:

$$\boxed{\,u_{j,k}^{\,n+1/2}=\frac{u_{j,k}^{\,n}+r_x\,u_{j-1,k}^{\,n+1/2}+\Delta t\,\cos(x_j y_k)}{1+r_x}\,},\qquad j=1,\dots,N_x.$$

Подсхема ②. Умножая на $\Delta t$ и группируя неизвестные слоя $n+1$:

$$(1+r_y)\,u_{j,k}^{\,n+1}-r_y\,u_{j,k+1}^{\,n+1}=u_{j,k}^{\,n+1/2}.$$

Идём от входной границы $y{=}1$ ($k=N_y$, где $u_{j,N_y}^{\,n+1}=\sin x_j$ известно) в сторону убывания $k$:

$$\boxed{\,u_{j,k}^{\,n+1}=\frac{u_{j,k}^{\,n+1/2}+r_y\,u_{j,k+1}^{\,n+1}}{1+r_y}\,},\qquad k=N_y-1,\dots,0.$$

6. Об устойчивости

Оба знаменателя $1+r_x>0$ и $1+r_y>0$ при любых $\Delta t,h_x,h_y$, а коэффициенты при «известных» узлах неотрицательны. Следовательно, неявная схема дробных шагов безусловно устойчива (ограничения на шаг по времени, в отличие от явной схемы, нет). Порядок аппроксимации — $O(\Delta t+h_x+h_y)$ (по пространству первый из-за односторонних upwind-разностей).

7. Порядок счёта (за один шаг $\Delta t$)

1. При $n=0$: $u^{0}_{j,k}=0$ во всех внутренних узлах.
2. Подсхема ①: для каждой строки $k$ пробегаем $j=1,\dots,N_x$ по формуле первого бокса, стартуя от $u_{0,k}^{\,n+1/2}=\cos y_k$ — получаем промежуточный слой $u^{\,n+1/2}$.
3. Подсхема ②: для каждого столбца $j$ пробегаем $k=N_y-1,\dots,0$ по формуле второго бокса, стартуя от $u_{j,N_y}^{\,n+1}=\sin x_j$ — получаем слой $u^{\,n+1}$.
4. Переходим к следующему $n$.

Условие

Привести к виду, удобному для метода простой итерации, краевую задачу

$$2\,\frac{du}{dx}+\frac{d^2u}{dx^2}+3x^2=0,\qquad u(x{=}0)=0,\quad u(x{=}1)=1.$$

Записать шаг итерации, итерационное соотношение, условие окончания итераций и начальное приближение.

1. Подготовка уравнения

Перепишем уравнение, расположив старшую производную первой:

$$\frac{d^2u}{dx^2}+2\,\frac{du}{dx}+3x^2=0.$$

Это линейное ОДУ 2-го порядка с граничными условиями первого рода (заданы значения функции на обоих концах). Идея метода простой итерации: записать разностную аппроксимацию, разрешить её относительно центрального узла $u_j$ и организовать последовательность приближений $u_j^{0}\to u_j^{1}\to\dots$, на каждом шаге пересчитывая $u_j$ через соседей, взятых с предыдущей итерации (поточечная итерация Якоби).

2. Сетка и разностная аппроксимация

Введём равномерную сетку $x_j=(j-1)h$, $j=1,\dots,N_x$, шаг $h=\dfrac{1}{N_x-1}$, так что $x_1=0$, $x_{N_x}=1$. Используем стандартные операторы курса со вторым порядком по $x$:

$$\frac{d^2u}{dx^2}\to\frac{u_{j+1}-2u_j+u_{j-1}}{h^2},\qquad \frac{du}{dx}\to\frac{u_{j+1}-u_{j-1}}{2h}$$

(центральная разность для первой производной допустима, так как член диффузии $u''$ обеспечивает диагональное доминирование). Заметим, что в самой разности первой производной узел $u_{j-1}$ входит со знаком минус; результирующий же коэффициент при $u_{j-1}$ после сложения с членом $u''$ станет $+(1-h)$ (см. п.3). Подставляя во внутренних узлах $j=2,\dots,N_x-1$:

$$\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}+2\cdot\frac{u_{j+1}-u_{j-1}}{2h}+3x_j^2=0.$$

3. Приведение к виду для прогонки (трёхдиагональная система)

Умножим на $h^2$ и соберём по узлам — получаем систему вида $c_j u_{j-1}+b_j u_j+a_j u_{j+1}=\xi_j$:

$$(1-h)\,u_{j-1}-2\,u_j+(1+h)\,u_{j+1}=-3h^2x_j^2,$$ $$\boxed{\;c_j=1-h,\qquad b_j=-2,\qquad a_j=1+h,\qquad \xi_j=-3h^2x_j^2.\;}$$

Достаточное условие сходимости прогонки (диагональное преобладание) $|a_j|+|c_j|\le|b_j|$ при $0 $$|a_j|+|c_j|=(1+h)+(1-h)=2=|b_j|.$$

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

4. Итерационное соотношение (шаг простой итерации)

Разрешим разностное уравнение относительно диагонального члена $u_j$ (коэффициент при $u_j$ равен $-2$):

$$2u_j=(1-h)\,u_{j-1}+(1+h)\,u_{j+1}+3h^2x_j^2.$$

Соседние значения берём с предыдущей итерации $s$, искомое — на новой итерации $s+1$ (верхний индекс — номер итерации):

$$\boxed{\;u_j^{\,s+1}=\frac{1+h}{2}\,u_{j+1}^{\,s}+\frac{1-h}{2}\,u_{j-1}^{\,s}+\frac{3h^2}{2}\,x_j^2,\qquad j=2,\dots,N_x-1.\;}$$

Это и есть выражение для шага итерации: один проход по всем внутренним узлам формирует новый слой приближения $u^{s+1}$ из старого $u^{s}$. Граничные узлы фиксированы условиями первого рода и на каждой итерации не меняются:

$$u_1^{\,s+1}=0,\qquad u_{N_x}^{\,s+1}=1.$$

5. Сходимость простой итерации

Матрица перехода имеет в каждой строке только два ненулевых внедиагональных коэффициента $\dfrac{1+h}{2}$ и $\dfrac{1-h}{2}$. Норма строки (достаточное условие сходимости $\|C\|_\infty\le1$):

$$\frac{|1+h|}{2}+\frac{|1-h|}{2}=\frac{(1+h)+(1-h)}{2}=1\qquad(0Сумма модулей коэффициентов равна $1$ — итерация принадлежит «граничному» (как у уравнения Лапласа) случаю; благодаря закреплённым значениям на границах (условия 1-го рода) спектральный радиус матрицы перехода строго меньше единицы, и процесс сходится к решению разностной задачи. (Численная проверка при $N_x=21$: спектральный радиус $\rho\approx0{,}986<1$; при точности $\varepsilon=10^{-8}$ процесс сходится примерно за $\sim10^{3}$ итераций — при $\varepsilon=10^{-4}$ заметно быстрее; расхождение с точным решением ОДУ $\sim5\cdot10^{-5}$, что отвечает порядку $O(h^2)$.)

6. Аналог прогоночных коэффициентов на левой границе ($\alpha_1,\beta_1$)

Если на каждой итерации вместо поточечного пересчёта применять прогонку для системы из п.3, прогоночное соотношение записывается как $u_j=\alpha_j u_{j+1}+\beta_j$ с

$$\alpha_j=\frac{-a_j}{b_j+c_j\alpha_{j-1}},\qquad \beta_j=\frac{\xi_j-c_j\beta_{j-1}}{b_j+c_j\alpha_{j-1}}.$$

Левое граничное условие 1-го рода $u_1=0$ означает $u_1=\alpha_1 u_2+\beta_1$ при любых $u_2$, откуда

$$\boxed{\;\alpha_1=0,\qquad \beta_1=0.\;}$$

(Аналогично правое условие $u_{N_x}=1$ даёт замыкание прогонки.) Для самого метода простой итерации эти коэффициенты не нужны — пересчёт идёт по формуле п.4.

7. Начальное приближение

Берём гладкую функцию, удовлетворяющую граничным условиям, — линейную интерполяцию между значениями на концах $u(0)=0$, $u(1)=1$:

$$\boxed{\;u_j^{\,0}=x_j\quad(\text{то есть }u_j^{0}=(j-1)h),\qquad u_1^{0}=0,\;u_{N_x}^{0}=1.\;}$$

Можно взять и $u_j^0\equiv0$ во внутренних узлах — сходимость от выбора старта не зависит, линейный профиль лишь ускоряет её.

8. Условие окончания итерационного процесса

Итерации прекращают, когда соседние приближения перестают различаться, — по норме разности:

$$\boxed{\;\max_{2\le j\le N_x-1}\big|u_j^{\,s+1}-u_j^{\,s}\big|<\varepsilon\;}$$

с заданной точностью $\varepsilon$ (например $\varepsilon=10^{-4}$). После выполнения неравенства массив $u_j^{\,s+1}$ принимается за численное решение краевой задачи.

Уравнение и условия

Дано трёхмерное (по координатам $x,y,z$) уравнение параболического типа с конвективными членами:

$$\frac{\partial u}{\partial t}-2\,\frac{\partial u}{\partial x}+3\,\frac{\partial u}{\partial z}=5\,\frac{\partial^2 u}{\partial y^2}+6\,\frac{\partial^2 u}{\partial z^2}+txyz,\qquad u=u(t,x,y,z),$$

на единичном кубе $x,y,z\in[0,1]$, $t\in[0,t_k]$, с граничными и начальным условиями

$$u(t,1,y,z)=yz,\qquad \begin{cases}u(t,x,0,z)=t,\\ u(t,x,1,z)=zx,\end{cases}\qquad \begin{cases}u(t,x,y,0)=t,\\ u(t,x,y,1)=xy,\end{cases}\qquad u(0,x,y,z)=xyz.$$

1. Канонический вид и выделение операторов (гл. 9, §10.1, семинар 9)

Приводим к стандартной форме нестационарной задачи с конвекцией (учебник, формула (9.22)):

$$\frac{\partial u}{\partial t}+v_x\,\frac{\partial u}{\partial x}+v_z\,\frac{\partial u}{\partial z}=\sigma_y\,\frac{\partial^2 u}{\partial y^2}+\sigma_z\,\frac{\partial^2 u}{\partial z^2}+f,$$ $$v_x=-2,\qquad v_z=+3,\qquad \sigma_y=5,\qquad \sigma_z=6,\qquad f=txyz.$$

Вторые производные присутствуют только по $y$ и по $z$; первые производные (конвекция) — по $x$ и по $z$. Это частный случай уравнения (9.25): часть производных отсутствует.

Стандартные операторы курса (вторая разность — знак $+$ у соседних узлов):

$$\Lambda_{yy}u_{j,k,m}=\frac{u_{j,k+1,m}-2u_{j,k,m}+u_{j,k-1,m}}{h_y^2}=O(h_y^2),\qquad \Lambda_{zz}u_{j,k,m}=\frac{u_{j,k,m+1}-2u_{j,k,m}+u_{j,k,m-1}}{h_z^2}=O(h_z^2).$$

2. Правило выбора односторонних разностей для конвекции (§10.1, семинар 7)

Конвективные члены стоят в левой части. Правило по знаку скорости: при $v>0$ — левая разность $\dfrac{u_i-u_{i-1}}{h}$, при $v<0$ — правая разность $\dfrac{u_{i+1}-u_i}{h}$. Обе односторонние, поэтому первого порядка по своему шагу.

• По $x$: $v_x=-2<0\ \Rightarrow$ правая разность $\dfrac{u_{j+1,k,m}-u_{j,k,m}}{h_x}$, $O(h_x)$. Единственное заданное условие по $x$ — правое $u(t,1,y,z)=yz$ (вход для счёта справа налево). Согласовано с правилом.
• По $z$ (конвекция): $v_z=+3>0\ \Rightarrow$ левая разность $\dfrac{u_{j,k,m}-u_{j,k,m-1}}{h_z}$, $O(h_z)$.

3. Структура схемы предиктор-корректор (гл. 9, §8.1)

В трёхмерной задаче схема предиктор-корректор строится так (учебник, (9.17)–(9.21)): интервал $\delta t=[t^n,t^{n+1}]$ делится пополам точкой $t^{n+1/2}$; первый полуинтервал $\Delta t/2$ расщепляется на равные части по числу пространственных направлений, и на каждой части записывается своя неявная подсхема предиктора, учитывающая операторы только одного направления. Совокупность подсхем-предикторов даёт значения на слое $t^{n+1/2}$. Затем корректор записывается на полном шаге $\Delta t$ как поправочное соотношение $\dfrac{u^{n+1}-u^{n}}{\Delta t}=(\text{все операторы при }t^{n+1/2})$; так как правая часть берётся при $t^{n+1/2}$, разность $\dfrac{u^{n+1}-u^{n}}{\Delta t}$ оказывается центральной относительно $t^{n+1/2}$ и даёт второй порядок по $t$. Всего схема состоит из четырёх подсхем: три предиктора + один корректор.

Здесь направлений три ($x$, $y$, $z$), поэтому $\Delta t/2$ делится на три части: $t^{n}\to t^{n+1/6}\to t^{n+1/3}\to t^{n+1/2}$. Каждая подсхема предиктора берёт операторы своего направления неявно (для направления $z$ это и диффузия, и конвекция — как в (9.24)); направление $x$ без второй производной даёт подсхему-аналог УрЧП первого порядка (как третья подсхема в (9.26)). Источник $f$ и «сквозные» члены вносятся в корректоре.

Предиктор (полушаг $\Delta t/2$, три подсхемы)

① по $x$ ($n\to n+1/6$): только конвекция по $x$ (правая разность, $v_x=-2$):

$$\frac{u_{j,k,m}^{\,n+1/6}-u_{j,k,m}^{\,n}}{\Delta t/2}-2\,\frac{u_{j+1,k,m}^{\,n+1/6}-u_{j,k,m}^{\,n+1/6}}{h_x}=0.$$

② по $y$ ($n+1/6\to n+1/3$): только диффузия по $y$ (неявно):

$$\frac{u_{j,k,m}^{\,n+1/3}-u_{j,k,m}^{\,n+1/6}}{\Delta t/2}=5\,\Lambda_{yy}\,u_{j,k,m}^{\,n+1/3}.$$

③ по $z$ ($n+1/3\to n+1/2$): диффузия и конвекция по $z$ (неявно, левая разность):

$$\frac{u_{j,k,m}^{\,n+1/2}-u_{j,k,m}^{\,n+1/3}}{\Delta t/2}+3\,\frac{u_{j,k,m}^{\,n+1/2}-u_{j,k,m-1}^{\,n+1/2}}{h_z}=6\,\Lambda_{zz}\,u_{j,k,m}^{\,n+1/2}.$$

Корректор (полный шаг $\Delta t$, все операторы при $t^{n+1/2}$)

Корректор записывается как поправочное соотношение (учебник, (9.20)): слева — разность между новым слоем $u^{n+1}$ и исходным слоем $u^{n}$ на полном шаге $\Delta t$; справа все операторы аппроксимированы относительно $t^{n+1/2}$. Именно такая запись делает левую разность центральной относительно $t^{n+1/2}$ (второй порядок по $t$):

$$\frac{u_{j,k,m}^{\,n+1}-u_{j,k,m}^{\,n}}{\Delta t}=5\,\Lambda_{yy}u_{j,k,m}^{\,n+1/2}+6\,\Lambda_{zz}u_{j,k,m}^{\,n+1/2}+2\,\frac{u_{j+1,k,m}^{\,n+1/2}-u_{j,k,m}^{\,n+1/2}}{h_x}-3\,\frac{u_{j,k,m}^{\,n+1/2}-u_{j,k,m-1}^{\,n+1/2}}{h_z}+f_{j,k,m}^{\,n+1/2},$$

где знаки конвекции получены переносом членов $-2\,\partial u/\partial x$ и $+3\,\partial u/\partial z$ из левой части в правую: $-(-2)=+2$ и $-(+3)=-3$; источник на полушаге $f_{j,k,m}^{\,n+1/2}=t_{n+1/2}\,x_j y_k z_m$, $\ t_{n+1/2}=t_n+\dfrac{\Delta t}{2}$.

4. Рекуррентные соотношения для подсхем

Подсхема ① (направление $x$). Вторых производных по $x$ нет — это аналог неявной схемы для УрЧП первого порядка, решается рекуррентным соотношением (а не прогонкой). Собирая члены:

$$u_{j,k,m}^{\,n+1/6}\Big(1+\frac{\Delta t}{h_x}\Big)-\frac{\Delta t}{h_x}\,u_{j+1,k,m}^{\,n+1/6}=u_{j,k,m}^{\,n}\ \Longrightarrow\ \boxed{\,u_{j,k,m}^{\,n+1/6}=\dfrac{u_{j,k,m}^{\,n}+\dfrac{\Delta t}{h_x}\,u_{j+1,k,m}^{\,n+1/6}}{1+\dfrac{\Delta t}{h_x}}\,}$$

(использовано $|v_x|\dfrac{\Delta t/2}{h_x}=\dfrac{2\cdot\Delta t/2}{h_x}=\dfrac{\Delta t}{h_x}$). Так как $v_x<0$, счёт идёт справа налево; правый узел $j=N_x$ задаётся условием $u(t,1,y,z)=yz$.

Подсхема ② (направление $y$, прогонка). Раскрывая $\Lambda_{yy}$ и приводя к виду $a_k u_{k-1}+b_k u_k+c_k u_{k+1}=\xi_k$:

$$a_k=c_k=-\frac{5\,\Delta t}{2h_y^2},\qquad b_k=1+\frac{5\,\Delta t}{h_y^2},\qquad \xi_{j,k,m}=u_{j,k,m}^{\,n+1/6}.$$

Прямой ход прогонки: $u_{j,k,m}^{\,n+1/3}=\alpha_{k+1}u_{j,k+1,m}^{\,n+1/3}+\beta_{k+1}$, где

$$\alpha_{k+1}=\frac{-c_k}{b_k+a_k\alpha_k},\qquad \beta_{k+1}=\frac{\xi_k-a_k\beta_k}{b_k+a_k\alpha_k}.$$

Граничные условия I рода по $y$: $u(t,x,0,z)=t\Rightarrow \alpha_1=0,\ \beta_1=t_{n+1/2}$; правая граница $u_{j,N_y,m}^{\,n+1/3}=z_m x_j$ (то есть $u(t,x,1,z)=zx$).

Подсхема ③ (направление $z$, прогонка). Неявны и диффузия $6\Lambda_{zz}$, и левая разность конвекции $+3\,\tfrac{u_m-u_{m-1}}{h_z}$. Перенося всё на слой $n+1/2$ и умножая на $\Delta t/2$, получаем $\tilde a_m u_{m-1}+\tilde b_m u_m+\tilde c_m u_{m+1}=\tilde\xi_m$ с

$$\tilde a_m=-\frac{\Delta t}{2}\Big(\frac{6}{h_z^2}+\frac{3}{h_z}\Big),\qquad \tilde c_m=-\frac{6}{h_z^2}\cdot\frac{\Delta t}{2}=-\frac{3\,\Delta t}{h_z^2},$$ $$\tilde b_m=1+\frac{\Delta t}{2}\Big(\frac{12}{h_z^2}+\frac{3}{h_z}\Big)=1+\frac{6\,\Delta t}{h_z^2}+\frac{3\,\Delta t}{2h_z},\qquad \tilde\xi_{j,k,m}=u_{j,k,m}^{\,n+1/3}.$$

Прямой ход прогонки: $u_{j,k,m}^{\,n+1/2}=\tilde\alpha_{m+1}u_{j,k,m+1}^{\,n+1/2}+\tilde\beta_{m+1}$, $\ \tilde\alpha_{m+1}=\dfrac{-\tilde c_m}{\tilde b_m+\tilde a_m\tilde\alpha_m}$, $\ \tilde\beta_{m+1}=\dfrac{\tilde\xi_m-\tilde a_m\tilde\beta_m}{\tilde b_m+\tilde a_m\tilde\alpha_m}$.

Граничные условия I рода по $z$: $u(t,x,y,0)=t\Rightarrow \tilde\alpha_1=0,\ \tilde\beta_1=t_{n+1/2}$; правая граница $u_{j,k,N_z}^{\,n+1/2}=x_j y_k$ (то есть $u(t,x,y,1)=xy$).

Корректор (рекуррентное соотношение). Корректор явный относительно слоя $n+1$ — пересчёт по исходному слою $n$ с пространственными операторами, взятыми по уже найденному слою $n+1/2$:

$$\boxed{\,u_{j,k,m}^{\,n+1}=u_{j,k,m}^{\,n}+\Delta t\Big[5\,\Lambda_{yy}u^{\,n+1/2}+6\,\Lambda_{zz}u^{\,n+1/2}+2\,\frac{u_{j+1,k,m}^{\,n+1/2}-u_{j,k,m}^{\,n+1/2}}{h_x}-3\,\frac{u_{j,k,m}^{\,n+1/2}-u_{j,k,m-1}^{\,n+1/2}}{h_z}+f_{j,k,m}^{\,n+1/2}\Big].\,}$$

На границах берутся точные значения условий варианта на слое $t^{n+1}$.

5. Устойчивость прогонок

Подсхема ② ($y$): $|a_k|+|c_k|=\dfrac{5\Delta t}{h_y^2}<1+\dfrac{5\Delta t}{h_y^2}=|b_k|$ — диагональное преобладание выполнено безусловно.

Подсхема ③ ($z$): $|\tilde a_m|+|\tilde c_m|=\dfrac{6\Delta t}{h_z^2}+\dfrac{3\Delta t}{2h_z}<1+\dfrac{6\Delta t}{h_z^2}+\dfrac{3\Delta t}{2h_z}=|\tilde b_m|$ — также выполнено (разность $=1$).

Подсхема ① ($x$, рекуррентная) и корректор устойчивы безусловно. В целом схема предиктор-корректор абсолютно устойчива (роль предиктора — обеспечение устойчивости, роль корректора — повышение порядка по времени).

6. Порядок аппроксимации (§10.2, §8.2)

• По времени: корректор записан как $\dfrac{u^{n+1}-u^{n}}{\Delta t}$ с правой частью при $t^{n+1/2}$, то есть центрирован относительно $t^{n+1/2}$ $\Rightarrow$ второй порядок $O(\Delta t^2)$.
• По $y$: только диффузия $\Lambda_{yy}$ с симметричной второй разностью $\Rightarrow$ второй порядок $O(h_y^2)$.
• По $z$: диффузия $\Lambda_{zz}$ сама по себе $O(h_z^2)$, но добавленный к этому направлению конвективный член $+3\,\partial u/\partial z$ аппроксимирован односторонней (левой) разностью первого порядка $O(h_z)$. Суммарный порядок по $z$ определяется наименьшим $\Rightarrow$ первый, $O(h_z)$.
• По $x$: только конвекция, односторонняя (правая) разность $\Rightarrow$ первый порядок $O(h_x)$.

Итоговый порядок аппроксимации схемы:

$$\psi=O\big(\Delta t^2,\ h_y^2,\ h_z,\ h_x\big).$$

Блок-схема алгоритма — предиктор из трёх последовательных подсхем (① рекуррентно по $x$, ② прогонка по $y$, ③ прогонка по $z$) с получением слоя $t^{n+1/2}$, затем корректор-пересчёт на $t^{n+1}$ — см. страницу «Блок-схемы».

Уравнение и условия

Двумерное параболическое уравнение с переменными коэффициентами диффузии, конвекцией и нелинейным стоком:

$$\frac{\partial u}{\partial t}-\frac{\partial u}{\partial x}-\frac{\partial u}{\partial y}=y\,\frac{\partial^2 u}{\partial x^2}+x\,\frac{\partial^2 u}{\partial y^2}-3u^2.$$

Граничные и начальное условия (все 1-го рода, Дирихле):

$$u(t,0,y)=ty,\quad u(t,1,y)=0;\qquad u(t,x,0)=tx,\quad u(t,x,1)=2;\qquad u(0,x,y)=0.$$

Приведём к канонической форме главы 9 (9.22), оставив первые производные слева:

$$\frac{\partial u}{\partial t}+v_1\frac{\partial u}{\partial x}+v_2\frac{\partial u}{\partial y}=y\,\frac{\partial^2u}{\partial x^2}+x\,\frac{\partial^2u}{\partial y^2}-3u^2,\qquad v_1=v_2=-1.$$

Коэффициенты диффузии переменные: $\sigma_x=y_k$, $\sigma_y=x_j$ (берутся в текущем узле). Скорости конвекции $v_1=v_2=-1<0$ — это и есть знаки перед $u_x,u_y$ в форме (9.22); именно они определяют выбор односторонних разностей. Нелинейный сток $-3u^2$.

Сетка и разностные операторы

$u_{j,k}^n=u(t^n,x_j,y_k)$, шаги $h_x,h_y,\Delta t$; $x_j=jh_x$, $y_k=kh_y$, $t^n=n\Delta t$. Вторые разности (центральные, 2-й порядок):

$$\Lambda_{xx}u_{j,k}=\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^{2}},\qquad \Lambda_{yy}u_{j,k}=\frac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^{2}}.$$

Выбор разности для первых производных (правило канона, lek9_10). Конечную разность выбирают по знаку скорости так, чтобы обеспечить безусловное диагональное преобладание и устойчивость прогонки: при $v<0$ — правая (forward) разность, при $v>0$ — левая (backward). У нас $v_1=v_2=-1<0$, поэтому берём правые разности:

$$\frac{\partial u}{\partial x}\approx\frac{u_{j+1,k}-u_{j,k}}{h_x},\qquad\frac{\partial u}{\partial y}\approx\frac{u_{j,k+1}-u_{j,k}}{h_y}.$$

Это аппроксимация первого порядка по конвектируемым координатам.

Идея схемы расщепления (метод дробных шагов)

Шаг $\Delta t$ разбивается промежуточным слоем $t^{n+1/2}$ на две подсхемы. По канону схемы расщепления (метод дробных шагов, гл. 9, разд. 6) в каждой подсхеме присутствует только один направленный оператор (диффузия вместе со «своей» конвекцией), неявно и с полным коэффициентом диффузии $\sigma$ (без $\sigma/2$); оператор другого направления в подсхему не входит. Система трёхдиагональна (прогонка), а за полный шаг каждый оператор срабатывает ровно один раз.

Нелинейный сток $-3u^2$ линеаризуем по уже вычисленному слою и относим (один раз) во вторую подсхему: $-3u^2\approx-3\big(u_{j,k}^{\,n+1/2}\big)^2$.

Подсхема ① ($n\to n+1/2$): неявно по $x$ (прогонка по $j$)

$$\frac{u_{j,k}^{\,n+1/2}-u_{j,k}^{\,n}}{\Delta t}+v_1\frac{u_{j+1,k}^{\,n+1/2}-u_{j,k}^{\,n+1/2}}{h_x}=y_k\,\Lambda_{xx}u_{j,k}^{\,n+1/2}.$$

Содержит только оператор по $x$ (диффузия с полным коэффициентом $y_k$ + конвекция, правая разность при $v_1<0$) на слое $n+1/2$; правая часть — $u^{n}$.

Подсхема ② ($n+1/2\to n+1$): неявно по $y$ (прогонка по $k$)

$$\frac{u_{j,k}^{\,n+1}-u_{j,k}^{\,n+1/2}}{\Delta t}+v_2\frac{u_{j,k+1}^{\,n+1}-u_{j,k}^{\,n+1}}{h_y}=x_j\,\Lambda_{yy}u_{j,k}^{\,n+1}-3\big(u_{j,k}^{\,n+1/2}\big)^2.$$

Содержит только оператор по $y$ (полный коэффициент $x_j$, правая разность при $v_2<0$) на слое $n+1$; нелинейный сток взят на слое $n+1/2$. Оператор по $x$ сюда не входит.

Приведение подсхемы ① к трёхдиагональному виду (прогонка по $j$)

Умножим на $\Delta t$ и перенесём неизвестные слоя $n+1/2$ влево. С $v_1=-1$ правая разность даёт $-\dfrac{\Delta t}{h_x}(u_{j+1,k}-u_{j,k})$; конвекция вкладывает $|v_1|\Delta t/h_x$ в диагональ $b_j$ и в внедиагональ $c_j$ (узел $j+1$), а в $a_j$ — ноль (только диффузия):

$$a_j\,u_{j-1,k}^{\,n+1/2}+b_j\,u_{j,k}^{\,n+1/2}+c_j\,u_{j+1,k}^{\,n+1/2}=\xi_{j,k},$$ $$\boxed{\;a_j=-\frac{y_k\,\Delta t}{h_x^{2}},\qquad b_j=1+\frac{\Delta t}{h_x}+\frac{2\,y_k\,\Delta t}{h_x^{2}},\qquad c_j=-\frac{\Delta t}{h_x}-\frac{y_k\,\Delta t}{h_x^{2}},\qquad \xi_{j,k}=u_{j,k}^{\,n}.\;}$$

Контроль знаков конвекции: при $v_1=-1<0$ берётся правая (forward) разность — единственный устойчивый выбор: левая разность при $v<0$ испортила бы знак $a_j$ и сняла бы безусловное преобладание. Её вклад $|v_1|\Delta t/h_x=\Delta t/h_x$ идёт в диагональ $b_j$ и в $c_j$ (узел $j+1$), а $a_j$ конвекцией не затрагивается. Коэффициент диффузии полный ($y_k\Delta t/h_x^2$, без делителя 2). Совпадает с канонической диагональю $b_j=1+|v_1|\Delta t/h_x+2\sigma\Delta t/h_x^2$ (lek9_10_1).

Сходимость прогонки для подсхемы ①

Достаточное условие — диагональное преобладание $|a_j|+|c_j|\le|b_j|$. Во внутренних узлах $y_k>0$; оба внедиагональных коэффициента отрицательны при любом шаге (это и есть преимущество односторонней разности):

$$|a_j|=\frac{y_k\Delta t}{h_x^{2}},\qquad |c_j|=\frac{\Delta t}{h_x}+\frac{y_k\Delta t}{h_x^{2}},\qquad |a_j|+|c_j|=\frac{\Delta t}{h_x}+\frac{2\,y_k\,\Delta t}{h_x^{2}}.$$ $$|b_j|-(|a_j|+|c_j|)=\Big(1+\frac{\Delta t}{h_x}+\frac{2y_k\Delta t}{h_x^2}\Big)-\Big(\frac{\Delta t}{h_x}+\frac{2y_k\Delta t}{h_x^2}\Big)=1>0.$$

Преобладание строгое и безусловное (никаких ограничений вида $h_x<2y_k$ не требуется). Прогонка устойчива.

Рекуррентное соотношение и $\alpha_1,\beta_1$ (подсхема ①, по $x$)

Прогонка в форме учебника (4.11) $u_{j,k}^{\,n+1/2}=\alpha_j\,u_{j+1,k}^{\,n+1/2}+\beta_j$ с коэффициентами (4.13):

$$\alpha_j=\frac{-a_j}{\,b_j+c_j\,\alpha_{j-1}\,},\qquad \beta_j=\frac{\xi_{j,k}-c_j\,\beta_{j-1}}{\,b_j+c_j\,\alpha_{j-1}\,}.$$

Левая граница $x=0$ — ГУ 1-го рода $u(t,0,y)=t\,y$: $u_{1,k}^{\,n+1/2}=t^{\,n+1/2}y_k$. Записывая её в виде $u_{1,k}=\alpha_1u_{2,k}+\beta_1$:

$$\boxed{\;\alpha_1=0,\qquad \beta_1=t^{\,n+1/2}\,y_k.\;}$$

Правая граница $x=1$ — ГУ 1-го рода $u(t,1,y)=0$: $u_{N_x,k}^{\,n+1/2}=0$, отсюда стартует обратный ход.

Приведение подсхемы ② к трёхдиагональному виду (прогонка по $k$)

Аналогично, $\sigma_y=x_j$, $v_2=-1<0$ ⇒ правая разность по $y$, полный коэффициент диффузии:

$$\tilde a_k\,u_{j,k-1}^{\,n+1}+\tilde b_k\,u_{j,k}^{\,n+1}+\tilde c_k\,u_{j,k+1}^{\,n+1}=\tilde\xi_{j,k},$$ $$\boxed{\;\tilde a_k=-\frac{x_j\,\Delta t}{h_y^{2}},\qquad \tilde b_k=1+\frac{\Delta t}{h_y}+\frac{2\,x_j\,\Delta t}{h_y^{2}},\qquad \tilde c_k=-\frac{\Delta t}{h_y}-\frac{x_j\,\Delta t}{h_y^{2}},\;}$$ $$\tilde\xi_{j,k}=u_{j,k}^{\,n+1/2}-3\,\Delta t\,\big(u_{j,k}^{\,n+1/2}\big)^2.$$

Оператор по $x$ в правую часть не входит; нелинейный сток включён в $\tilde\xi$ через известный слой $n+1/2$.

Сходимость прогонки для подсхемы ②

При $x_j>0$ оба внедиагональных коэффициента отрицательны при любом шаге, и

$$|\tilde a_k|+|\tilde c_k|=\frac{\Delta t}{h_y}+\frac{2\,x_j\,\Delta t}{h_y^{2}},\qquad |\tilde b_k|-(|\tilde a_k|+|\tilde c_k|)=1>0.$$

Преобладание строгое и безусловное — прогонка устойчива.

Рекуррентное соотношение и $\hat\alpha_1,\hat\beta_1$ (подсхема ②, по $y$)

$u_{j,k}^{\,n+1}=\hat\alpha_k\,u_{j,k+1}^{\,n+1}+\hat\beta_k$,

$$\hat\alpha_k=\frac{-\tilde a_k}{\,\tilde b_k+\tilde c_k\,\hat\alpha_{k-1}\,},\qquad \hat\beta_k=\frac{\tilde\xi_{j,k}-\tilde c_k\,\hat\beta_{k-1}}{\,\tilde b_k+\tilde c_k\,\hat\alpha_{k-1}\,}.$$

Нижняя граница $y=0$ — ГУ 1-го рода $u(t,x,0)=t\,x$: $u_{j,1}^{\,n+1}=t^{\,n+1}x_j$, т.е.

$$\boxed{\;\hat\alpha_1=0,\qquad \hat\beta_1=t^{\,n+1}\,x_j.\;}$$

Верхняя граница $y=1$ — ГУ 1-го рода $u(t,x,1)=2$: $u_{j,N_y}^{\,n+1}=2$, отсюда обратный ход.

Порядок аппроксимации и алгоритм

Из-за наличия первых производных и канонической односторонней (против потока) разности схема имеет порядок $O(\Delta t,\;h_x,\;h_y)$: первый по времени (последовательное расщепление операторов) и первый по пространству по обеим конвектируемым координатам (вторые разности дали бы $h^2$, но односторонняя аппроксимация конвекции понижает порядок до $h$). Алгоритм одного шага $n\to n+1$:

1. Начальное условие $u^{0}\equiv0$.
2. Подсхема ①: для каждого $k$ — прогонка по $x$ со стартом $\alpha_1=0,\ \beta_1=t^{\,n+1/2}y_k$ и правой границей $u_{N_x,k}=0$ $\;\Rightarrow\;u^{\,n+1/2}$.
3. Подсхема ②: для каждого $j$ — прогонка по $y$ со стартом $\hat\alpha_1=0,\ \hat\beta_1=t^{\,n+1}x_j$ и верхней границей $u_{j,N_y}=2$; сток взят на $u^{\,n+1/2}$ $\;\Rightarrow\;u^{\,n+1}$.
4. Переход $n\to n+1$ и повтор.

Уравнение и данные

$$\dfrac{\partial u}{\partial t}-\dfrac{\partial u}{\partial x}=\dfrac{\partial u}{\partial y}-e^{txy},$$

с условиями

$$u(t,0,y)=1,\quad u(t,1,y)=e^{y};\qquad u(t,x,0)=1,\quad u(t,x,1)=e^{x};\qquad u(0,x,y)=1.$$

1. Тип уравнения и канонический вид

В уравнении присутствуют только первые производные по $t$, $x$, $y$ — это двумерное уравнение 1-го порядка (уравнение переноса). Приводим его к стандартному виду $\dfrac{\partial u}{\partial t}+v_1\dfrac{\partial u}{\partial x}+v_2\dfrac{\partial u}{\partial y}=f$, перенося все пространственные производные влево:

$$\dfrac{\partial u}{\partial t}-\dfrac{\partial u}{\partial x}-\dfrac{\partial u}{\partial y}=-\,e^{txy}\quad\Rightarrow\quad v_1=-1,\;\; v_2=-1,\;\; f(t,x,y)=-\,e^{txy}.$$

2. Выбор «соответствующих» граничных условий

Для уравнения 1-го порядка по каждому пространственному направлению ставится ровно одно граничное условие — на входной (по потоку) границе. Сторону входа определяет знак скорости переноса, а разностный оператор берётся «против потока» (upwind), то есть со стороны, откуда приходит информация.

• Направление $x$: $v_1=-1<0$ — характеристики направлены справа налево, поток входит через границу $x=1$. Значит используется ГУ $u(t,1,y)=e^{y}$, а разность берётся правая: $\dfrac{u_{j+1,k}-u_{j,k}}{h_x}$.
• Направление $y$: $v_2=-1<0$ — поток входит через границу $y=1$. Используется ГУ $u(t,x,1)=e^{x}$, разность правая: $\dfrac{u_{j,k+1}-u_{j,k}}{h_y}$.

Условия на выходных границах $u(t,0,y)=1$ и $u(t,x,0)=1$ в самой разностной схеме не используются (иначе задача была бы переопределена); они нужны лишь для контроля/анализа. Начальное условие $u(0,x,y)=1$ задаёт слой $n=0$.

Введём равномерную сетку $x_j=jh_x$ ($j=0,\dots,N_x$), $y_k=kh_y$ ($k=0,\dots,N_y$), $t_n=n\Delta t$, и обозначим $u^n_{j,k}\approx u(t_n,x_j,y_k)$. Для краткости далее $h_x=h_y=h$.

3. Неявная разностная схема методом дробных шагов

В методе дробных шагов (расщепление по направлениям) переход $n\to n+1$ разбивается на два полушага, и на каждом полушаге неявно аппроксимируется только один оператор переноса; источник $f$ относим ко второму полушагу.

Подсхема ① ($n\to n+1/2$): неявный перенос по $x$. Используем правую разность по $x$ (вход $x=1$):

$$\dfrac{u^{n+1/2}_{j,k}-u^{n}_{j,k}}{\Delta t}=\dfrac{u^{n+1/2}_{j+1,k}-u^{n+1/2}_{j,k}}{h}.$$

Подсхема ② ($n+1/2\to n+1$): неявный перенос по $y$ с источником. Правая разность по $y$ (вход $y=1$); источник берём на новом слое $f^{n+1}_{j,k}=-e^{t_{n+1}x_jy_k}$:

$$\dfrac{u^{n+1}_{j,k}-u^{n+1/2}_{j,k}}{\Delta t}=\dfrac{u^{n+1}_{j,k+1}-u^{n+1}_{j,k}}{h}-e^{\,t_{n+1}x_jy_k}.$$

Сумма двух подсхем (при $f$ на втором полушаге) аппроксимирует исходное $\dfrac{\partial u}{\partial t}-\dfrac{\partial u}{\partial x}-\dfrac{\partial u}{\partial y}=-e^{txy}$ с первым порядком по $\Delta t$ и по $h$.

4. Рекуррентные соотношения для каждой подсхемы

Каждая подсхема — двухдиагональная (в уравнение входят лишь два соседних узла), поэтому решается прямой подстановкой (бегущим счётом), прогонка не нужна. Диагональный коэффициент $(1+r)$ по модулю строго больше внедиагонального $|-r|=r$ при любых $\Delta t,h>0$ — есть диагональное преобладание, бегущий счёт устойчив.

Подсхема ① по $x$. Введём $r_x=\dfrac{\Delta t}{h}$. Умножим на $\Delta t$ и соберём неизвестные слоя $n+1/2$:

$$u^{n+1/2}_{j,k}-u^{n}_{j,k}=r_x\big(u^{n+1/2}_{j+1,k}-u^{n+1/2}_{j,k}\big)\;\Rightarrow\;(1+r_x)\,u^{n+1/2}_{j,k}-r_x\,u^{n+1/2}_{j+1,k}=u^{n}_{j,k}.$$

Поскольку известно граничное значение при $j=N_x$ ($u^{n+1/2}_{N_x,k}=e^{y_k}$), счёт ведём от границы $x=1$ в сторону убывания $j$ ($j=N_x-1,\dots,1$):

$$\boxed{\,u^{n+1/2}_{j,k}=\dfrac{u^{n}_{j,k}+r_x\,u^{n+1/2}_{j+1,k}}{1+r_x}\,},\qquad r_x=\dfrac{\Delta t}{h}.$$

Подсхема ② по $y$. Введём $r_y=\dfrac{\Delta t}{h}$. Умножим на $\Delta t$:

$$u^{n+1}_{j,k}-u^{n+1/2}_{j,k}=r_y\big(u^{n+1}_{j,k+1}-u^{n+1}_{j,k}\big)-\Delta t\,e^{\,t_{n+1}x_jy_k}.$$

Соберём неизвестные слоя $n+1$:

$$(1+r_y)\,u^{n+1}_{j,k}-r_y\,u^{n+1}_{j,k+1}=u^{n+1/2}_{j,k}-\Delta t\,e^{\,t_{n+1}x_jy_k}.$$

Известно граничное значение при $k=N_y$ ($u^{n+1}_{j,N_y}=e^{x_j}$), поэтому счёт ведём от границы $y=1$ в сторону убывания $k$ ($k=N_y-1,\dots,1$):

$$\boxed{\,u^{n+1}_{j,k}=\dfrac{u^{n+1/2}_{j,k}-\Delta t\,e^{\,t_{n+1}x_jy_k}+r_y\,u^{n+1}_{j,k+1}}{1+r_y}\,},\qquad r_y=\dfrac{\Delta t}{h}.$$

5. Устойчивость и порядок аппроксимации

Обе подсхемы неявные: знаменатели $1+r_x>0$ и $1+r_y>0$ положительны при любых $\Delta t,h>0$, а множитель перехода по модулю не превосходит единицы. Поэтому схема метода дробных шагов безусловно (абсолютно) устойчива — ограничения на шаг по времени нет (в отличие от явной схемы, где потребовалось бы условие типа Куранта $\dfrac{|v_1|\Delta t}{h_x}+\dfrac{|v_2|\Delta t}{h_y}\le 1$).

Порядок аппроксимации схемы — $O(\Delta t)+O(h)$ (первый порядок по времени за счёт расщепления и неявной аппроксимации и первый по пространству за счёт односторонних разностей).

Уравнение и граничные условия

Дано стационарное обыкновенное дифференциальное уравнение 2-го порядка

$$-\frac{du}{dx}+2x\,\frac{d^2u}{dx^2}=5u,\qquad u=u(x),\;x\in[0,1],$$ $$\frac{du}{dx}(x{=}0)=u(x{=}0),\qquad \frac{du}{dx}(x{=}1)=2\,u(x{=}1).$$

Оба граничных условия — 3-го рода (связывают значение функции и её производной на границе).

1. Приведение к виду для метода установления

Канонический вид ОДУ 2-го порядка курса (гл. 10.1, формула 10.11):

$$v\,\frac{du}{dx}=\sigma\,\frac{d^{2}u}{dx^{2}}-k\,u+f(x),\qquad \sigma>0.$$

Правила приведения: вторую производную ставим в правую часть с положительным коэффициентом $\sigma$, первую производную — в левую часть. Перепишем уравнение: $-u'+2x\,u''-5u=0$, домножим на $(-1)$ и перенесём $u'$ влево:

$$u'=2x\,u''-5u\;\;\Longrightarrow\;\; \underbrace{1}_{v}\cdot u'=\underbrace{2x}_{\sigma}\,u''-\underbrace{5}_{k}\,u+\underbrace{0}_{f}.$$ $$\boxed{v=1,\qquad \sigma=\sigma(x)=2x\ge0,\qquad k=5>0,\qquad f(x)=0.}$$

Коэффициент при второй производной переменный: $\sigma_j=2x_j$. На отрезке $[0,1]$ он неотрицателен, причём $\sigma=0$ при $x=0$. Так как при $u''$ обращается в нуль коэффициент диффузии, прямую прогонку для стационарной задачи применять неудобно; по условию используем метод установления.

Вводим фиктивную производную по «времени» (в левую часть, со знаком «плюс»), превращая стационарную задачу в нестационарную; искомая функция становится функцией двух переменных $u(x)\to\widetilde u(x,t)$:

$$\frac{\partial\widetilde u}{\partial t}+v\,\frac{\partial\widetilde u}{\partial x}=\sigma\,\frac{\partial^{2}\widetilde u}{\partial x^{2}}-k\,\widetilde u+f,$$ $$\frac{\partial\widetilde u}{\partial t}+\frac{\partial\widetilde u}{\partial x}=2x\,\frac{\partial^{2}\widetilde u}{\partial x^{2}}-5\,\widetilde u.$$

Это уравнение параболического типа. Граничные условия берутся из исходной задачи и от времени не зависят, поэтому при $t\to\infty$ решение «устанавливается»: $\widetilde u(x,t)\to u(x)$, $\partial\widetilde u/\partial t\to0$, и предел является решением исходной стационарной задачи. Так как $v=1>0$, для аппроксимации $\partial u/\partial x$ берём левую конечную разность $\dfrac{u_j-u_{j-1}}{h}$.

2. Неявная разностная схема

Вводим сетку по координате $x_j=(j-1)h$, $j=1,\dots,N$, $h=\dfrac{1}{N-1}$, и по «времени» (итерациям) с шагом $\Delta t$. Верхним индексом $n$ нумеруем итерации (временные слои). Записываем неявную схему (все пространственные операторы — на новом слое $n{+}1$), вторую производную аппроксимируем центральной разностью, конвективный член — левой разностью ($v>0$):

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+v\,\frac{u_{j}^{n+1}-u_{j-1}^{n+1}}{h}=\sigma_j\,\frac{u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}-k\,u_{j}^{n+1}+f,$$

с $v=1$, $\sigma_j=2x_j$, $k=5$, $f=0$. Неявная схема абсолютно устойчива, поэтому шаг $\Delta t$ можно брать грубым (произвольным) — это ускоряет установление. Порядок аппроксимации $O(\Delta t,h)$ (по «времени» — неявный Эйлер $O(\Delta t)$, $u''$ — центральная $O(h^2)$, конвективный член — левая «противопотоковая» разность $O(h)$, поэтому суммарно $O(\Delta t+h)$).

3. Приведение к виду для прогонки

Домножим на $\Delta t$ и сгруппируем неизвестные слоя $n{+}1$ по узлам $u_{j-1}^{n+1},\,u_{j}^{n+1},\,u_{j+1}^{n+1}$. Получаем трёхдиагональную систему

$$a_j\,u_{j+1}^{n+1}+b_j\,u_{j}^{n+1}+c_j\,u_{j-1}^{n+1}=\xi_j^{n},$$

с коэффициентами (общий вид):

$$a_j=-\sigma_j\frac{\Delta t}{h^{2}},\quad b_j=1+v\frac{\Delta t}{h}+2\sigma_j\frac{\Delta t}{h^{2}}+k\,\Delta t,\quad c_j=-v\frac{\Delta t}{h}-\sigma_j\frac{\Delta t}{h^{2}},\quad \xi_j^{n}=u_j^{n}+\Delta t\,f.$$

Подставляя $v=1$, $\sigma_j=2x_j$, $k=5$, $f=0$:

$$\boxed{a_j=-\,2x_j\frac{\Delta t}{h^{2}},\qquad b_j=1+\frac{\Delta t}{h}+4x_j\frac{\Delta t}{h^{2}}+5\,\Delta t,\qquad c_j=-\frac{\Delta t}{h}-2x_j\frac{\Delta t}{h^{2}},\qquad \xi_j^{n}=u_j^{n}.}$$

4. Проверка сходимости прогонки

Достаточное условие устойчивости прогонки — диагональное преобладание $|a_j|+|c_j|<|b_j|$. Все коэффициенты $a_j,c_j$ отрицательны, $b_j>0$, поэтому

$$|a_j|+|c_j|=2x_j\frac{\Delta t}{h^{2}}+\Big(\frac{\Delta t}{h}+2x_j\frac{\Delta t}{h^{2}}\Big)=\frac{\Delta t}{h}+4x_j\frac{\Delta t}{h^{2}},$$ $$|b_j|=1+\frac{\Delta t}{h}+4x_j\frac{\Delta t}{h^{2}}+5\,\Delta t.$$

Тогда

$$|b_j|-\big(|a_j|+|c_j|\big)=1+5\,\Delta t>0\quad\Rightarrow\quad |a_j|+|c_j|<|b_j|.$$

Неравенство строгое: запас $1+5\Delta t$ даёт фиктивная производная по времени (вклад $+1$) вместе с членом реакции $+k\,\Delta t=5\Delta t$ ($k>0$). Значит, на каждой итерации прогонка устойчива и сходится при любых $\Delta t,h>0$.

5. Итерационное (прогоночное) соотношение

На каждой итерации трёхдиагональная система решается методом прогонки. Итерационным выражением служит прогоночное соотношение (обратный ход, гл. 4.2.2):

$$u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j,\qquad j=N-1,\dots,1,$$

с прогоночными коэффициентами (прямой ход, $j=2,\dots,N-1$):

$$\alpha_j=\frac{-a_j}{b_j+c_j\,\alpha_{j-1}},\qquad \beta_j=\frac{\xi_j^{n}-c_j\,\beta_{j-1}}{b_j+c_j\,\alpha_{j-1}}.$$

6. Коэффициенты $\alpha_1,\ \beta_1$ (из левого ГУ)

Левое ГУ 3-го рода $\dfrac{du}{dx}(0)=u(0)$ аппроксимируем правой разностью в узле $j=1$:

$$\frac{u_2^{n+1}-u_1^{n+1}}{h}=u_1^{n+1}\;\Rightarrow\;u_2^{n+1}=(1+h)\,u_1^{n+1}\;\Rightarrow\;u_1^{n+1}=\frac{1}{1+h}\,u_2^{n+1}+0.$$

Сравнивая с $u_1^{n+1}=\alpha_1\,u_2^{n+1}+\beta_1$, получаем

$$\boxed{\alpha_1=\frac{1}{1+h},\qquad \beta_1=0.}$$

7. Решение на правой границе $u_N^{n+1}$ (из правого ГУ)

Правое ГУ 3-го рода $\dfrac{du}{dx}(1)=2\,u(1)$ аппроксимируем левой разностью в узле $j=N$ и подставляем $u_{N-1}^{n+1}=\alpha_{N-1}u_N^{n+1}+\beta_{N-1}$:

$$\frac{u_N^{n+1}-u_{N-1}^{n+1}}{h}=2\,u_N^{n+1}\;\Rightarrow\;u_N^{n+1}-\alpha_{N-1}u_N^{n+1}-\beta_{N-1}=2h\,u_N^{n+1},$$ $$\boxed{u_N^{n+1}=\frac{\beta_{N-1}}{\,1-2h-\alpha_{N-1}\,}.}$$

8. Начальное приближение и условие окончания итераций

Начальное (нулевое) приближение — нужно из-за фиктивной производной по времени; так как $f\equiv0$, удобно взять

$$u_j^{0}=f(x_j)=0,\qquad j=1,\dots,N$$

(в силу абсолютной устойчивости неявной схемы установление происходит при любом гладком старте). Итерации продолжают до установления — пока норма разности двух соседних приближений не станет меньше заданной точности $\varepsilon$:

$$\big\|u^{n+1}-u^{n}\big\|=\sqrt{h\sum_{j=1}^{N}\big(u_j^{n+1}-u_j^{n}\big)^{2}}\le\varepsilon.$$

Примечание. В данном варианте правая часть и ГУ однородны ($f\equiv0$, ГУ связывают $u$ и $u'$ без свободного члена), поэтому установившееся решение — тождественный нуль $u(x)\equiv0$; нулевой старт совпадает с ответом сразу. На корректность самой методики (приведение, знаки коэффициентов, диагональное преобладание, формулы прогонки и ГУ) это не влияет — все соотношения выше остаются верными для общего случая ненулевых $f$.

Алгоритм (метод установления, неявная схема)

1. Задать сетку $x_j=(j-1)h$, шаг итерации $\Delta t$, точность $\varepsilon$.
2. Задать нулевое приближение $u_j^{0}=0$.
3. Итерация $n\to n+1$ (прогонка):
   • из левого ГУ: $\alpha_1=\dfrac{1}{1+h}$, $\beta_1=0$;
   • прямой ход $j=2,\dots,N-1$: вычислить $a_j,b_j,c_j,\xi_j^{n}$ и $\alpha_j,\beta_j$;
   • из правого ГУ: $u_N^{n+1}=\dfrac{\beta_{N-1}}{1-2h-\alpha_{N-1}}$;
   • обратный ход $j=N-1,\dots,1$: $u_j^{n+1}=\alpha_j u_{j+1}^{n+1}+\beta_j$.
4. Проверить $\big\|u^{n+1}-u^{n}\big\|\le\varepsilon$. Если не выполнено — положить $u^{n}:=u^{n+1}$ и вернуться к п. 3; иначе решение установилось: $u_j\approx u_j^{n+1}$.

Постановка

Нестационарное трёхмерное по пространству уравнение $u=u(t,x,y,z)$:

$$\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}-2\frac{\partial u}{\partial y}-8\frac{\partial u}{\partial z}=7t\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial z^2}\right)+t^2,$$

с условиями

$$\begin{cases}u(t,0,y,z)=0,\\ u(t,1,y,z)=t,\end{cases}\qquad u(t,x,1,z)=zx,\qquad \begin{cases}u(t,x,y,0)=0,\\ u(t,x,y,1)=t^2,\end{cases}\qquad u(0,x,y,z)=y.$$

Вторые производные есть только по $x$ и $z$ — это параболические (диффузионные) направления. По $y$ присутствует только первая производная — чисто конвективное направление.

1. Канонический вид и коэффициенты

Сравниваем с общим видом уравнения с первыми производными (§9.10, (9.22)):

$$\frac{\partial u}{\partial t}+v_1\frac{\partial u}{\partial x}+v_2\frac{\partial u}{\partial y}+v_3\frac{\partial u}{\partial z}=\sigma\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial z^2}\right)-ku+f.$$ $$\sigma=7t,\qquad v_1=+1,\quad v_2=-2,\quad v_3=-8,\qquad k=0,\qquad f=t^2.$$

Правило выбора разностей для конвекции (против потока, §9.10.1). Одностороннюю разность выбирают так, чтобы усилить диагональ (обеспечить безусловное диагональное преобладание, M-матрицу). При $v>0$ — левая разность $\dfrac{u_{j}-u_{j-1}}{h}$ и требуется левое ГУ; при $v<0$ — правая разность $\dfrac{u_{j+1}-u_{j}}{h}$ и требуется правое ГУ:

• $v_1=+1>0$ (по $x$) — левая разность $\dfrac{u_{j}-u_{j-1}}{h_x}$;
• $v_2=-2<0$ (по $y$) — правая разность $\dfrac{u_{l+1}-u_{l}}{h_y}$;
• $v_3=-8<0$ (по $z$) — правая разность $\dfrac{u_{m+1}-u_{m}}{h_z}$.

Сеточные операторы ($j,l,m$ — индексы по $x,y,z$):

$$\Lambda_{xx}u=\frac{u_{j+1}-2u_{j}+u_{j-1}}{h_x^2},\qquad \Lambda_{zz}u=\frac{u_{m+1}-2u_{m}+u_{m-1}}{h_z^2},$$ $$\Lambda_{x}u=\frac{u_{j}-u_{j-1}}{h_x},\quad \Lambda_{y}u=-2\,\frac{u_{l+1}-u_{l}}{h_y},\quad \Lambda_{z}u=-8\,\frac{u_{m+1}-u_{m}}{h_z}.$$

2. Идея схемы предиктор-корректор (§9.8)

Шаг $\Delta t$ делится точкой $t^{n+1/2}=t^n+\Delta t/2$ на две половины. Предиктор расщепляет первую половину $\Delta t/2$ методом дробных шагов — по одному дробному шагу на каждое пространственное направление; в каждом дробном шаге «своё» направление берётся неявно (что и обеспечивает абсолютную устойчивость), остальные направления на этом шаге не участвуют. Знаменатель каждой подсхемы предиктора равен $\Delta t/2$. Результат предиктора — промежуточные значения $u^{n+1/2}$.

Чисто конвективное направление $y$ (§9.10.2). По $y$ нет второй производной, поэтому соответствующая подсхема — аналог одномерного уравнения первого порядка. Согласно §9.10.2 (схема (9.26)) её берут неявно и решают не прогонкой, а рекуррентным соотношением; такая подсхема «также абсолютно устойчива». Брать конвекцию по $y$ явно нельзя — явный upwind лишь условно устойчив ($|v_2|\Delta t/h_y\le1$). Поэтому предиктор состоит из трёх подсхем: ① неявно по $x$ (прогонка), ② неявно по $y$ (рекуррентно), ③ неявно по $z$ (прогонка). Каждый дифференциальный оператор учитывается ровно один раз — в своей подсхеме; свободный член $f$ — также один раз (в первой подсхеме, как $\xi^{n}=u^n+\tfrac{\Delta t}{2}f^n$, ср. §9.10.1).

Корректор (④) — единый полношаговый ($n\to n+1$) пересчёт с правой частью, аппроксимированной относительно $t^{n+1/2}$; он поднимает порядок по времени до второго (§9.8.2, §9.9, (9.21)). Обозначим $\tau=\Delta t/2$.

3. Предиктор

Подсхема ① ($n\to n+1/6$): неявно по $x$

Учитываются только операторы по $x$ (диффузия и конвекция) — неявно; свободный член $f^n$ учитывается здесь (один раз):

$$\frac{u_{j,l,m}^{\,n+1/6}-u_{j,l,m}^{\,n}}{\Delta t/2}+\Lambda_{x}u^{\,n+1/6}-7t\,\Lambda_{xx}u^{\,n+1/6}=f^{\,n},\qquad f^{\,n}=(t^n)^2.$$

На новом слое только $u_{j-1},u_{j},u_{j+1}$ — подсхема трёхдиагональна по $x$, решается прогонкой $u_{j}^{\,n+1/6}=\alpha_{j+1}u_{j+1}^{\,n+1/6}+\beta_{j+1}$. Так как $v_1>0$ (левая разность), конвективный вклад идёт в $a_j$ (множитель при $u_{j-1}$) и усиливает диагональ $b_j$. Умножая на $\tau=\Delta t/2$ и приводя к виду $a_j u_{j-1}+b_j u_{j}+c_j u_{j+1}=\xi_j$ ($t=t^{n+1/6}$):

$$a_j=-\frac{\Delta t}{2}\!\left(\frac{1}{h_x}+\frac{7t}{h_x^2}\right),\qquad b_j=1+\frac{\Delta t}{2}\!\left(\frac{1}{h_x}+\frac{2\cdot7t}{h_x^2}\right),\qquad c_j=-\frac{\Delta t}{2}\,\frac{7t}{h_x^2},$$ $$\xi_j=u_{j,l,m}^{\,n}+\frac{\Delta t}{2}\,f^{\,n},\qquad f^{\,n}=(t^n)^2.$$ $$\alpha_{j+1}=\frac{-c_j}{b_j+a_j\alpha_j},\qquad \beta_{j+1}=\frac{\xi_j-a_j\beta_j}{b_j+a_j\alpha_j},\qquad u_{j}^{\,n+1/6}=\alpha_{j+1}u_{j+1}^{\,n+1/6}+\beta_{j+1}.$$

Подсхема ② ($n+1/6\to n+1/3$): неявно по $y$ (рекуррентно)

Учитывается только конвекция по $y$, $v_2=-2<0$ (правая разность), неявно:

$$\frac{u_{j,l,m}^{\,n+1/3}-u_{j,l,m}^{\,n+1/6}}{\Delta t/2}+\Lambda_{y}u^{\,n+1/3}=0,\qquad \Lambda_{y}u^{\,n+1/3}=-2\,\frac{u_{l+1}^{\,n+1/3}-u_{l}^{\,n+1/3}}{h_y}.$$

Это аналог одномерного уравнения первого порядка — решается рекуррентным соотношением (§9.10.2). Разрешая относительно $u_l$ (через правый сосед, т.к. $v_2<0$; $\Delta t/2$ и множитель $2$ сокращаются):

$$\boxed{\;u_{l}^{\,n+1/3}=\frac{u_{l}^{\,n+1/6}+\dfrac{\Delta t}{h_y}\,u_{l+1}^{\,n+1/3}}{1+\dfrac{\Delta t}{h_y}} =\frac{h_y\,u_{l}^{\,n+1/6}+\Delta t\,u_{l+1}^{\,n+1/3}}{h_y+\Delta t}.\;}$$

Рекурсия идёт справа налево — от заданного правого ГУ на грани $y=1$ к меньшим $l$. Подсхема безусловно (абсолютно) устойчива благодаря неявности.

Подсхема ③ ($n+1/3\to n+1/2$): неявно по $z$

Учитываются только операторы по $z$ (диффузия и конвекция) — неявно:

$$\frac{u_{j,l,m}^{\,n+1/2}-u_{j,l,m}^{\,n+1/3}}{\Delta t/2}+\Lambda_{z}u^{\,n+1/2}-7t\,\Lambda_{zz}u^{\,n+1/2}=0.$$

Трёхдиагональна по $z$, прогонка $u_{m}^{\,n+1/2}=\tilde\alpha_{m+1}u_{m+1}^{\,n+1/2}+\tilde\beta_{m+1}$. Так как $v_3=-8<0$ (правая разность), конвективный вклад идёт в $\tilde c_m$ (при $u_{m+1}$) и усиливает диагональ $\tilde b_m$ ($t=t^{n+1/2}$):

$$\tilde a_m=-\frac{\Delta t}{2}\,\frac{7t}{h_z^2},\qquad \tilde b_m=1+\frac{\Delta t}{2}\!\left(\frac{8}{h_z}+\frac{2\cdot7t}{h_z^2}\right),\qquad \tilde c_m=-\frac{\Delta t}{2}\!\left(\frac{8}{h_z}+\frac{7t}{h_z^2}\right),$$ $$\tilde\xi_m=u_{j,l,m}^{\,n+1/3}.$$ $$\tilde\alpha_{m+1}=\frac{-\tilde c_m}{\tilde b_m+\tilde a_m\tilde\alpha_m},\qquad \tilde\beta_{m+1}=\frac{\tilde\xi_m-\tilde a_m\tilde\beta_m}{\tilde b_m+\tilde a_m\tilde\alpha_m},\qquad u_{m}^{\,n+1/2}=\tilde\alpha_{m+1}u_{m+1}^{\,n+1/2}+\tilde\beta_{m+1}.$$

Результат подсхем ①–③ — значения $u^{\,n+1/2}$ (предиктор завершён).

4. Корректор ④ ($n\to n+1$): повышение порядка по времени

Единый полношаговый пересчёт по всему $\Delta t$ с правой частью на полуслое $u^{\,n+1/2}$ (§9.8.2, §9.9, (9.20)–(9.21)). Здесь учитываются все операторы — каждый ровно один раз, и источник один раз:

$$\frac{u_{j,l,m}^{\,n+1}-u_{j,l,m}^{\,n}}{\Delta t} =7t\,\Lambda_{xx}u^{\,n+1/2}+7t\,\Lambda_{zz}u^{\,n+1/2}-\Lambda_{x}u^{\,n+1/2}-\Lambda_{y}u^{\,n+1/2}-\Lambda_{z}u^{\,n+1/2}+f^{\,n+1/2},$$

где конвективные операторы перенесены в правую часть из левой (отсюда знак «минус» перед $\Lambda_x,\Lambda_y,\Lambda_z$), $t=t^{n+1/2}$, $f^{\,n+1/2}=(t^{n+1/2})^2$. Рекуррентное соотношение корректора (9.21):

$$\boxed{\;u_{j,l,m}^{\,n+1}=u_{j,l,m}^{\,n}+\Delta t\Big(7t\,\Lambda_{xx}u^{\,n+1/2}+7t\,\Lambda_{zz}u^{\,n+1/2}-\Lambda_{x}u^{\,n+1/2}-\Lambda_{y}u^{\,n+1/2}-\Lambda_{z}u^{\,n+1/2}+f^{\,n+1/2}\Big).\;}$$

Все операторы — на известном полуслое $u^{\,n+1/2}$, поэтому корректор полностью явный (одно рекуррентное соотношение); абсолютная устойчивость всей схемы обеспечивается неявным предиктором (§9.8.2).

5. Начальное и граничные условия, старт

Начальное условие: $u^{0}_{j,l,m}=y_l$ (т.е. $u(0,x,y,z)=y$).

• Подсхема ① (по $x$): левое ГУ $u(t,0,y,z)=0\Rightarrow \alpha_1=0,\ \beta_1=0$; на правом конце $u_{J}=t$ (на нужном слое).
• Подсхема ② (по $y$): $v_2=-2<0$ — требуется только одно правое ГУ на грани $y=1$: $u(t,x,1,z)=zx$; от него стартует рекурсия справа налево. Других ГУ по $y$ не требуется (§9.10.2).
• Подсхема ③ (по $z$): левое ГУ $u(t,x,y,0)=0\Rightarrow \tilde\alpha_1=0,\ \tilde\beta_1=0$; на правом конце $u_{M}=t^2$.

6. Сходимость прогонок (диагональное преобладание)

Для подсхемы ① (по $x$): внедиагональные коэффициенты отрицательны безусловно, и

$$|a_j|+|c_j|=\frac{\Delta t}{2}\!\left(\frac{1}{h_x}+\frac{2\cdot7t}{h_x^2}\right)<1+\frac{\Delta t}{2}\!\left(\frac{1}{h_x}+\frac{2\cdot7t}{h_x^2}\right)=b_j,\qquad b_j-(|a_j|+|c_j|)=1>0.$$

Для подсхемы ③ (по $z$): $\ \tilde b_m-(|\tilde a_m|+|\tilde c_m|)=1>0$. Диагональное преобладание выполняется безусловно (за счёт единицы от временной производной), конвекция лишь усиливает диагональ. Рекуррентная подсхема ② по $y$ устойчива безусловно (знаменатель $1+\Delta t/h_y>1$). Никаких ограничений на число Пекле/Куранта не требуется.

7. Порядок аппроксимации

Корректор аппроксимирован относительно $t^{n+1/2}$ — временной оператор второго порядка $O(\Delta t^2)$. Диффузия $\Lambda_{xx},\Lambda_{zz}$ — центральные вторые разности $O(h_x^2),O(h_z^2)$; конвекция по всем трём координатам аппроксимирована односторонними (upwind) разностями — первый порядок. По $x$ и $z$ суммарную погрешность по пространству определяет худший (первый) порядок. Итог:

$$\psi=O\big(\Delta t^2,\ h_x,\ h_y,\ h_z\big).$$

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

Ответ

• Предиктор (половина $\Delta t/2$, три дробных шага, $\tau=\Delta t/2$; каждый оператор и источник — ровно один раз): ① неявно по $x$ (прогонка) с $a_j=-\tfrac{\Delta t}{2}\big(\tfrac{1}{h_x}+\tfrac{7t}{h_x^2}\big)$, $b_j=1+\tfrac{\Delta t}{2}\big(\tfrac{1}{h_x}+\tfrac{2\cdot7t}{h_x^2}\big)$, $c_j=-\tfrac{\Delta t}{2}\tfrac{7t}{h_x^2}$, $\xi_j=u^{n}+\tfrac{\Delta t}{2}(t^n)^2$; ② неявно по $y$ — рекуррентно $u_l=\dfrac{h_y u_l^{old}+\Delta t\,u_{l+1}^{new}}{h_y+\Delta t}$ ($v_2=-2<0$, правое ГУ $u(t,x,1,z)=zx$); ③ неявно по $z$ (прогонка) с $\tilde a_m=-\tfrac{\Delta t}{2}\tfrac{7t}{h_z^2}$, $\tilde b_m=1+\tfrac{\Delta t}{2}\big(\tfrac{8}{h_z}+\tfrac{2\cdot7t}{h_z^2}\big)$, $\tilde c_m=-\tfrac{\Delta t}{2}\big(\tfrac{8}{h_z}+\tfrac{7t}{h_z^2}\big)$.
• Корректор ④ — полношаговое явное рекуррентное соотношение $u^{n+1}=u^{n}+\Delta t\big(7t\Lambda_{xx}u^{n+1/2}+7t\Lambda_{zz}u^{n+1/2}-\Lambda_x u^{n+1/2}-\Lambda_y u^{n+1/2}-\Lambda_z u^{n+1/2}+f^{n+1/2}\big)$.
• Прогонки и рекурсия устойчивы безусловно; диагональное преобладание $b-(|a|+|c|)=1>0$. Порядок аппроксимации $O(\Delta t^2,\ h_x,\ h_y,\ h_z)$; схема абсолютно устойчива.