Теоретические вопросы: подробный разбор

Зачем нужна классификация

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

Общий вид уравнения

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

$$a_{11}\frac{\partial^{2}u}{\partial x^{2}}+2a_{12}\frac{\partial^{2}u}{\partial x\partial y}+a_{22}\frac{\partial^{2}u}{\partial y^{2}}=f\!\left(x,y,u,u_{x}',u_{y}'\right).$$

Дискриминант и три типа

Тип уравнения определяется по знаку величины (дискриминанта):

$$D=a_{12}^{2}-a_{11}a_{22}.$$

В зависимости от знака $D$ дифференциальные уравнения в частных производных 2-го порядка относят к уравнениям:

• при $D>0$ — гиперболического типа;
• при $D=0$ — параболического типа;
• при $D<0$ — эллиптического типа.

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

Стационарные и нестационарные процессы

Эллиптические уравнения описывают стационарные процессы (нет производной $\dfrac{\partial u}{\partial t}$, нет изменения во времени). Параболические и гиперболические уравнения описывают нестационарные процессы (есть $\dfrac{\partial u}{\partial t}$, есть изменение во времени).

Правила для многомерных уравнений

Принадлежность многомерных дифференциальных уравнений в частных производных 2-го порядка к тому или иному типу определяют по следующим правилам:

• Правило 1. Если в уравнении присутствуют производные 2-го порядка по всем независимым переменным и знаки перед ними одинаковые — уравнение относят к эллиптическому типу, например: $$\sigma_{x}\frac{\partial^{2}u}{\partial x^{2}}+\sigma_{y}\frac{\partial^{2}u}{\partial y^{2}}+\sigma_{z}\frac{\partial^{2}u}{\partial z^{2}}=f\!\left(t,x,y,z,u,u_{x}',u_{y}',u_{z}'\right).$$
• Правило 2. Если в уравнении отсутствует производная 2-го порядка хотя бы по одной из независимых переменных — уравнение относят к параболическому типу, например: $$\frac{\partial u}{\partial t}=\sigma_{x}\frac{\partial^{2}u}{\partial x^{2}}+\sigma_{y}\frac{\partial^{2}u}{\partial y^{2}}+\sigma_{z}\frac{\partial^{2}u}{\partial z^{2}}+f\!\left(t,x,y,z,u,u_{x}',u_{y}',u_{z}'\right).$$

Примеры определения типа

Пример 1. Уравнение теплопроводности (параболическое):

$$\frac{\partial u}{\partial t}=a_{11}\frac{\partial^{2}u}{\partial x^{2}}+f(x,t),\qquad a_{11}=1,\ a_{12}=0,\ a_{22}=0;$$ $$D=0-1\cdot 0=0\ \Rightarrow\ \text{параболическое}.$$

Пример 2. Уравнение Лапласа (эллиптическое):

$$\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}=0,\qquad a_{11}=1,\ a_{12}=0,\ a_{22}=1;$$ $$D=0^{2}-1\cdot 1=-1<0\ \Rightarrow\ \text{эллиптическое}.$$

Пример 3. Уравнение Пуассона (эллиптическое):

$$\lambda\left(\frac{\partial^{2}T}{\partial x^{2}}+\frac{\partial^{2}T}{\partial y^{2}}\right)=q(x,y),\qquad a_{11}=\lambda,\ a_{12}=0,\ a_{22}=\lambda;$$ $$D=0^{2}-\lambda\cdot\lambda=-\lambda^{2}<0\ \Rightarrow\ \text{эллиптическое}.$$

Пример 4. Уравнение конвекции-диффузии в стационарном режиме (эллиптическое):

$$v\frac{\partial C}{\partial x}=D_{L}\frac{\partial^{2}C}{\partial x^{2}}+D_{R}\frac{\partial^{2}C}{\partial r^{2}}+\frac{D_{R}}{r}\frac{\partial C}{\partial r}-kC,$$ $$a_{11}=D_{L},\ a_{12}=0,\ a_{22}=D_{R}\ \Rightarrow\ D=-D_{L}\cdot D_{R}<0\ \Rightarrow\ \text{эллиптическое}.$$

Пример 5. Трубчатый реактор в нестационарном режиме (параболическое):

$$\frac{\partial C}{\partial t}+v\frac{\partial C}{\partial x}=D_{L}\frac{\partial^{2}C}{\partial x^{2}}-kC,\qquad a_{11}=D_{L},\ a_{12}=0,\ a_{22}=0;$$ $$D=0-D_{L}\cdot 0=0\ \Rightarrow\ \text{параболическое}.$$

Дальнейшие примеры разбора типа уравнений и обезразмеривания приведены в Семинаре 0 и Семинаре 1. Дифференциальные уравнения 3-го и более высоких порядков в курсе не рассматриваются.

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

Для решения дифференциальных уравнений численными методами требуются дополнительные условия. Если искомая функция (концентрация, температура и т.д.) является функцией времени $u=u(t)$, требуются начальные условия, характеризующие значение функции в момент времени, принятый за начальный:

$$u(t=0)=u^{0}.$$

Если искомая функция также является функцией пространственных координат $u=u(t,x)$, то начальные условия характеризуют её распределение в пространстве в начальный момент времени:

$$u(t=0,x)=u^{0}(x).$$

В этом случае помимо начальных условий требуются ещё и граничные условия, характеризующие значение функции $u(t,x)$ на границе изучаемой системы с внешней средой для любого момента времени. Если функция зависит от нескольких пространственных координат, граничные условия задают по каждой из них. Количество граничных условий по каждой координате определяется порядком старшей производной функции $u$ по этой координате в уравнении (подробнее — Начальные и граничные условия).

Уравнение, на примере которого рассматривают классификацию

Классификацию граничных условий рассматривают на примере уравнения теплового баланса трубчатого реактора с продольным перемешиванием:

$$\rho C_{T}\left[\frac{\partial T}{\partial t}+v\frac{\partial T}{\partial x}\right]=\lambda\frac{\partial^{2}T}{\partial x^{2}}+\Delta H\,w,$$

где $w$ — скорость реакции; $\Delta H$ — тепловой эффект реакции; $C_{T},\rho,T$ — теплоёмкость, плотность и температура смеси; $v$ — линейная скорость потока; $\lambda$ — коэффициент теплопроводности; $x$ — координата по длине реактора. Начальное условие характеризует распределение температуры по длине реактора в начальный момент времени: $T(t=0,x)=T^{0}(x)$.

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

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

$$T(t,x=0)=\varphi_{1}(t),\qquad T(t,x=l)=\varphi_{2}(t).$$

Граничные условия 2-го рода (Неймана)

Задают изменение функции (производную по координате) на границах реактора для любого момента времени:

$$\frac{\partial T}{\partial x}(t,x=0)=\varphi_{1}(t),\qquad \frac{\partial T}{\partial x}(t,x=l)=\varphi_{2}(t).$$

Граничные условия 3-го рода (Робена)

Определяют закон свободного теплообмена с окружающей средой на границах реактора для любого момента времени:

$$\frac{\partial T}{\partial x}(t,x=0)=\alpha\{T(t,x=0)-T_{cp}\},$$ $$\frac{\partial T}{\partial x}(t,x=l)=\alpha\{T(t,x=l)-T_{cp}\},$$

где $\alpha$ — коэффициент теплоотдачи; $l$ — длина реактора; $T_{cp}$ — температура окружающей среды.

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

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

$$T(t,x=0)=\varphi_{1}(t),\qquad \frac{\partial T}{\partial x}(t,x=l)=\alpha\{T(t,x=l)-T_{cp}\}.$$

Обобщённая запись ГУ 3-го рода (для аппроксимации)

Граничные условия 3-го рода можно записать в обобщённом виде:

$$\frac{\partial u}{\partial x}(t,x=a)=\varphi_{1}(t)\,u(t,x=a)+\psi_{1}(t),$$ $$\frac{\partial u}{\partial x}(t,x=b)=\varphi_{2}(t)\,u(t,x=b)+\psi_{2}(t).$$

Физический смысл выбора рода ГУ (пример)

Для уравнения поперечной диффузии $\dfrac{\partial C}{\partial t}=D_r\dfrac{\partial^{2}C}{\partial r^{2}}-kC$ (параболический тип) род ГУ выбирают из физики процесса:

• в центре трубы ($r=0$) градиент концентрации равен нулю — условие II рода: $\frac{\partial C}{\partial r}(t,r=0)=0$;
• у стенки ($r=R$), где поток вещества отсутствует, — условие I рода: $C(t,r=R)=0$;
• условие III рода применяется на границе раздела фаз, где идёт массо- или теплопередача: $D\,\frac{\partial C}{\partial x}(t,x=l)=\beta[C(t,x=l)-C_{\text{ср}}(t)]$, где $\beta$ — коэффициент массопередачи.

Подробные примеры — в Семинаре 0 и Семинаре 2. О том, как граничные условия I–III рода записывают в разностном виде, см. вопрос об аппроксимации НУ и ГУ.

Постановка

Обыкновенными называют дифференциальные уравнения, в которых искомая функция зависит от одной переменной (времени или координаты). Рассмотрим примеры математических моделей химических реакторов, в которых протекает простая необратимая реакция типа

$$n\mathrm{X}\to\mathrm{P},$$

скорость которой определяется формулой $w=kc^{n}$, где $k$ — константа скорости реакции, $c$ — концентрация вещества X. Типы уравнений и обозначения курса описаны в Типах дифференциальных уравнений.

1) Проточный реактор идеального смешения (ОДУ 1-го порядка)

Модель включает материальный и тепловой балансы по концентрации компонента X и температуре:

$$V\frac{dc}{dt}=v_{q}(c_{0}-c)-wV,\qquad V\rho C_{T}\frac{dT}{dt}=v_{q}(\rho_{0}C_{T0}T_{0}-\rho C_{T}T)+\Delta H\,wV,$$

где $v_{q}$ — объёмный расход поступающего раствора; $\Delta H$ — тепловой эффект реакции; $V$ — рабочий объём реактора; $C_{T},\rho,T$ — теплоёмкость, плотность и температура смеси; индекс (0) — значение на входе. Модель состоит из двух обыкновенных дифференциальных уравнений 1-го порядка, которые дополняют начальными условиями (задача Коши):

$$c(t=0)=c^{0},\qquad T(t=0)=T^{0}.$$

2) Трубчатый реактор идеального вытеснения, стационар (ОДУ 1-го порядка)

$$v\frac{dc}{dx}=-w,\qquad v\rho C_{T}\frac{dT}{dx}=\Delta H\,w,$$

где $v$ — линейная скорость потока; $x$ — координата по длине реактора. Модель состоит из двух обыкновенных дифференциальных уравнений 1-го порядка, которые дополняют граничными условиями:

$$c(x=0)=c_{0},\qquad T(x=0)=T_{0}.$$

3) Трубчатый реактор с продольным перемешиванием, стационар (ОДУ 2-го порядка)

$$v\frac{dc}{dx}=D_{L}\frac{d^{2}c}{dx^{2}}-w,\qquad v\rho C_{T}\frac{dT}{dx}=\lambda\frac{d^{2}T}{dx^{2}}+\Delta H\,w,$$

где $D_{L},\lambda$ — коэффициенты диффузии и теплопроводности. Модель состоит из двух обыкновенных дифференциальных уравнений 2-го порядка; так как старшая производная по координате — второго порядка, требуются по два граничных условия на каждое уравнение:

$$\begin{cases}c(x=0)=c_{0}\\ c(x=l)=c_{l}\end{cases}\qquad \begin{cases}T(x=0)=T_{0}\\ T(x=l)=T_{l}\end{cases}$$

здесь $l$ — длина реактора.

Вывод

Таким образом, ОДУ возникают в моделях реакторов, когда искомая функция (концентрация, температура) зависит лишь от одной переменной: от времени $t$ (реактор смешения, нестационар) или от координаты $x$ (трубчатые реакторы в стационарном режиме). Модели с ОДУ 1-го порядка по времени дополняют начальными условиями, по координате — граничными; модели с ОДУ 2-го порядка требуют двух граничных условий по координате. Пример ОДУ 2-го порядка с разбором (движение тела в атмосфере) приведён в Семинаре 1.

Постановка

Дифференциальными уравнениями в частных производных описывают процессы, в которых искомая функция зависит от двух и более переменных (например, $c=c(t,x)$). Рассмотрим математические модели химических реакторов с простой необратимой реакцией $nX\to P$ (скорость $w=kc^{n}$). Тип каждого уравнения определяют по дискриминанту $D=a_{12}^{2}-a_{11}a_{22}$ (см. Классификацию уравнений 2-го порядка).

4) Трубчатый реактор, нестационар (УрЧП 1-го порядка)

Материальный баланс по концентрации компонента X:

$$\frac{\partial c}{\partial t}+v\frac{\partial c}{\partial x}=-w,\qquad c=c(t,x).$$

Уравнение дополняют начальным и граничным условиями:

$$c(t=0,x)=c^{0}(x),\qquad c(t,x=0)=\varphi(t).$$

5) Трубчатый реактор с продольным перемешиванием, нестационар (параболическое)

$$\frac{\partial c}{\partial t}+v\frac{\partial c}{\partial x}=D_{L}\frac{\partial^{2}c}{\partial x^{2}}-w,\qquad c=c(t,x).$$

Определим тип уравнения:

$$a_{11}=D_{L},\ a_{22}=0,\ a_{12}=0\ \Rightarrow\ D=0-D_{L}\cdot 0=0,$$

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

$$c(t=0,x)=c^{0}(x),\qquad c(t,x=0)=\varphi_{1}(t),\qquad c(t,x=l)=\varphi_{2}(t).$$

6) Реактор с продольным и поперечным перемешиванием, нестационар (параболическое, многомерное)

$$\frac{\partial c}{\partial t}+v\frac{\partial c}{\partial x}=D_{L}\frac{\partial^{2}c}{\partial x^{2}}+D_{R}\frac{\partial^{2}c}{\partial r^{2}}+\frac{D_{R}}{r}\frac{\partial c}{\partial r}-w,$$

где $r$ — координата по радиусу; $D_{L},D_{R}$ — коэффициенты диффузии в продольном и поперечном направлениях. Уравнение многомерное, $c=c(t,x,r)$. Так как отсутствует производная 2-го порядка по времени, уравнение параболического типа. Его дополняют начальным и граничными условиями:

$$c(t=0,x,r)=c^{0}(x,r),$$ $$\begin{cases}c(t,x=0,r)=\varphi_{1}(t,r)\\ c(t,x=l,r)=\varphi_{2}(t,r)\end{cases}\qquad \begin{cases}c(t,x,r=0)=\psi_{1}(t,x)\\ c(t,x,r=R)=\psi_{2}(t,x)\end{cases}$$

здесь $R$ — радиус реактора.

7) Реактор с продольным и поперечным перемешиванием, стационар (эллиптическое)

$$v\frac{\partial c}{\partial x}=D_{L}\frac{\partial^{2}c}{\partial x^{2}}+D_{R}\frac{\partial^{2}c}{\partial r^{2}}+\frac{D_{R}}{r}\frac{\partial c}{\partial r}-w,\qquad c=c(x,r).$$

Определим тип уравнения:

$$a_{11}=D_{L},\ a_{22}=D_{R},\ a_{12}=0\ \Rightarrow\ D=0-D_{L}\cdot D_{R}<0,$$

следовательно, уравнение эллиптического типа. Его дополняют граничными условиями:

$$\begin{cases}c(t,x=0,r)=\varphi_{1}(t,r)\\ c(t,x=l,r)=\varphi_{2}(t,r)\end{cases}\qquad \begin{cases}c(t,x,r=0)=\psi_{1}(t,x)\\ c(t,x,r=R)=\psi_{2}(t,x)\end{cases}$$

Вывод

УрЧП возникают, когда искомая функция зависит от нескольких переменных. Нестационарные модели реакторов дают уравнения параболического (есть $\partial/\partial t$, нет второй производной по времени) типа и УрЧП 1-го порядка, а стационарные модели с диффузией по двум координатам — уравнения эллиптического типа. Дополнительные примеры обезразмеривания и определения типа таких уравнений — в Семинаре 1 и Семинаре 0.

Что такое аппроксимация

В основе изучаемых методов численного решения лежит преобразование дифференциальной задачи в разностную, называемое аппроксимацией. Прежде чем аппроксимировать целое уравнение, рассматривают аппроксимацию простейших дифференциальных операторов — производных первого и второго порядков. Интервал изменения переменной $x\in[a,b]$ разбивают на $n$ равных частей; вводят обозначения: $j$ — номер точки деления, $u(x_{j})=u_{j}$ — значение функции в точке $x_{j}$, $x_{j+1}-x_{j}=\Delta x=h$ — шаг (см. Разностную аппроксимацию производной 1-го порядка).

Аппроксимация $\partial u/\partial x$ (производная 1-го порядка)

Производную $\dfrac{du}{dx}\big|_{x_{j}}$ аппроксимируют тремя разностными операторами:

• правая конечная разность: $\displaystyle \lambda_{x}^{+}u=\frac{u_{j+1}-u_{j}}{h}$;
• левая конечная разность: $\displaystyle \lambda_{x}^{-}u=\frac{u_{j}-u_{j-1}}{h}$;
• центральная конечная разность: $\displaystyle \lambda_{x}^{0}u=\frac{u_{j+1}-u_{j-1}}{2h}$.

Также аппроксимацию можно задать линейной комбинацией правой и левой разностей:

$$\lambda_{x}^{\sigma}u=\sigma\frac{u_{j+1}-u_{j}}{h}+(1-\sigma)\frac{u_{j}-u_{j-1}}{h},\qquad 0\le\sigma\le 1,$$

при $\sigma=0$ — левая разность, при $\sigma=1$ — правая, при $\sigma=1/2$ — центральная.

Аппроксимация $\partial u/\partial t$ (производная по времени)

Производную по времени в точке $(t^{n},x_{j})$ разностной сетки аппроксимируют конечной разностью, стабилизируя значение координаты $x$ в точке $j$. В стандартной (явной) схеме используют правую конечную разность по времени:

$$\left.\frac{\partial u}{\partial t}\right|_{t^{n},x_{j}}\ \rightarrow\ \frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t},$$

где $\Delta t$ — шаг по времени, $u_{j}^{n}=u(t^{n},x_{j})$. Это та же правая конечная разность, что и для $\partial u/\partial x$, но по переменной $t$.

Аппроксимация производной 2-го порядка (в т.ч. $\partial^{2}u/\partial t^{2}$)

Поскольку первая производная $\dfrac{du}{dx}=w(x)$ — функция той же переменной, вторую производную представляют как первую производную функции $w$. Последовательно применяя правую и левую разности, получают разностный оператор для второй производной:

$$\left.\frac{d^{2}u}{dx^{2}}\right|_{x_{j}}\ \rightarrow\ \lambda_{xx}u=\frac{\dfrac{u_{j+1}-u_{j}}{h}-\dfrac{u_{j}-u_{j-1}}{h}}{h}=\frac{u_{j+1}-2u_{j}+u_{j-1}}{h^{2}}.$$

Тот же приём («разность разностей») применяется к оператору второй производной по любой координате. Важное замечание по курсу: разностные схемы здесь строятся двухслойными по времени (используются только слои $n$ и $n+1$), поэтому вторая производная по времени $\dfrac{\partial^{2}u}{\partial t^{2}}$ как самостоятельный разностный оператор в материалах курса не вводится — оператор второй производной разбирается именно для пространственной координаты ($\lambda_{xx}$, см. раздел 2.2.3 учебника).

Понятие порядка аппроксимации

Чтобы выяснить точность аппроксимации, значения $u_{j+1},u_{j-1}$ раскладывают в ряд Тейлора относительно точки $x_{j}$:

$$u_{j+1}=u_{j}+u_{j}'h+u_{j}''\frac{h^{2}}{2!}+u_{j}'''\frac{h^{3}}{3!}+\dots,$$ $$u_{j-1}=u_{j}-u_{j}'h+u_{j}''\frac{h^{2}}{2!}-u_{j}'''\frac{h^{3}}{3!}+\dots$$

Подставляя в правую разность:

$$\lambda_{x}^{+}u=\frac{u_{j+1}-u_{j}}{h}=u_{j}'+u_{j}''\frac{h}{2!}+u_{j}'''\frac{h^{2}}{3!}+\dots$$

Первое слагаемое — истинная производная, остальные составляют ошибку аппроксимации. При $h\to 0$ наибольший вклад в ошибку вносит слагаемое с наименьшей степенью $h$. Порядок аппроксимации по $h$ равен этой наименьшей степени $h$ в ошибке: чем выше порядок, тем точнее аппроксимация и меньше её ошибка.

Отсюда порядки операторов:

• правая конечная разность: $\lambda_{x}^{+}u=u_{j}'+O(h)$ — первый порядок;
• левая конечная разность: $\lambda_{x}^{-}u=u_{j}'+O(h)$ — первый порядок;
• центральная конечная разность: $\lambda_{x}^{0}u=u_{j}'+O(h^{2})$ — второй порядок (точнее левой и правой);
• оператор второй производной: $\lambda_{xx}u=u_{j}''+u_{j}^{IV}\dfrac{h^{2}}{12}=u_{j}''+O(h^{2})$ — второй порядок.

Порядок аппроксимации разностной схемы

Разностную схему составляют из отдельных разностных операторов; каждый оператор имеет свой порядок аппроксимации, поэтому и схема имеет порядок аппроксимации — по каждой независимой переменной отдельно. Для явной схемы уравнения параболического типа $\dfrac{\partial u}{\partial t}=\sigma\dfrac{\partial^{2}u}{\partial x^{2}}+f$:

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

разложение в ряд Тейлора около точки $(t^{n},x_{j})$ даёт:

$$\left.\frac{\partial u}{\partial t}\right|^{n}_{j}+O(\Delta t)=\sigma\left.\frac{\partial^{2}u}{\partial x^{2}}\right|^{n}_{j}+O(h^{2})+f.$$

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

$$O(\Delta t)+O(h^{2})\qquad\text{или}\qquad O(\Delta t,h^{2}).$$

Неявная схема имеет тот же порядок аппроксимации. Примеры определения порядка аппроксимации конкретных схем (через ряд Тейлора) разобраны в Семинаре 2; подробности — в учебнике: Понятие порядка аппроксимации, Аппроксимация производной 2-го порядка, Порядок аппроксимации разностной схемы.

Постановка: что аппроксимируем

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

$$\frac{\partial u}{\partial t}=\sigma\frac{\partial^{2} u}{\partial x^{2}}+f(u,t,x);\quad u(t{=}0,x)=\xi(x);\quad\begin{cases}u(t,x{=}a)=\varphi_{1}(t)\\ u(t,x{=}b)=\varphi_{2}(t)\end{cases}$$

Производная по времени аппроксимируется правой разностью, производная второго порядка по координате — центральной разностью. Принципиальное различие явной и неявной схем — на каком временно́м слое (n или n+1) берётся пространственный оператор.

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

Пространственный оператор записывается на известном n-м шаге по времени (учебник 4.1.1):

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}=\sigma\frac{u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}+f\left(u_{j}^{n},t^{n},x_{j}\right)$$

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

$$u_{j}^{n+1}=u_{j}^{n}+\frac{\sigma\Delta t}{h^{2}}\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)+f\left(u_{j}^{n},t^{n},x_{j}\right)\Delta t$$

Такое равенство называют рекуррентным соотношением: оно позволяет рассчитывать значение функции в узле напрямую через известные значения в других (соседних) узлах. Значения $u_{j}^{0}$ берутся из начального условия, краевые узлы $u_{1}^{n+1},u_{N}^{n+1}$ — из аппроксимации граничных условий, остальные узлы $j=2,\dots,N-1$ — по рекуррентной формуле.

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

Пространственный оператор записывается на искомом (n+1)-м шаге по времени (учебник 4.2.1):

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}=\sigma\frac{u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}+f\left(t^{n},x_{j}\right)$$

Здесь в одно уравнение входят три неизвестных значения нового слоя — $u_{j-1}^{n+1},u_{j}^{n+1},u_{j+1}^{n+1}$. Поэтому выразить $u_{j}^{n+1}$ напрямую (как в явной схеме) нельзя: для всех внутренних узлов получается система линейных уравнений с трёхдиагональной матрицей. Её решают специальным методом — методом прогонки (см. вопрос 9).

Сравнение свойств

• Порядок аппроксимации у обеих схем одинаковый: $O\left(\Delta t,h^{2}\right)$ — первый порядок по времени и второй по координате.
• Устойчивость. Явная схема условно устойчива: её устойчивость зависит от ограничения на шаги сетки $\dfrac{\Delta t}{h^{2}}\le\dfrac{1}{2\sigma}$ (см. вопрос 7). Неявная схема абсолютно устойчива — устойчива при любых $\Delta t$ и $h$ (см. вопрос 8).
• Метод решения. Явная схема имеет очень простой метод решения (прямая подстановка в рекуррентную формулу) — это её достоинство. Неявная схема требует более сложного метода прогонки, но зато не накладывает ограничений на шаг.

Таким образом, выбор между схемами — это компромисс: простота расчёта (явная) против свободы выбора шага и устойчивости (неявная). Блок-схемы алгоритмов обеих схем приведены на странице Блок-схемы, а пошаговые примеры записи и решения — в Семинаре 4.

Идея спектрального (гармонического) метода

Для исследования устойчивости разностных схем в курсе используется гармонический (спектральный) метод анализа (учебник 3.2). Погрешность решения $z_{j}^{n}$ удовлетворяет той же разностной схеме, что и сама функция; погрешность представляют в виде комплексной гармоники:

$$z_{j}^{n}=\lambda^{n}e^{i\alpha j}$$

где $\lambda$ — собственное число оператора перехода, $i$ — мнимая единица, $\alpha$ — фаза гармоники. Необходимое условие устойчивости (норма погрешности не должна расти, $\|z^{n+1}\|\le\|z^{n}\|$) принимает вид условия на собственные числа:

$$|\lambda|\le 1$$

Шаг 1. Отбрасываем свободный член

Сравнивая явную схему параболического уравнения с разностной схемой для погрешности её решения, видим, что они совпадают по структуре. Значит, наличие свободного члена $f(t,x)$ не влияет на устойчивость (при условии, что он не содержит искомую функцию). Поэтому исследуем однородную явную схему (учебник 3.3):

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}=\sigma\frac{u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}$$

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

Подставляем $u_{j}^{n}=\lambda^{n}e^{i\alpha j}$:

$$\frac{\lambda^{n+1}e^{i\alpha j}-\lambda^{n}e^{i\alpha j}}{\Delta t}=\sigma\frac{\lambda^{n}e^{i\alpha(j+1)}-2\lambda^{n}e^{i\alpha j}+\lambda^{n}e^{i\alpha(j-1)}}{h^{2}}$$

Делим обе части на общий множитель $\lambda^{n}e^{i\alpha j}$:

$$\frac{\lambda-1}{\Delta t}=\sigma\frac{e^{i\alpha}-2+e^{-i\alpha}}{h^{2}}$$

Шаг 3. Тригонометрические преобразования

Используем формулу Эйлера $e^{\pm i\alpha}=\cos\alpha\pm i\sin\alpha$, откуда числитель $e^{i\alpha}-2+e^{-i\alpha}=2\cos\alpha-2$. С помощью тождества $\cos\alpha=1-2\sin^{2}\dfrac{\alpha}{2}$ получаем ключевое соотношение (используется в каждом примере, семинар 3):

$$e^{i\alpha}-2+e^{-i\alpha}=2\cos\alpha-2=-4\sin^{2}\frac{\alpha}{2}$$

Тогда уравнение принимает вид и позволяет выразить $\lambda$:

$$\frac{\lambda-1}{\Delta t}=\sigma\frac{-4\sin^{2}\dfrac{\alpha}{2}}{h^{2}}\quad\Rightarrow\quad\lambda=1-\frac{4\sigma\Delta t}{h^{2}}\sin^{2}\frac{\alpha}{2}$$

Шаг 4. Применяем условие |λ| ≤ 1

Подставляем выражение для $\lambda$ в необходимое условие устойчивости:

$$|\lambda|\le 1\quad\Rightarrow\quad-1\le 1-\frac{4\sigma\Delta t}{h^{2}}\sin^{2}\frac{\alpha}{2}\le 1$$

Правое неравенство $1-\dfrac{4\sigma\Delta t}{h^{2}}\sin^{2}\dfrac{\alpha}{2}\le 1$ выполняется автоматически при $\sigma>0$ (вычитаемое неотрицательно). Рассмотрим левое неравенство:

$$1-\frac{4\sigma\Delta t}{h^{2}}\sin^{2}\frac{\alpha}{2}\ge-1\quad\Rightarrow\quad\frac{\Delta t}{h^{2}}\sin^{2}\frac{\alpha}{2}\le\frac{1}{2\sigma}$$

Чтобы условие выполнялось при любой фазе $\alpha$, задаём для $\sin^{2}\dfrac{\alpha}{2}$ максимальное значение, равное 1, и переходим к более строгому условию:

$$\boxed{\frac{\Delta t}{h^{2}}\le\frac{1}{2\sigma}}$$

Вывод

Полученное неравенство — условие устойчивости явной разностной схемы, аппроксимирующей одномерное уравнение параболического типа. Поскольку устойчивость зависит от ограничения, накладываемого на выбор интервалов деления сетки (соотношение $\Delta t$ и $h$), такую схему называют условно устойчивой (учебник 3.4). Практическое следствие: при уменьшении шага по пространству шаг по времени приходится резко уменьшать (например, при $\sigma=10$, $h=10^{-1}$ получается $\Delta t\le 5\cdot10^{-4}$), что делает явную схему вычислительно дорогой — это разобрано на примерах в Семинаре 3.

Исследуемая схема и метод

Рассматривается неявная разностная схема параболического уравнения (учебник 3.5 = раздел 3.4 методички):

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}=\sigma\frac{u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}+f\left(t^{n},x_{j}\right)$$

Как и для явной схемы, разностная схема для погрешности по структуре совпадает с однородной неявной схемой, поэтому свободный член $f(t,x)$ отбрасываем (он не влияет на устойчивость, если не содержит искомую функцию) и представляем решение в виде гармоники $u_{j}^{n}=\lambda^{n}e^{i\alpha j}$ — тот же спектральный метод, что и в вопросе 7.

Шаг 1–2. Подстановка гармоники

Подставляем $u_{j}^{n}=\lambda^{n}e^{i\alpha j}$. Принципиальное отличие от явной схемы: пространственный оператор записан на слое (n+1), поэтому в его трёх членах появляется множитель $\lambda^{n+1}$:

$$\frac{\lambda^{n+1}e^{i\alpha j}-\lambda^{n}e^{i\alpha j}}{\Delta t}=\sigma\frac{\lambda^{n+1}e^{i\alpha(j+1)}-2\lambda^{n+1}e^{i\alpha j}+\lambda^{n+1}e^{i\alpha(j-1)}}{h^{2}}$$

Делим обе части на $\lambda^{n}e^{i\alpha j}$. В правой части остаётся дополнительный множитель $\lambda$:

$$\frac{\lambda-1}{\Delta t}=\sigma\lambda\frac{e^{i\alpha}-2+e^{-i\alpha}}{h^{2}}$$

Шаг 3. Тригонометрическое тождество

Применяем то же тождество, что и для явной схемы:

$$e^{i\alpha}-2+e^{-i\alpha}=2\cos\alpha-2=-4\sin^{2}\frac{\alpha}{2}$$

Подставляем:

$$\frac{\lambda-1}{\Delta t}=\sigma\lambda\frac{-4\sin^{2}\dfrac{\alpha}{2}}{h^{2}}$$

Шаг 4. Выражаем λ

В отличие от явной схемы, здесь $\lambda$ входит в обе части. Группируем члены с $\lambda$ в левой части:

$$\lambda-1=-\frac{4\sigma\Delta t}{h^{2}}\sin^{2}\frac{\alpha}{2}\,\lambda\quad\Rightarrow\quad\lambda\left(1+\frac{4\sigma\Delta t}{h^{2}}\sin^{2}\frac{\alpha}{2}\right)=1$$

Откуда:

$$\boxed{\lambda=\frac{1}{1+\dfrac{4\sigma\Delta t}{h^{2}}\sin^{2}\dfrac{\alpha}{2}}}$$

Шаг 5. Анализ |λ| и вывод

В знаменателе стоит единица плюс заведомо неотрицательная величина (при $\sigma>0$, $\Delta t>0$, $h>0$ и $\sin^{2}\dfrac{\alpha}{2}\ge 0$). Следовательно, знаменатель всегда $\ge 1$, а значит:

$$0<\lambda\le 1\quad\Rightarrow\quad|\lambda|\le 1$$

Необходимое условие устойчивости $|\lambda|\le 1$ выполняется при любых значениях $\Delta t$ и $h$ — никаких ограничений на шаги сетки не возникает. Такие схемы, устойчивость которых не зависит от выбора интервалов деления на разностной сетке, называют абсолютно устойчивыми.

Сопоставление с явной схемой

В явной схеме появлялась разность $\lambda=1-\dfrac{4\sigma\Delta t}{h^{2}}\sin^{2}\dfrac{\alpha}{2}$, которая могла стать $<-1$ при большом шаге — отсюда условие $\dfrac{\Delta t}{h^{2}}\le\dfrac{1}{2\sigma}$. В неявной схеме та же величина оказалась в знаменателе дроби, поэтому $\lambda$ заведомо не превосходит 1. Это — главное преимущество неявной схемы; платой за него служит более сложный метод решения (метод прогонки, вопрос 9). Абсолютная устойчивость неявной схемы на конкретных примерах проверена в Семинаре 3.

Зачем нужен метод прогонки

Неявная разностная схема параболического уравнения содержит в каждом уравнении три неизвестных значения нового слоя $u_{j-1}^{n+1},u_{j}^{n+1},u_{j+1}^{n+1}$, поэтому, в отличие от явной схемы, выразить значение в узле напрямую нельзя. Для внутренних узлов получается система с трёхдиагональной матрицей, которую решает метод прогонки (учебник 4.2.2).

Шаг 1. Приведение к виду, удобному для прогонки

Группируем в левой части члены, содержащие значения на (n+1)-м слое, в правой — все остальные:

$$-\frac{\sigma\Delta t}{h^{2}}u_{j+1}^{n+1}+\left(1+\frac{2\sigma\Delta t}{h^{2}}\right)u_{j}^{n+1}-\frac{\sigma\Delta t}{h^{2}}u_{j-1}^{n+1}=u_{j}^{n}+\Delta t\,f\left(t^{n},x_{j}\right)$$

Вводим обозначения для коэффициентов трёхдиагональной системы:

$$a_{j}=-\frac{\sigma\Delta t}{h^{2}};\quad b_{j}=1+\frac{2\sigma\Delta t}{h^{2}};\quad c_{j}=-\frac{\sigma\Delta t}{h^{2}};\quad\xi_{j}^{n}=u_{j}^{n}+\Delta t\,f\left(t^{n},x_{j}\right)$$

и получаем компактную запись:

$$a_{j}u_{j+1}^{n+1}+b_{j}u_{j}^{n+1}+c_{j}u_{j-1}^{n+1}=\xi_{j}^{n}$$

Шаг 2. Рекуррентное прогоночное соотношение

Чтобы связать неизвестные между собой, вводят дополнительное условие в виде линейной зависимости, справедливой для всех $j=1,\dots,N-1$:

$$u_{j}^{n+1}=\alpha_{j}u_{j+1}^{n+1}+\beta_{j}$$

Это рекуррентное прогоночное соотношение, а $\alpha_{j},\beta_{j}$ — прогоночные коэффициенты.

Шаг 3. Вывод формул для α_j, β_j

Запишем то же соотношение для индекса $(j-1)$:

$$u_{j-1}^{n+1}=\alpha_{j-1}u_{j}^{n+1}+\beta_{j-1}$$

и подставим его в преобразованную схему вместо $u_{j-1}^{n+1}$:

$$a_{j}u_{j+1}^{n+1}+b_{j}u_{j}^{n+1}+c_{j}\alpha_{j-1}u_{j}^{n+1}+c_{j}\beta_{j-1}=\xi_{j}^{n}$$

Выражаем $u_{j}^{n+1}$, собирая коэффициент при нём $\left(b_{j}+c_{j}\alpha_{j-1}\right)$:

$$u_{j}^{n+1}=-\frac{a_{j}}{b_{j}+c_{j}\alpha_{j-1}}\,u_{j+1}^{n+1}+\frac{\xi_{j}^{n}-c_{j}\beta_{j-1}}{b_{j}+c_{j}\alpha_{j-1}}$$

Сравнивая это выражение с рекуррентным прогоночным соотношением $u_{j}^{n+1}=\alpha_{j}u_{j+1}^{n+1}+\beta_{j}$, по членам получаем формулы для прогоночных коэффициентов:

$$\boxed{\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}}}$$

Эти формулы позволяют рассчитать коэффициенты на j-м шаге по координате, если известны их значения на $(j-1)$-м шаге. Это и есть прямой ход прогонки (слева направо). Чтобы запустить рекурсию, нужны начальные значения $\alpha_{1},\beta_{1}$.

Шаг 4. Начальные коэффициенты α₁, β₁ (из левого ГУ)

Запишем прогоночное соотношение для $j=1$: $u_{1}^{n+1}=\alpha_{1}u_{2}^{n+1}+\beta_{1}$ и сопоставим с левым граничным условием. Для ГУ 1-го рода $u_{1}^{n+1}=\varphi_{1}(t^{n+1})$ сравнение даёт:

$$\alpha_{1}=0,\qquad\beta_{1}=\varphi_{1}\left(t^{n+1}\right)$$

Для ГУ 2-го рода $\dfrac{\partial u}{\partial x}(t,a)=\varphi_{1}(t)$ аппроксимация $\dfrac{u_{2}^{n+1}-u_{1}^{n+1}}{h}=\varphi_{1}(t^{n+1})$ даёт $u_{1}^{n+1}=u_{2}^{n+1}-h\varphi_{1}(t^{n+1})$, откуда:

$$\alpha_{1}=1,\qquad\beta_{1}=-h\varphi_{1}\left(t^{n+1}\right)$$

Для ГУ 3-го рода аналогично получают $\alpha_{1}=\dfrac{1}{1+h\varphi_{1}(t^{n+1})}$, $\beta_{1}=-\dfrac{h\psi_{1}(t^{n+1})}{1+h\varphi_{1}(t^{n+1})}$. Методика определения не меняется, меняются лишь формулы.

Шаг 5. Решение на правой границе и обратный ход

Прогоночное соотношение даёт $u_{j}^{n+1}$, только если известно соседнее справа значение $u_{j+1}^{n+1}$. Поэтому сначала находят крайнюю правую точку. Для ПГУ 1-го рода: $u_{N}^{n+1}=\varphi_{2}(t^{n+1})$. Для ПГУ 2-го рода подставляют $u_{N-1}^{n+1}=\alpha_{N-1}u_{N}^{n+1}+\beta_{N-1}$ в аппроксимацию $\dfrac{u_{N}^{n+1}-u_{N-1}^{n+1}}{h}=\varphi_{2}(t^{n+1})$ и получают:

$$u_{N}^{n+1}=\frac{h\varphi_{2}\left(t^{n+1}\right)+\beta_{N-1}}{1-\alpha_{N-1}}$$

После этого выполняют обратный ход (справа налево, $j=N-1,\dots,1$) по формуле $u_{j}^{n+1}=\alpha_{j}u_{j+1}^{n+1}+\beta_{j}$, рассчитывая все значения нового слоя. Полная блок-схема алгоритма прогонки приведена в учебнике 4.2.6 и на странице Блок-схемы.

Условие устойчивости (сходимости) прогонки — диагональное преобладание

Абсолютная устойчивость самой неявной схемы ещё не гарантирует сходимость её метода решения. Достаточным условием сходимости метода прогонки является диагональное преобладание в трёхдиагональной матрице — модуль диагонального коэффициента должен превосходить сумму модулей внедиагональных:

$$\boxed{|a_{j}|+|c_{j}|<|b_{j}|}$$

Для рассматриваемой неявной схемы это условие выполняется всегда, так как:

$$|a_{j}|+|c_{j}|=\frac{2\sigma\Delta t}{h^{2}}<1+\frac{2\sigma\Delta t}{h^{2}}=|b_{j}|$$

(добавляется единица из коэффициента $b_{j}$). Поэтому метод прогонки для неявной схемы параболического уравнения всегда устойчив. На практике при свободном члене вида $-ku$ коэффициент $b_{j}$ ещё увеличивается на $k\Delta t$, что лишь усиливает диагональное преобладание — это разобрано в Семинаре 4 (пример 3 с прогонкой).

Постановка задачи

Рассматривается одномерное дифференциальное уравнение параболического типа с начальным и граничными условиями (учебник, гл. 4.3.1):

$$\frac{\partial u}{\partial t}=\sigma\frac{\partial^{2} u}{\partial x^{2}}+f(t, x),\quad u(t{=}0,x)=\xi(x);\quad \begin{cases}u(t,x{=}a)=\varphi_1(t)\\ u(t,x{=}b)=\varphi_2(t)\end{cases}$$

Неявной называют разностную схему, в которой пространственные производные аппроксимируются на новом, (n+1)-м шаге по времени, так что схема содержит несколько неизвестных значений функции на новом слое и решается совместно (методом прогонки). Ниже приведены три неявные схемы для этого уравнения: классическая неявная схема и схема Кранка-Николсона (различающиеся порядком аппроксимации по времени), а также упомянута схема Саульева.

1. Классическая неявная разностная схема — первый порядок по времени

Обе производные второго порядка по координате аппроксимируются на (n+1)-м шаге по времени (учебник, гл. 4.2.1):

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

• Порядок аппроксимации: $O(\Delta t,\,h^{2})$ — первый по времени, второй по координате (такой же, как у явной схемы).
• Устойчивость: абсолютно устойчива (доказано спектральным методом в гл. 3.4) — погрешность не возрастает при любом выборе шагов $\Delta t$ и $h$.
• Метод решения: шаблон содержит три неизвестных на $(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$:

$$-\frac{\sigma\Delta t}{h^{2}}u_{j+1}^{n+1}+\Bigl(1+2\frac{\sigma\Delta t}{h^{2}}\Bigr)u_{j}^{n+1}-\frac{\sigma\Delta t}{h^{2}}u_{j-1}^{n+1}=u_{j}^{n}+\Delta t\,f(t^{n},x_j)$$ $$a_{j}=c_{j}=-\frac{\sigma\Delta t}{h^{2}},\quad b_{j}=1+2\frac{\sigma\Delta t}{h^{2}},\quad \xi_{j}^{n}=u_{j}^{n}+\Delta t\,f(t^{n},x_j)$$

Достаточное условие сходимости прогонки $|a_j|+|c_j|<|b_j|$ выполняется при любых $\Delta t,h$:

$$|a_j|+|c_j|=2\frac{\sigma\Delta t}{h^{2}}<1+2\frac{\sigma\Delta t}{h^{2}}=|b_j|$$

2. Разностная схема Кранка-Николсона — второй порядок по времени

Вывод схемы

Идея состоит в том, чтобы вторую производную по координате представить в виде суммы двух половин и аппроксимировать одну половину на n-м шаге, а другую — на (n+1)-м шаге по времени (учебник, гл. 4.3.1):

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

Эта схема называется разностной схемой Кранка-Николсона в честь авторов, создавших её.

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

Деление производной $\partial^{2}u/\partial x^{2}$ пополам с аппроксимацией одной половины на n-м, а другой на (n+1)-м шаге указывает, что в целом эта производная аппроксимируется относительно точки (n+1/2). Тогда конечная разность по времени по отношению к точке $(n+1/2)$ является центральной, имеющей второй порядок аппроксимации. Следовательно, схема аппроксимирует уравнение со вторым порядком и по времени, и по координате:

$$O(\Delta t^{2})+O(h^{2})\quad\text{или}\quad O(\Delta t^{2},h^{2})$$

Порядок аппроксимации схемы Кранка-Николсона выше, чем у явной и неявной схем, т.е. результаты будут более точными.

Устойчивость (спектральный метод)

Отбрасываем $f(t^n,x_j)$ (не влияет на устойчивость) и подставляем гармонику $u_j^n=\lambda^n e^{i\alpha j}$. После деления на $\lambda^n e^{i\alpha j}$ и использования $e^{i\alpha}-2+e^{-i\alpha}=-4\sin^2\tfrac{\alpha}{2}$ получаем:

$$\frac{\lambda-1}{\Delta t}=-\frac{\sigma}{2h^{2}}4\sin^{2}\frac{\alpha}{2}-\frac{\sigma\lambda}{2h^{2}}4\sin^{2}\frac{\alpha}{2}$$ $$\lambda=\frac{1-\dfrac{2\sigma\Delta t}{h^{2}}\sin^{2}\dfrac{\alpha}{2}}{1+\dfrac{2\sigma\Delta t}{h^{2}}\sin^{2}\dfrac{\alpha}{2}}$$

При $\sigma>0$ числитель по абсолютному значению меньше знаменателя, поэтому необходимое условие устойчивости $|\lambda|\le1$ выполняется при любых $\Delta t$ и $h$ — схема Кранка-Николсона абсолютно устойчива.

Метод решения

Шаблон содержит три неизвестных на $(n+1)$-м слое, поэтому схема решается методом прогонки (учебник, гл. 4.3.3). Приведение к прогоночному виду:

$$-\frac{\sigma\Delta t}{2h^{2}}u_{j+1}^{n+1}+\Bigl(1+\frac{\sigma\Delta t}{h^{2}}\Bigr)u_{j}^{n+1}-\frac{\sigma\Delta t}{2h^{2}}u_{j-1}^{n+1}=u_{j}^{n}+\frac{\sigma\Delta t}{2h^{2}}\bigl(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\bigr)+f(t^{n},x_j)\Delta t$$ $$a_{j}=c_{j}=-\frac{\sigma\Delta t}{2h^{2}},\quad b_{j}=1+\frac{\sigma\Delta t}{h^{2}}$$

Достаточное условие сходимости прогонки выполняется при любых $\Delta t,h$:

$$|a_j|+|c_j|=\frac{\sigma\Delta t}{h^{2}}<1+\frac{\sigma\Delta t}{h^{2}}=|b_j|$$

На семинаре 5 показана полная реализация схемы Кранка-Николсона по 8-пунктовому алгоритму: запись схемы, приведение к прогоночному виду, проверка сходимости прогонки, вывод $\alpha_j,\beta_j$, получение $\alpha_1,\beta_1$ из левого ГУ и $U_{N_x}^{n+1}$ из правого ГУ (семинар 5).

3. Разностная схема Саульева — второй порядок по времени, но без прогонки

Ещё одна схема второго порядка по времени — схема Саульева. Вторую производную записывают как $\lambda_{xx}u=\dfrac{\frac{u_{j+1}-u_j}{h}-\frac{u_j-u_{j-1}}{h}}{h}$ и аппроксимируют первую дробь на n-м шаге, а вторую — на (n+1)-м шаге, получая первую ступень; затем аналогично записывают вторую ступень со сдвигом на $(n+2)$ (учебник, гл. 4.4.1):

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

В совокупности обе ступени дают схему Саульева. Рассматривая их полусумму и раскладывая $u_j^{n}, u_j^{n+2}$ в ряд Тейлора относительно $(t^{n+1},x_j)$, получают, что разностный оператор по времени относительно точки $(n+1)$ является центральной разностью — схема имеет порядок $O(\Delta t^{2},h^{2})$. Схема абсолютно устойчива (принимается без доказательства), а её главное достоинство — она не требует прогонки: каждая ступень содержит лишь две неизвестные на новом слое и решается явным рекуррентным соотношением (вперёд для первой ступени, назад для второй). Однако точность зависит от отношения $\Delta t/h$ (хорошие результаты при $\Delta t\sim h^2$), а оценить решение можно только на $(n+2)$-м шаге, поэтому интервал по времени делят на чётное число шагов.

Сравнительная сводка (гл. 4.5)

Схема	Порядок	Устойчивость	Метод решения
Неявная	$O(\Delta t,h^{2})$	абсолютная	прогонка
Кранка-Николсона	$O(\Delta t^{2},h^{2})$	абсолютная	прогонка
Саульева	$O(\Delta t^{2},h^{2})$ при $\Delta t\sim h^{2}$	абсолютная	две явные ступени (без прогонки)

Вывод: в качестве неявной схемы первого порядка по времени используется классическая неявная схема, а для повышения порядка по времени до второго применяют схему Кранка-Николсона (полунеявную, решаемую прогонкой) или схему Саульева (источники: учебник, гл. 4.5; семинар 5).

Постановка задачи

Дифференциальные уравнения в частных производных 1-го порядка часто встречаются в моделях химико-технологических процессов: баланс по концентрации в реакторе идеального вытеснения $\frac{\partial c}{\partial t}+v\frac{\partial c}{\partial x}=-kc$ и уравнение баланса числа частиц при кристаллизации $\frac{\partial f}{\partial t}+\eta\frac{\partial f}{\partial l}=0$. Для простоты их записывают в общем виде (учебник, гл. 5.1):

$$\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}=f(t,x)$$

Параметр $v$ может быть как положительным, так и отрицательным, но не равным нулю (при $v=0$ это уже обыкновенное дифференциальное уравнение). Уравнение дополняется начальным условием $u(t{=}0,x)=\xi(x)$ и одним граничным условием 1-го рода; будет ли оно левым или правым — определяется методом решения.

Четыре разностные схемы (гл. 5.2)

Производную по времени аппроксимируют правой конечной разностью. Производную по координате $\partial u/\partial x$ можно аппроксимировать правой или левой конечной разностью, со стабилизацией значения $t$ либо на n-м (явная), либо на (n+1)-м (неявная) шаге. Это даёт четыре схемы (учебник, гл. 5.2):

1) Явная, правая разность (5.2):

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

2) Явная, левая разность (5.3):

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

3) Неявная, правая разность (5.4):

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

4) Неявная, левая разность (5.5):

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

Все четыре схемы имеют первый порядок аппроксимации и по времени, и по координате: $O(\Delta t,h)$.

Устойчивость (краткие результаты)

Исследование спектральным методом (подробности — в вопросе 12) даёт:

• Явная, правая разность (5.2): условно устойчива при $v<0$, условие $-1\le v\frac{\Delta t}{h}<0$; при $v>0$ — неустойчива.
• Явная, левая разность (5.3): условно устойчива при $v>0$, условие $0
• Неявная, правая разность (5.4): абсолютно устойчива при $v<0$.
• Неявная, левая разность (5.5): абсолютно устойчива при $v>0$.

Методы решения

Явные схемы

Шаблон содержит одну неизвестную на $(n+1)$-м слое, поэтому она прямо выражается рекуррентным соотношением.

Для правой разности (5.2) — значения считаются справа налево, а на правой границе берётся правое ГУ $u_{N}^{n+1}=\varphi_2(t^{n+1})$:

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

Для левой разности (5.3) — слева направо, а на левой границе берётся левое ГУ $u_{1}^{n+1}=\varphi_1(t^{n+1})$:

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

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

Шаблон содержит две неизвестные на $(n+1)$-м слое, но прогонка не нужна — значения вычисляются последовательным рекуррентным соотношением, используя уже найденное соседнее значение на новом слое.

Неявная правая разность (5.4) — нужно соседнее справа значение $u_{j+1}^{n+1}$, расчёт ведётся справа налево ($j=N-1,\dots,1$), стартуя с правого ГУ:

$$u_{j}^{n+1}=\frac{u_{j}^{n}+v\frac{\Delta t}{h}u_{j+1}^{n+1}+\Delta t\,f(t^{n},x_j)}{1+v\frac{\Delta t}{h}},\qquad u_{N}^{n+1}=\varphi_2(t^{n+1})$$

Неявная левая разность (5.5) — нужно соседнее слева значение $u_{j-1}^{n+1}$, расчёт ведётся слева направо ($j=2,\dots,N$), стартуя с левого ГУ:

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

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

Правило выбора схемы по знаку параметра $v$

Определяющим фактором является знак параметра $v$ при производной первого порядка (учебник, гл. 5.8):

Чтобы разностная схема была устойчива (условно — для явной, абсолютно — для неявной), при положительном $v$ для аппроксимации первой производной по координате следует использовать левую конечную разность, при отрицательном $v$ — правую конечную разность. Кроме того, для решения схемы при $v>0$ потребуется левое граничное условие, при $v<0$ — правое.

Геометрически это соответствует направлению переноса (потока): разность берётся «против потока» (с той стороны, откуда приходит информация), и граничное условие задаётся на входной границе.

Важное замечание: правило применимо, только если производная первого порядка по координате находится в левой части уравнения (вид (5.1)). Если она в правой части, её необходимо сначала перенести в левую и только потом применять правило.

Знак $v$	Разность по $x$	Граничное условие	Рекомендуемая схема
$v>0$	левая	левое	неявная с левой разностью (5.5), абс. устойчива
$v<0$	правая	правое	неявная с правой разностью (5.4), абс. устойчива

Практический разбор всех случаев (явные/неявные схемы с левой/правой разностью, выбор по знаку $v$, нахождение максимального $\Delta t$, блок-схемы) приведён на семинаре 6.

Метод доказательства — спектральный (метод гармоник)

Устойчивость исследуется спектральным методом для уравнения $\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x}=f(t,x)$. Свободный член $f(t^n,x_j)$ отбрасывается (он не влияет на устойчивость), а решение представляется в виде гармоники $u_j^n=\lambda^n e^{i\alpha j}$. Необходимое условие устойчивости — модуль собственного числа оператора перехода не превосходит единицы:

$$|\lambda|\le1$$

Так как $\lambda$ — комплексное число, условие $|\lambda|\le1$ означает, что собственные числа должны лежать внутри или на границе круга радиуса 1 с центром в начале координат комплексной плоскости. Доказательство сводится к геометрическому анализу: для каждой схемы получают $\lambda$ (или $1/\lambda$) как точку на некоторой окружности и сравнивают её положение с единичным кругом (учебник, гл. 5.3.1).

1. Явная схема с правой разностью (5.2)

Подставляя гармонику и деля на $\lambda^n e^{i\alpha j}$:

$$\frac{\lambda-1}{\Delta t}+v\frac{e^{i\alpha}-1}{h}=0\quad\Rightarrow\quad \lambda=1+v\frac{\Delta t}{h}-v\frac{\Delta t}{h}e^{i\alpha}$$

Случай $v<0$. Обозначим $r=-v\frac{\Delta t}{h}>0$, тогда $\lambda=1-r+re^{i\alpha}$. Это окружность с центром в точке $(1-r,0)$ и радиусом $|re^{i\alpha}|=\sqrt{r^2\cos^2\alpha+r^2\sin^2\alpha}=r$. Сравнивая с единичным кругом, получаем три случая:

• при $r<1$ — окружность целиком внутри единичного круга (устойчиво);
• при $r=1$ — совпадает с границей (на пределе устойчивости);
• при $r>1$ — выходит за пределы круга (неустойчиво).

Следовательно, при $v<0$ схема устойчива при $r=-v\frac{\Delta t}{h}\le1$.

Случай $v>0$. Обозначим $q=v\frac{\Delta t}{h}>0$, тогда $\lambda=1+q-qe^{i\alpha}$ — окружность с центром в $(1+q,0)$ и радиусом $q$. Эта окружность находится вне единичного круга при любом $q$, поэтому при $v>0$ схема неустойчива.

Вывод: явная схема с правой разностью условно устойчива, условие $-1\le v\frac{\Delta t}{h}<0$ (т.е. только при $v<0$).

2. Явная схема с левой разностью (5.3)

Аналогично (учебник, гл. 5.4.1):

$$\frac{\lambda-1}{\Delta t}+v\frac{1-e^{-i\alpha}}{h}=0\quad\Rightarrow\quad \lambda=1-v\frac{\Delta t}{h}+v\frac{\Delta t}{h}e^{-i\alpha}$$

Случай $v<0$. $q=-v\frac{\Delta t}{h}>0$, $\lambda=1+q-qe^{-i\alpha}$ — окружность с центром $(1+q,0)$, радиусом $q$, лежит вне единичного круга при любом $q$ → неустойчиво.

Случай $v>0$. $r=v\frac{\Delta t}{h}>0$, $\lambda=1-r+re^{-i\alpha}$ — окружность с центром $(1-r,0)$, радиусом $r$. Те же три случая: при $r<1$ внутри, при $r=1$ на границе, при $r>1$ вне круга. Условие устойчивости $r=v\frac{\Delta t}{h}\le1$.

Вывод: явная схема с левой разностью условно устойчива при $00$).

3. Неявная схема с правой разностью (5.4)

Здесь удобнее выражать величину, обратную $\lambda$ (учебник, гл. 5.5.1):

$$\frac{\lambda-1}{\Delta t}+v\frac{\lambda e^{i\alpha}-\lambda}{h}=0\quad\Rightarrow\quad \frac{1}{\lambda}=1-v\frac{\Delta t}{h}+v\frac{\Delta t}{h}e^{i\alpha}$$

Условие устойчивости преобразуется: $|\lambda|\le1\;\Leftrightarrow\;\bigl|\tfrac{1}{\lambda}\bigr|\ge1$. То есть величины $1/\lambda$ должны лежать вне или на границе единичного круга.

Случай $v<0$. $r=-v\frac{\Delta t}{h}>0$, $\frac{1}{\lambda}=1+r-re^{i\alpha}$ — окружность с центром $(1+r,0)$, радиусом $r$. Она находится вне единичного круга при любом $r$, поэтому условие $|1/\lambda|\ge1$ выполняется всегда — схема абсолютно устойчива при $v<0$.

4. Неявная схема с левой разностью (5.5)

Аналогично (учебник, гл. 5.6.1):

$$\frac{\lambda-1}{\Delta t}+v\frac{\lambda-\lambda e^{-i\alpha}}{h}=0\quad\Rightarrow\quad \frac{1}{\lambda}=1+v\frac{\Delta t}{h}-v\frac{\Delta t}{h}e^{-i\alpha}$$

Случай $v>0$. $r=v\frac{\Delta t}{h}>0$, $\frac{1}{\lambda}=1+r-re^{-i\alpha}$ — окружность с центром $(1+r,0)$, радиусом $r$, лежит вне единичного круга при любом $r$ → условие $|1/\lambda|\ge1$ выполняется всегда. Схема абсолютно устойчива при $v>0$.

Геометрическая суть доказательства

Все четыре доказательства опираются на один приём: собственное число (или обратное к нему) представляется в виде $\text{центр}\pm r\,e^{\pm i\alpha}$, что есть точка на окружности радиуса $r$ с центром на вещественной оси. Устойчивость сводится к взаимному расположению этой окружности и единичного круга:

• для явных схем требуется $|\lambda|\le1$ → окружность $\lambda$ должна лежать внутри единичного круга (получается только при правильном выборе разности и при ограничении на $\Delta t$, т.е. условная устойчивость);
• для неявных схем требуется $|1/\lambda|\ge1$ → окружность $1/\lambda$ должна лежать вне единичного круга; центр в точке $(1+r,0)$ всегда смещён вправо настолько, что окружность целиком вне круга при любом $r$ → абсолютная устойчивость.

Итог: устойчивость определяется знаком $v$ и выбором конечной разности «против потока». Геометрический анализ устойчивости с тремя случаями расположения окружностей ($r<1$, $r=1$, $r>1$) на комплексной плоскости подробно разобран на семинаре 6 (где также показано, что неверный выбор разности — например, правая при $v>0$ — даёт окружность с центром $(1+r,0)$ вне единичного круга, т.е. неустойчивую явную схему).

Постановка задачи

Рассматривается двумерное дифференциальное уравнение параболического типа, не содержащее первых производных по координатам (учебник, гл. 7.2):

$$\frac{\partial u}{\partial t}=\sigma\Bigl(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\Bigr)+f(t,x,y),\qquad \sigma>0$$

с начальным условием $u(t{=}0,x,y)=\xi(x,y)$ и двумя граничными условиями по каждой из координат $x$ и $y$. Вводится трёхмерная разностная сетка: $n$ — номер по $t$, $j$ — по $x$, $k$ — по $y$; шаги $\Delta t$, $h_x$, $h_y$; $u_{j,k}^{n}=u(t^n,x_j,y_k)$.

Построение явной схемы (гл. 7.3)

Производная по времени аппроксимируется правой конечной разностью, обе вторые производные по координатам — на n-м (старом) шаге по времени разностным оператором второго порядка:

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

Подставляя в уравнение, получаем явную разностную схему:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=\sigma\Bigl(\frac{u_{j+1,k}^{n}-2u_{j,k}^{n}+u_{j-1,k}^{n}}{h_x^{2}}+\frac{u_{j,k+1}^{n}-2u_{j,k}^{n}+u_{j,k-1}^{n}}{h_y^{2}}\Bigr)+f_{j,k}^{n}$$

Порядок аппроксимации: $O(\Delta t,h_x^{2},h_y^{2})$ — первый по времени, второй по каждой координате.

Исследование устойчивости (гл. 7.4.1)

Применяем спектральный метод. Отбрасываем $f_{j,k}^{n}$ и представляем решение в виде двумерной гармоники (учебник, гл. 7.4.1):

$$u_{j,k}^{n}=\lambda^{n}e^{i\alpha j}e^{i\beta k},\qquad \alpha\in[0,2\pi],\;\beta\in[0,2\pi]$$

Подставляя и деля на $\lambda^{n}e^{i\alpha j}e^{i\beta k}$:

$$\frac{\lambda-1}{\Delta t}=\sigma\Bigl(\frac{e^{i\alpha}-2+e^{-i\alpha}}{h_x^{2}}+\frac{e^{i\beta}-2+e^{-i\beta}}{h_y^{2}}\Bigr)$$

Используя $e^{i\alpha}-2+e^{-i\alpha}=-4\sin^2\tfrac{\alpha}{2}$ (и аналогично по $\beta$), получаем:

$$\frac{\lambda-1}{\Delta t}=-\frac{\sigma}{h_x^{2}}4\sin^{2}\frac{\alpha}{2}-\frac{\sigma}{h_y^{2}}4\sin^{2}\frac{\beta}{2}$$ $$\lambda=1-4\sigma\frac{\Delta t}{h_x^{2}}\sin^{2}\frac{\alpha}{2}-4\sigma\frac{\Delta t}{h_y^{2}}\sin^{2}\frac{\beta}{2}$$

Применяем необходимое условие устойчивости $|\lambda|\le1$:

$$-1\le 1-4\sigma\frac{\Delta t}{h_x^{2}}\sin^{2}\frac{\alpha}{2}-4\sigma\frac{\Delta t}{h_y^{2}}\sin^{2}\frac{\beta}{2}\le1$$

Правое неравенство выполняется автоматически (вычитаются неотрицательные величины). Рассматриваем левое:

$$1-4\sigma\frac{\Delta t}{h_x^{2}}\sin^{2}\frac{\alpha}{2}-4\sigma\frac{\Delta t}{h_y^{2}}\sin^{2}\frac{\beta}{2}\ge-1\;\Rightarrow\;\frac{\Delta t}{h_x^{2}}\sin^{2}\frac{\alpha}{2}+\frac{\Delta t}{h_y^{2}}\sin^{2}\frac{\beta}{2}\le\frac{1}{2\sigma}$$

Задавая для $\sin^2\tfrac{\alpha}{2}$ и $\sin^2\tfrac{\beta}{2}$ максимальное значение, равное 1, переходим к более строгому условию, справедливому для любых $\alpha,\beta$.

Условие устойчивости

Общее условие устойчивости явной схемы для двумерного уравнения (формула (7.4)):

$$\boxed{\;\frac{\Delta t}{h_x^{2}}+\frac{\Delta t}{h_y^{2}}\le\frac{1}{2\sigma}\quad\Leftrightarrow\quad \Delta t\le\frac{1}{\dfrac{2\sigma}{h_x^{2}}+\dfrac{2\sigma}{h_y^{2}}}\;}$$

Если шаги по координатам равны $h_x=h_y=h$, условие принимает простой вид:

$$\frac{\Delta t}{h^{2}}\le\frac{1}{4\sigma}\qquad\Bigl(\text{т.е. }\Delta t\le\frac{h^{2}}{4\sigma}\Bigr)$$

Сравнение с одномерным случаем. Для одномерного уравнения условие устойчивости явной схемы было $\frac{\Delta t}{h^{2}}\le\frac{1}{2\sigma}$. Видно, что увеличение размерности на порядок (с 1D до 2D) приводит к уменьшению вдвое максимально допустимого значения $\Delta t$, при котором явная схема остаётся устойчивой. Таким образом, явная схема для двумерного уравнения является условно устойчивой.

Метод решения (гл. 7.4.2)

Шаблон явной схемы содержит лишь одну неизвестную величину на $(n+1)$-м слое, поэтому она прямо выражается рекуррентным соотношением:

$$u_{j,k}^{n+1}=u_{j,k}^{n}+\sigma\frac{\Delta t}{h_x^{2}}\bigl(u_{j+1,k}^{n}-2u_{j,k}^{n}+u_{j-1,k}^{n}\bigr)+\sigma\frac{\Delta t}{h_y^{2}}\bigl(u_{j,k+1}^{n}-2u_{j,k}^{n}+u_{j,k-1}^{n}\bigr)+\Delta t\,f_{j,k}^{n}$$

или в компактной операторной форме:

$$u_{j,k}^{n+1}=u_{j,k}^{n}+\sigma\Delta t\,\lambda_{xx}u_{j,k}^{n}+\sigma\Delta t\,\lambda_{yy}u_{j,k}^{n}+\Delta t\,f_{j,k}^{n}$$

Это соотношение позволяет рассчитать все внутренние значения на $(n+1)$-м слое по известным значениям на n-м слое. Граничные значения $u_{1,k}^{n+1}$, $u_{N_x,k}^{n+1}$, $u_{j,1}^{n+1}$, $u_{j,N_y}^{n+1}$ определяются из граничных условий (для ГУ 1-го рода — непосредственно):

$$\begin{cases}u_{1,k}^{n+1}=\varphi_1(t^{n+1},y_k)\\ u_{N_x,k}^{n+1}=\varphi_2(t^{n+1},y_k)\end{cases}\qquad \begin{cases}u_{j,1}^{n+1}=\psi_1(t^{n+1},x_j)\\ u_{j,N_y}^{n+1}=\psi_2(t^{n+1},x_j)\end{cases}$$

Расчёт организуется двойным циклом по $j$ и $k$ для всех внутренних узлов на каждом шаге по времени (алгоритм — учебник, гл. 7.4.3).

Практический разбор устойчивости явной схемы в 2D (с условием $\Delta t\le\frac{h^2}{4\sigma}$ при $h_x=h_y=h$) и пример её реализации приведены на семинаре 8.

Вопрос 14. Схема расщепления (метод дробных шагов) для двумерного уравнения параболического типа

1. Постановка задачи и исходная неявная схема

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

$$\frac{\partial u}{\partial t}=\sigma\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)+f(t,x,y);\quad \sigma>0 \tag{7.1}$$

Аппроксимация обеих вторых производных на $(n+1)$-м шаге даёт неявную разностную схему:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=\sigma\left(\frac{u_{j+1,k}^{n+1}-2u_{j,k}^{n+1}+u_{j-1,k}^{n+1}}{h_{x}^{2}}+\frac{u_{j,k+1}^{n+1}-2u_{j,k}^{n+1}+u_{j,k-1}^{n+1}}{h_{y}^{2}}\right)+f_{j,k}^{n} \tag{7.3}$$

Эта схема абсолютно устойчива (доказывается спектральным методом — собственные числа оператора перехода удовлетворяют необходимому условию устойчивости при любых $\Delta t, h_x, h_y$), но её разностный шаблон содержит пять неизвестных значений функции на $(n+1)$-м шаге. Поэтому без дополнительных преобразований она неразрешима (метод прогонки в чистом виде неприменим). Подробнее: 7.5. Характеристика неявной разностной схемы.

2. Методика записи уравнений схемы (метод дробных шагов)

Метод дробных шагов позволяет представить неявную схему (7.3) в виде двух подсхем, каждая из которых решается методом прогонки. Для этого интервал $\Delta t$ между точками $t^{n}$ и $t^{n+1}$ делится пополам; промежуточная точка обозначается $t^{n+1/2}$ (рис. 7.5).

Первая подсхема записывается на первом полушаге интервала $\Delta t$ и учитывает только производную второго порядка по координате $x$ (неявная по $x$):

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}=\sigma\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+f_{j,k}^{n} \tag{7.7}$$

Вторая подсхема записывается на втором полушаге $\Delta t$ и учитывает только производную второго порядка по координате $y$ (неявная по $y$):

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}=\sigma\frac{u_{j,k+1}^{n+1}-2u_{j,k}^{n+1}+u_{j,k-1}^{n+1}}{h_{y}^{2}} \tag{7.8}$$

Подробный вывод: 7.6. Схема расщепления.

3. Проверка аппроксимации (сложение подсхем) и порядок аппроксимации

Складывая подсхемы (7.7) и (7.8), получаем соотношение, отличающееся от неявной схемы (7.3) лишь тем, что вторая производная по $x$ аппроксимирована не на $(n+1)$-м шаге, а на шаге $(n+1/2)$:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=\sigma\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+\sigma\frac{u_{j,k+1}^{n+1}-2u_{j,k}^{n+1}+u_{j,k-1}^{n+1}}{h_{y}^{2}}+f_{j,k}^{n}$$

Следовательно, исходное уравнение (7.1) аппроксимируется последовательным решением двух подсхем, называемых в совокупности схемой расщепления. Порядок аппроксимации схемы расщепления такой же, как у неявной схемы (7.3) — первый по времени и второй по каждой из координат:

$$O\left(\Delta t,\,h_{x}^{2},\,h_{y}^{2}\right)$$

Свободный член $f$ может быть учтён не в первой, а во второй подсхеме — его вид при этом не изменится.

4. Характеристика первой подсхемы и метод решения (прогонка по $x$)

Первая подсхема (7.7) — аналог одномерной неявной схемы: абсолютно устойчива, решается методом прогонки. Приводим её к трёхдиагональному виду (4.10):

$$-\sigma\frac{\Delta t}{h_{x}^{2}}u_{j+1,k}^{n+1/2}+\left(1+2\sigma\frac{\Delta t}{h_{x}^{2}}\right)u_{j,k}^{n+1/2}-\sigma\frac{\Delta t}{h_{x}^{2}}u_{j-1,k}^{n+1/2}=u_{j,k}^{n}+\Delta t\,f_{j,k}^{n}$$

Коэффициенты:

$$a_{j}=c_{j}=-\sigma\frac{\Delta t}{h_{x}^{2}},\quad b_{j}=1+2\sigma\frac{\Delta t}{h_{x}^{2}},\quad \xi_{j,k}^{n}=u_{j,k}^{n}+\Delta t\,f_{j,k}^{n}$$

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

$$|a_{j}|+|c_{j}|=2\sigma\frac{\Delta t}{h_{x}^{2}}<1+2\sigma\frac{\Delta t}{h_{x}^{2}}=|b_{j}|$$

Рекуррентное прогоночное соотношение и прогоночные коэффициенты:

$$u_{j,k}^{n+1/2}=\alpha_{j}\,u_{j+1,k}^{n+1/2}+\beta_{j} \tag{7.9}$$

$$\alpha_{j}=-\frac{a_{j}}{b_{j}+c_{j}\alpha_{j-1}},\quad \beta_{j}=\frac{\xi_{j,k}^{n}-c_{j}\beta_{j-1}}{b_{j}+c_{j}\alpha_{j-1}} \tag{7.10}$$

Прогонка ведётся по координате $x$ (индекс $j$); коэффициенты $\alpha_1,\beta_1$ и решение на правой границе берутся из граничных условий по $x$. Так как соотношения содержат переменную $k$, задаётся внешний цикл $k=2,\ldots,N_{y}-1$, т.е. на первом полушаге прогонка применяется $N_{y}-2$ раза. Результат — значения $u^{n+1/2}_{j,k}$, нужные для второй подсхемы. (Так как по отдельности подсхемы не аппроксимируют (7.1), оценить погрешность $u^{n+1/2}$ невозможно — близость к истинным значениям гарантируется лишь для $u^{n+1}$.) См. 7.6.2. Характеристика первой подсхемы.

5. Характеристика второй подсхемы и метод решения (прогонка по $y$)

Вторая подсхема (7.8) также абсолютно устойчива и решается прогонкой. Трёхдиагональный вид:

$$-\sigma\frac{\Delta t}{h_{y}^{2}}u_{j,k+1}^{n+1}+\left(1+2\sigma\frac{\Delta t}{h_{y}^{2}}\right)u_{j,k}^{n+1}-\sigma\frac{\Delta t}{h_{y}^{2}}u_{j,k-1}^{n+1}=u_{j,k}^{n+1/2}$$

Коэффициенты:

$$\tilde a_{k}=\tilde c_{k}=-\sigma\frac{\Delta t}{h_{y}^{2}},\quad \tilde b_{k}=1+2\sigma\frac{\Delta t}{h_{y}^{2}},\quad \tilde\xi_{j,k}^{n+1/2}=u_{j,k}^{n+1/2}$$

Условие сходимости: $|\tilde a_{k}|+|\tilde c_{k}|=2\sigma\frac{\Delta t}{h_{y}^{2}}<|\tilde b_{k}|$. Рекуррентное соотношение и коэффициенты:

$$u_{j,k}^{n+1}=\widetilde\alpha_{k}\,u_{j,k+1}^{n+1}+\widetilde\beta_{k} \tag{7.11}$$

$$\widetilde\alpha_{k}=-\frac{\tilde a_{k}}{\tilde b_{k}+\tilde c_{k}\widetilde\alpha_{k-1}},\quad \widetilde\beta_{k}=\frac{\tilde\xi_{j,k}^{n+1/2}-\tilde c_{k}\widetilde\beta_{k-1}}{\tilde b_{k}+\tilde c_{k}\widetilde\alpha_{k-1}} \tag{7.12}$$

Прогонка ведётся по координате $y$ (индекс $k$); $\widetilde\alpha_1,\widetilde\beta_1$ и решение на правой границе берутся из граничных условий по $y$. Внешний цикл $j=2,\ldots,N_{x}-1$ — прогонка применяется $N_{x}-2$ раза. Результат — значения $u(t,x,y)$ на $(n+1)$-м шаге. См. 7.6.3. Характеристика второй подсхемы.

6. Алгоритм решения

1. Задание начальных условий: цикл по $j$ и $k$, $u_{j,k}^{0}=\xi(x_j,y_k)$.
2. Цикл по времени $n$.
3. Первый полушаг (по $x$): постановка граничных условий по $x$ для шага $n+1/2$; внешний цикл по $k=2,\ldots,N_y-1$; для каждого $k$ — прямой и обратный ход прогонки (7.9)-(7.10) по $j$ → получаем $u_{j,k}^{n+1/2}$.
4. Второй полушаг (по $y$): постановка граничных условий по $y$ для шага $n+1$; внешний цикл по $j=2,\ldots,N_x-1$; для каждого $j$ — прогонка (7.11)-(7.12) по $k$ → получаем $u_{j,k}^{n+1}$.
5. Переход к следующему шагу по времени.

Блок-схема приведена на рис. 7.6 (см. 7.6.4. Алгоритм решения и страницу с блок-схемами). Схема расщепления — наиболее простой способ интерпретации неявной схемы (7.3).

7. Разбор на семинаре

На семинаре 8 (Пример 2) метод расщепления разобран на уравнении $\frac{\partial u}{\partial t}+\frac{\partial u}{\partial y}=0{,}7\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+1$. Подшаг 1 (прогонка по $x$): приведение к виду $u_{j+1,k}^{n+1/2}\!\left[\frac{-0{,}7\Delta t}{h_x^2}\right]+u_{j,k}^{n+1/2}\!\left[1+\frac{1{,}4\Delta t}{h_x^2}+\frac{\Delta t}{h_x}\right]+u_{j-1,k}^{n+1/2}\!\left[\frac{-\Delta t}{h_x}-\frac{0{,}7\Delta t}{h_x^2}\right]=u_{j,k}^{n}+1$, проверка $|a_j|+|c_j|\le|b_j|$, поиск $\alpha_1,\beta_1$ из ЛГУ по $x$. Подшаг 2 (прогонка по $y$): аналогично с коэффициентами по $h_y$ и определением $\widetilde\alpha_1,\widetilde\beta_1$ из граничного условия 2-го рода по $y$.

Итог

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

Вопрос 15. Схема переменных направлений для двумерного уравнения параболического типа

1. Назначение схемы

Схема переменных направлений — это ещё один способ интерпретации абсолютно устойчивой, но неразрешимой напрямую неявной разностной схемы (7.3), аппроксимирующей уравнение

$$\frac{\partial u}{\partial t}=\sigma\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)+f(t,x,y);\quad \sigma>0 \tag{7.1}$$

В отличие от схемы расщепления, она позволяет повысить порядок аппроксимации по времени (с первого до второго). Подробно: 7.7. Схема переменных направлений.

2. Методика записи уравнений схемы

Интервал $\Delta t$ делится пополам точкой $t^{n+1/2}$. Схема записывается в виде двух подсхем:

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}=\frac{\sigma}{2}\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+\frac{\sigma}{2}\frac{u_{j,k+1}^{n}-2u_{j,k}^{n}+u_{j,k-1}^{n}}{h_{y}^{2}}$$

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}=\frac{\sigma}{2}\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+\frac{\sigma}{2}\frac{u_{j,k+1}^{n+1}-2u_{j,k}^{n+1}+u_{j,k-1}^{n+1}+f_{j,k}^{n+1/2}}{h_{y}^{2}} \tag{7.13}$$

Первая подсхема аппроксимирует производную по времени на первом полушаге $\Delta t$, является неявной по координате $x$ (уровень $n+1/2$) и явной по координате $y$ (уровень $n$). Вторая подсхема аппроксимирует производную по времени на втором полушаге $\Delta t$, является неявной по координате $y$ (уровень $n+1$) и явной по координате $x$ (уровень $n+1/2$). Каждая из подсхем (как и в схеме расщепления) абсолютно устойчива и решается методом прогонки.

Две особенности при записи схемы (7.13), которые нужно обязательно учитывать:

1. коэффициенты перед разностными операторами, аппроксимирующими $\partial^2 u/\partial x^2$ и $\partial^2 u/\partial y^2$, делятся пополам (множитель $\sigma/2$);
2. свободный член записывается во второй подсхеме и аппроксимируется на шаге $(n+1/2)$.

3. Проверка аппроксимации и порядок

Складывая обе подсхемы и используя обозначения $\lambda_{xx}u_{j,k}=\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{h_x^2}$, $\lambda_{yy}u_{j,k}=\frac{u_{j,k+1}-2u_{j,k}+u_{j,k-1}}{h_y^2}$ (7.6), получаем:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=\sigma\lambda_{xx}u_{j,k}^{n+1/2}+\frac{\sigma}{2}\left(\lambda_{yy}u_{j,k}^{n+1}+\lambda_{yy}u_{j,k}^{n}\right)+f_{j,k}^{n+1/2}$$

Правая часть аппроксимирована относительно точки $t^{n+1/2}$. Значит, разностный оператор по времени в левой части является центральной конечной разностью, имеющей второй порядок аппроксимации. Поэтому схема переменных направлений имеет порядок

$$O\left(\Delta t^{2},\,h_{x}^{2},\,h_{y}^{2}\right)$$

и является более точной по сравнению со схемой расщепления (которая даёт лишь $O(\Delta t, h_x^2, h_y^2)$).

4. Метод решения (прогонки по направлениям)

Алгоритм решения аналогичен схеме расщепления: первый полушаг — прогонка по $x$, второй полушаг — прогонка по $y$. Коэффициенты трёхдиагонального вида (4.10):

Первая подсхема (прогонка по $x$, индекс $j$):

$$a_{j}=c_{j}=-\frac{\sigma}{2}\frac{\Delta t}{h_{x}^{2}},\quad b_{j}=1+\sigma\frac{\Delta t}{h_{x}^{2}},\quad \xi_{j,k}^{n}=u_{j,k}^{n}+\frac{\sigma}{2}\Delta t\,\lambda_{yy}u_{j,k}^{n}$$

В правую часть входит явно вычисляемый по известному $n$-му слою член $\frac{\sigma}{2}\Delta t\,\lambda_{yy}u_{j,k}^{n}$.

Вторая подсхема (прогонка по $y$, индекс $k$):

$$\tilde a_{k}=\tilde c_{k}=-\frac{\sigma}{2}\frac{\Delta t}{h_{y}^{2}},\quad \tilde b_{k}=1+\sigma\frac{\Delta t}{h_{y}^{2}},$$

$$\tilde\xi_{j,k}^{n+1/2}=u_{j,k}^{n+1/2}+\frac{\sigma}{2}\Delta t\,\lambda_{xx}u_{j,k}^{n+1/2}+\Delta t\,f_{j,k}^{n+1/2}$$

Для обеих подсхем достаточное условие сходимости прогонки (диагональное преобладание) выполняется, так как $|a_j|+|c_j|=\sigma\frac{\Delta t}{h_x^2}<1+\sigma\frac{\Delta t}{h_x^2}=|b_j|$ (аналогично для $y$). Прогоночные соотношения и коэффициенты — те же, что в схеме расщепления (7.9)-(7.12); $\alpha_1,\beta_1$ и решение на правых границах определяются из граничных условий по соответствующей координате.

5. Алгоритм решения

1. Задание начальных условий $u_{j,k}^{0}=\xi(x_j,y_k)$.
2. Цикл по времени $n$.
3. Первый полушаг: ГУ по $x$ на $(n+1/2)$; для каждого $k=2,\ldots,N_y-1$ — вычисление правой части $\xi_{j,k}^{n}$ (с явным членом по $y$) и прогонка по $j$ → $u_{j,k}^{n+1/2}$.
4. Второй полушаг: ГУ по $y$ на $(n+1)$; для каждого $j=2,\ldots,N_x-1$ — вычисление правой части $\tilde\xi_{j,k}^{n+1/2}$ (с явным членом по $x$ и свободным членом $f^{n+1/2}$) и прогонка по $k$ → $u_{j,k}^{n+1}$.
5. Переход к следующему шагу по времени.

Блок-схемы — на странице алгоритмов.

6. Разбор на семинаре

На семинаре 9 схема разобрана подробно. Структура подсхем: ① неявная по $x$ (на уровне $n+1/2$), явная по $y$ (на $n$); ② неявная по $y$ (на $n+1$), явная по $x$ (на $n+1/2$). Пример 1: $\frac{\partial u}{\partial t}=10\frac{\partial^2 u}{\partial x^2}+8\frac{\partial^2 u}{\partial y^2}+t+xy$. Подсхема ① в прогоночном виде: $u_{j+1,k}^{n+1/2}\!\left[\frac{-5\Delta t}{h_x^2}\right]+u_{jk}^{n+1/2}\!\left[1+\frac{10\Delta t}{h_x^2}\right]+u_{j-1,k}^{n+1/2}\!\left[\frac{-5\Delta t}{h_x^2}\right]=u_{jk}^{n}+4\Delta t\,\Lambda_{yy}^{n}$; проверка сходимости $\left|\frac{-5\Delta t}{h_x^2}\right|+\left|\frac{-5\Delta t}{h_x^2}\right|<\left|1+\frac{10\Delta t}{h_x^2}\right|$. Подсхема ②: $u_{jk+1}^{n+1}\!\left[\frac{-4\Delta t}{h_y^2}\right]+u_{jk}^{n+1}\!\left[1+\frac{8\Delta t}{h_y^2}\right]+u_{jk-1}^{n+1}\!\left[\frac{-4\Delta t}{h_y^2}\right]=u_{jk}^{n+1/2}+5\Delta t\,\Lambda_{xx}^{n+1/2}+\Delta t\,f^{n+1/2}_{jk}$. Видно, что коэффициент перед $\Lambda$ делится пополам, а свободный член входит во вторую подсхему на уровне $n+1/2$.

Итог

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

Вопрос 16. Схема предиктор-корректор для двумерного уравнения параболического типа

1. Назначение схемы

Схема предиктор-корректор — ещё одна интерпретация неявной разностной схемы (7.3), которая (как и схема переменных направлений) позволяет повысить порядок аппроксимации по времени до второго. Это наиболее сложная интерпретация неявной схемы, аппроксимирующей уравнение

$$\frac{\partial u}{\partial t}=\sigma\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)+f(t,x,y);\quad \sigma>0 \tag{7.1}$$

Подробно: 7.9.1. Методика записи уравнений схемы.

2. Методика записи уравнений схемы — особое деление интервала

Схема требует особого способа расщепления интервала $\Delta t$ (рис. 7.7): интервал $\Delta t$ между $t^{n}$ и $t^{n+1}$ делится пополам (промежуточная точка $t^{n+1/2}$); затем интервал $\Delta t/2$ между $t^{n}$ и $t^{n+1/2}$ снова делится пополам (промежуточная точка $t^{n+1/4}$).

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

Предиктор, первая подсхема (на первом полушаге $\Delta t/2$, учитывает только производную по $x$, неявная по $x$):

$$\frac{u_{j,k}^{n+1/4}-u_{j,k}^{n}}{\Delta t/2}=\sigma\frac{u_{j+1,k}^{n+1/4}-2u_{j,k}^{n+1/4}+u_{j-1,k}^{n+1/4}}{h_{x}^{2}} \tag{7.15}$$

Предиктор, вторая подсхема (на втором полушаге $\Delta t/2$, учитывает только производную по $y$, неявная по $y$):

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n+1/4}}{\Delta t/2}=\sigma\frac{u_{j,k+1}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j,k-1}^{n+1/2}}{h_{y}^{2}} \tag{7.16}$$

Результат последовательного решения подсхем (7.15), (7.16), называемых в совокупности предиктором, — значения $u(t,x,y)$ на шаге $(n+1/2)$. Для завершения расчёта на всём интервале $\Delta t$ используется поправочное соотношение — корректор:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}=\sigma\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+\sigma\frac{u_{j,k+1}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j,k-1}^{n+1/2}}{h_{y}^{2}}+f_{j,k}^{n+1/2} \tag{7.17}$$

Корректор берёт значения, посчитанные предиктором на шаге $n+1/2$, и по ним делает полный шаг от $n$ к $n+1$.

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

Каждая из подсхем предиктора (7.15), (7.16) абсолютно устойчива и решается методом прогонки (по $x$ и по $y$ соответственно). Коэффициенты трёхдиагонального вида (4.10) — обратите внимание, что из-за шага $\Delta t/2$ возникает множитель $\sigma/2$:

Первая подсхема предиктора (7.15), прогонка по $x$:

$$a_{j}=c_{j}=-\frac{\sigma}{2}\frac{\Delta t}{h_{x}^{2}},\quad b_{j}=1+\sigma\frac{\Delta t}{h_{x}^{2}},\quad \xi_{j,k}^{n}=u_{j,k}^{n}$$

Вторая подсхема предиктора (7.16), прогонка по $y$:

$$\tilde a_{k}=\tilde c_{k}=-\frac{\sigma}{2}\frac{\Delta t}{h_{y}^{2}},\quad \tilde b_{k}=1+\sigma\frac{\Delta t}{h_{y}^{2}},\quad \tilde\xi_{j,k}^{n+1/4}=u_{j,k}^{n+1/4}$$

Для обеих подсхем достаточное условие сходимости прогонки выполняется. Прогоночные соотношения и коэффициенты $\alpha_j,\beta_j$ / $\widetilde\alpha_k,\widetilde\beta_k$ имеют тот же вид (7.9)-(7.12), что и в схеме расщепления.

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

$$u_{j,k}^{n+1}=u_{j,k}^{n}+\sigma\Delta t\,\lambda_{xx}u_{j,k}^{n+1/2}+\sigma\Delta t\,\lambda_{yy}u_{j,k}^{n+1/2}+\Delta t\,f_{j,k}^{n+1/2} \tag{7.18}$$

См. 7.9.2. Характеристика подсхем. Метод решения.

4. Порядок аппроксимации и роли частей схемы

Правая часть корректора (7.17) аппроксимирована относительно точки $t^{n+1/2}$, поэтому разностный оператор по времени в левой части — центральная конечная разность второго порядка. Следовательно, схема предиктор-корректор имеет порядок

$$O\left(\Delta t^{2},\,h_{x}^{2},\,h_{y}^{2}\right)$$

что делает её более точной по сравнению со схемой расщепления. Распределение ролей:

• Предиктор (7.15), (7.16) — обеспечивает абсолютную устойчивость всей схемы;
• Корректор (7.17) — обеспечивает повышение порядка аппроксимации по времени (до второго).

5. Алгоритм решения

1. Задание начальных условий: цикл по $j=1,\ldots,N_x$ и $k=1,\ldots,N_y$, $u_{j,k}^{0}=\xi(x_j,y_k)$.
2. Цикл по времени $n$.
3. Предиктор, шаг по $x$ ($n\!\to\!n+1/4$): ГУ по $x$; прогонка (7.15) по $j$ в цикле по $k$ → $u_{j,k}^{n+1/4}$.
4. Предиктор, шаг по $y$ ($n+1/4\!\to\!n+1/2$): ГУ по $y$; прогонка (7.16) по $k$ в цикле по $j$ → $u_{j,k}^{n+1/2}$.
5. Корректор ($n\!\to\!n+1$): по рекуррентному соотношению (7.18) вычисляются все $u_{j,k}^{n+1}$ (граничные значения — из ГУ на шаге $n+1$).
6. Переход к следующему шагу по времени.

Блок-схема приведена на рис. 7.8 (см. 7.9.3. Алгоритм решения и страницу алгоритмов).

6. Разбор на семинаре

На семинаре 9 (Пример 3) схема разобрана на уравнении $\frac{\partial u}{\partial t}=3\frac{\partial^2 u}{\partial x^2}+4\frac{\partial^2 u}{\partial y^2}+e^{t}$. ① Предиктор по $x$: $\frac{u_{jk}^{n+1/4}-u_{jk}^{n}}{\Delta t/2}=3\frac{u_{j+1,k}^{n+1/4}-2u_{jk}^{n+1/4}+u_{j-1,k}^{n+1/4}}{h_x^2}$. ② Предиктор по $y$: $\frac{u_{jk}^{n+1/2}-u_{jk}^{n+1/4}}{\Delta t/2}=4\frac{u_{jk+1}^{n+1/2}-2u_{jk}^{n+1/2}+u_{jk-1}^{n+1/2}}{h_y^2}$. ③ Корректор: $\frac{u_{jk}^{n+1}-u_{jk}^{n+1/2}}{\Delta t}=3\Lambda_{xx}^{n+1/2}+4\Lambda_{yy}^{n+1/2}+e^{(n+1/2)\Delta t}$, рекуррентное соотношение $u_{jk}^{n+1}=u_{jk}^{n+1/2}+3\Delta t\,\Lambda_{xx}^{n+1/2}+4\Delta t\,\Lambda_{yy}^{n+1/2}+e^{(n+1/2)\Delta t}$, порядок $O(\Delta t^2, h_x^2, h_y^2)$. (На оси времени отметки $n$ — по $x$, $n+1/4$ — по $y$, $n+1/2$ — конец предиктора, $n+1$ — корректор.)

Итог

Схема состоит из трёх подсхем: предиктор делает «предварительный» расчёт на половину шага двумя прогонками (по $x$ на четверти $\Delta t/2$ и по $y$ на следующей четверти) и обеспечивает устойчивость; корректор по полученным значениям на $t^{n+1/2}$ явно (рекуррентно) пересчитывает решение на весь шаг $\Delta t$, что благодаря центрированию относительно $t^{n+1/2}$ даёт второй порядок по времени $O(\Delta t^2, h_x^2, h_y^2)$.

Вопрос 17. Схема расщепления (метод дробных шагов) для двумерных уравнений в частных производных 1-го порядка

1. Постановка задачи

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

$$\frac{\partial u}{\partial t}+v_{1}\frac{\partial u}{\partial x}+v_{2}\frac{\partial u}{\partial y}=f(t,x,y);\quad t\in[0,t_k],\ x\in[a,b],\ y\in[c,d] \tag{8.1}$$

Знаки коэффициентов $v_1, v_2$ определяют выбор конечной разности (правило выбора разности по направлению потока). Для случая $v_1>0, v_2>0$ (аппроксимация обеих производных по координатам левой разностью) неявная разностная схема имеет вид:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}+v_{1}\frac{u_{j,k}^{n+1}-u_{j-1,k}^{n+1}}{h_{x}}+v_{2}\frac{u_{j,k}^{n+1}-u_{j,k-1}^{n+1}}{h_{y}}=f_{j k}^{n} \tag{8.8}$$

(Аналогично записываются схемы 8.9-8.11 для других сочетаний знаков $v_1, v_2$.) Каждая из неявных схем (8.8)-(8.11) имеет первый порядок аппроксимации по времени и по каждой координате — $O(\Delta t, h_x, h_y)$ — и абсолютно устойчива (доказывается спектральным методом). Подробно: 8.2.1. Характеристика, 8.2.2. Исследование устойчивости.

2. Идея метода дробных шагов для уравнений 1-го порядка

Для решения неявных схем (8.8)-(8.11) применяется метод дробных шагов, подробно рассмотренный для параболических уравнений в разделе 7.6. Суть та же: интервал $\Delta t$ расщепляется пополам (точка $t^{n+1/2}$, рис. 8.5), что позволяет представить неявную схему в виде двух подсхем с более простым методом решения. Рассмотрим на примере схемы (8.8). Подробно: 8.2.3. Метод решения с использованием схемы расщепления.

3. Методика записи уравнений подсхем

Преобразуем неявную схему (8.8) в схему расщепления:

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}+v_{1}\frac{u_{j,k}^{n+1/2}-u_{j-1,k}^{n+1/2}}{h_{x}}=f_{j k}^{n}$$

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}+v_{2}\frac{u_{j,k}^{n+1}-u_{j,k-1}^{n+1}}{h_{y}}=0 \tag{8.12}$$

В первой подсхеме производная по времени аппроксимирована на первом полушаге $\Delta t$ — она неявная по координате $x$ (учитывает только $\partial u/\partial x$). Во второй подсхеме производная по времени аппроксимирована на втором полушаге $\Delta t$ — она неявная по координате $y$ (учитывает только $\partial u/\partial y$). Свободный член $f$ записан в первой подсхеме.

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

Складывая обе подсхемы, получаем соотношение, отличающееся от исходной схемы (8.8) только тем, что производная $\partial u/\partial x$ аппроксимирована не на $(n+1)$-м шаге, а на шаге $(n+1/2)$:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}+v_{1}\frac{u_{j,k}^{n+1/2}-u_{j-1,k}^{n+1/2}}{h_{x}}+v_{2}\frac{u_{j,k}^{n+1}-u_{j,k-1}^{n+1}}{h_{y}}=f_{j k}^{n}$$

Следовательно, схема расщепления (8.12) имеет тот же порядок аппроксимации, что и неявная схема (8.8) — первый по времени и по каждой координате:

$$O\left(\Delta t,\,h_{x},\,h_{y}\right)$$

5. Метод решения подсхем (рекуррентные соотношения)

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

$$u_{j,k}^{n+1/2}=\frac{u_{j,k}^{n}+v_{1}\dfrac{\Delta t}{h_{x}}u_{j-1,k}^{n+1/2}+\Delta t\,f_{j k}^{n}}{1+v_{1}\dfrac{\Delta t}{h_{x}}};\qquad u_{j,k}^{n+1}=\frac{u_{j,k}^{n+1/2}+v_{2}\dfrac{\Delta t}{h_{y}}u_{j,k-1}^{n+1}}{1+v_{2}\dfrac{\Delta t}{h_{y}}} \tag{8.13}$$

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

$$u_{1,k}^{n+1/2}=\varphi\left(t^{n+1/2},y_{k}\right),\qquad u_{j,1}^{n+1}=\psi\left(t^{n+1},x_{j}\right)$$

(При $v_1<0$ и/или $v_2<0$ используется правая конечная разность и задаётся правое граничное условие по соответствующей координате; вид рекуррентных соотношений меняется, а циклы в алгоритме идут в обратном порядке: $j=N_x-1,\ldots,1$ и/или $k=N_y-1,\ldots,1$. Порядок аппроксимации и устойчивость сохраняются.)

6. Алгоритм решения

1. Задание начальных условий $u_{j,k}^{0}=\xi(x_j,y_k)$.
2. Цикл по времени $n$.
3. Первый полушаг (по $x$): задать левое ГУ $u_{1,k}^{n+1/2}=\varphi(t^{n+1/2},y_k)$; в двойном цикле (по $k$, затем по $j=2,\ldots,N_x-1$ в порядке возрастания) по первому рекуррентному соотношению (8.13) вычислить $u_{j,k}^{n+1/2}$.
4. Второй полушаг (по $y$): задать ГУ $u_{j,1}^{n+1}=\psi(t^{n+1},x_j)$; в двойном цикле (по $j$, затем по $k=2,\ldots,N_y-1$) по второму рекуррентному соотношению (8.13) вычислить $u_{j,k}^{n+1}$.
5. Переход к следующему шагу по времени.

Блок-схема — рис. 8.6 (см. 8.2.4. Алгоритм решения с использованием схемы расщепления и страницу алгоритмов).

7. Разбор на семинаре

На семинаре 8 (Пример 2) расщепление разобрано на уравнении с конвективным членом $\frac{\partial u}{\partial t}+\frac{\partial u}{\partial y}=0{,}7\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+1$, где из-за наличия второй производной по $x$ первая подсхема решается уже прогонкой (тот же метод дробных шагов, но для уравнения, содержащего и первые, и вторые производные — см. также 8.4. Решение двумерных параболических уравнений с первыми производными).

Сравнение со схемами повышенного порядка

Для уравнений 1-го порядка существуют и более точные интерпретации неявной схемы (8.8): схема переменных направлений и схема предиктор-корректор, обе дают $O(\Delta t^2, h_x, h_y)$ (см. 8.2.5, 8.2.6, 8.3. Сравнительная характеристика). Схема расщепления (8.12) — самая простая: $O(\Delta t, h_x, h_y)$, абсолютно устойчива, решается двумя рекуррентными соотношениями (8.13).

Итог

Метод дробных шагов превращает неразрешимую напрямую неявную схему (8.8) для уравнения 1-го порядка в две последовательные подсхемы: первая неявна по $x$ (рекуррентный расчёт $n\!\to\!n+1/2$ по координате $x$), вторая неявна по $y$ (рекуррентный расчёт $n+1/2\!\to\!n+1$ по координате $y$). В отличие от параболического случая, где подсхемы трёхдиагональны и требуют прогонки, здесь каждая подсхема разрешается прямым рекуррентным соотношением (8.13). Схема абсолютно устойчива и имеет порядок $O(\Delta t, h_x, h_y)$.

Постановка задачи

Рассматривается двумерное (многомерное) дифференциальное уравнение в частных производных первого порядка в общем виде (учебник, гл. 8.1.1):

$$\frac{\partial u}{\partial t}+v_{1}\frac{\partial u}{\partial x}+v_{2}\frac{\partial u}{\partial y}=f(t,x,y);\quad t\in[0,t_k],\ x\in[a,b],\ y\in[c,d].$$

Параметры $v_1$ и $v_2$ (скорости переноса по координатам) могут быть как положительными, так и отрицательными.

Характеристика явной схемы

Производная по времени аппроксимируется правой конечной разностью (значение на $(n{+}1)$-м слое выражается через $n$-й), а первые производные по координатам — по правилу выбора конечной разности (раздел 5.8): при $v>0$ — левая разность, при $v<0$ — правая. Это даёт четыре явные схемы в зависимости от знаков $v_1,v_2$. Для случая $v_1>0,\ v_2>0$:

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n}}{\Delta t}+v_{1}\frac{u_{j,k}^{n}-u_{j-1,k}^{n}}{h_{x}}+v_{2}\frac{u_{j,k}^{n}-u_{j,k-1}^{n}}{h_{y}}=f_{jk}^{n}.$$

Остальные три схемы ($v_1>0,v_2<0$; $v_1<0,v_2>0$; $v_1<0,v_2<0$) отличаются заменой соответствующей разности на правую. Все характерные значения функции на $n$-м слое известны — в правой части неизвестных нет.

Каждая из схем имеет первый порядок аппроксимации и по времени, и по каждой координате: $O(\Delta t,h_x,h_y)$.

Для решения требуется задать начальное условие $u(t{=}0,x,y)=\xi(x,y)$ и граничные условия (1-го рода). При $v_1>0$ нужно левое ГУ по $x$: $u(t,x{=}a,y)=\varphi(t,y)$; при $v_2>0$ — нижнее ГУ по $y$: $u(t,x,y{=}c)=\psi(t,x)$ (для других знаков границы соответственно меняются — см. 8.1.1).

Исследование устойчивости (спектральный метод)

Свободный член $f$ на устойчивость не влияет, поэтому его отбрасывают и подставляют решение в виде гармоники (8.1.2):

$$u_{j,k}^{n}=\lambda^{n}e^{i\alpha j}e^{i\beta k},\quad \alpha\in[0,2\pi],\ \beta\in[0,2\pi].$$

После подстановки в схему ($v_1>0,v_2>0$) и деления на $\lambda^{n}e^{i\alpha j}e^{i\beta k}$ получаем множитель перехода:

$$\lambda=1-v_{1}\frac{\Delta t}{h_{x}}-v_{2}\frac{\Delta t}{h_{y}}+v_{1}\frac{\Delta t}{h_{x}}e^{-i\alpha}+v_{2}\frac{\Delta t}{h_{y}}e^{-i\beta}.$$

Вводя обозначения $r_{1}=v_{1}\dfrac{\Delta t}{h_{x}}>0$, $r_{2}=v_{2}\dfrac{\Delta t}{h_{y}}>0$, $r=r_1+r_2$, имеем $\lambda=1-r+r_{1}e^{-i\alpha}+r_{2}e^{-i\beta}$. Рассматривая простейший случай $\alpha=\beta$:

$$\lambda=1-r+re^{-i\alpha}.$$

Это окружность на комплексной плоскости с центром в точке $(1-r,\,0)$ и радиусом $|re^{i\alpha}|=r$. По необходимому условию устойчивости собственные числа оператора перехода должны лежать внутри или на границе единичного круга с центром в начале координат. Окружность лежит внутри единичного круга при $r<1$, выходит за него при $r>1$ и касается границы при $r=1$.

Условие устойчивости

Отсюда явная схема для случая $v_1>0,v_2>0$ устойчива при $r=r_1+r_2=v_{1}\dfrac{\Delta t}{h_{x}}+v_{2}\dfrac{\Delta t}{h_{y}}\le1$. Анализ остальных трёх схем по той же методике даёт обобщённое (по модулям скоростей) условие устойчивости для всех четырёх явных схем:

$$\boxed{\;|v_{1}|\frac{\Delta t}{h_{x}}+|v_{2}|\frac{\Delta t}{h_{y}}\le1\;}$$

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

Метод решения

Разностный шаблон содержит единственное неизвестное — значение функции на $(n{+}1)$-м слое. Выражая его, получаем рекуррентное соотношение (8.1.3):

$$u_{j,k}^{n+1}=u_{j,k}^{n}-v_{1}\frac{\Delta t}{h_{x}}\left(u_{j,k}^{n}-u_{j-1,k}^{n}\right)-v_{2}\frac{\Delta t}{h_{y}}\left(u_{j,k}^{n}-u_{j,k-1}^{n}\right)+\Delta t\,f_{jk}^{n}.$$

Оно позволяет рассчитать все внутренние значения на $(n{+}1)$-м слое по известным значениям на $n$-м слое; значения на границах берутся из ГУ: $u_{1,k}^{n+1}=\varphi(t^{n+1},y_k)$, $u_{j,1}^{n+1}=\psi(t^{n+1},x_j)$. Алгоритм (блок-схема) приведён в разделе «Блок-схемы».

Итог

• Явная схема для двумерного УрЧП 1-го порядка имеет порядок $O(\Delta t,h_x,h_y)$ и решается явным рекуррентным соотношением (без СЛАУ).
• Главный недостаток — условная устойчивость: $|v_{1}|\dfrac{\Delta t}{h_{x}}+|v_{2}|\dfrac{\Delta t}{h_{y}}\le1$, ограничивающая шаг по времени.

См. также семинар 8 и семинар 9 (двумерные схемы) и сравнительную характеристику схем (8.3).

Суть метода установления

Рассматривается обыкновенное дифференциальное уравнение 2-го порядка (учебник, гл. 10.1):

$$v\frac{du}{dx}=\sigma\frac{d^{2}u}{dx^{2}}-ku+f(x);\quad x\in[a,b],\ \sigma>0,\ v>0,$$

с граничными условиями 1-го рода $u(a)=\varphi_1$, $u(b)=\varphi_2$. При $k>0$ соответствующая разностная схема решается методом прогонки. Однако при $k=0$ (а тем более $k<0$) достаточное условие сходимости прогонки нарушается ($|a_j|+|c_j|=|b_j|$, нет строгого преобладания диагонали), и метод прогонки неприменим. Для этого случая используют метод установления (10.2).

Приведение стационарной задачи к нестационарной

В стационарное уравнение (при $k=0$) добавляют фиктивную производную по времени, превращая искомую функцию $u(x)$ в функцию двух переменных $\widetilde u(x,t)$:

$$v\frac{du}{dx}=\sigma\frac{d^{2}u}{dx^{2}}+f(x)\ \rightarrow\ \frac{\partial\widetilde u}{\partial t}+v\frac{\partial\widetilde u}{\partial x}=\sigma\frac{\partial^{2}\widetilde u}{\partial x^{2}}+f(x).$$

Полученное уравнение — одномерное параболическое, методы решения которого известны (гл. 4, 6).

Почему процесс сходится (условие установления)

Если начальное условие $\widetilde u(x,t{=}0)$ отличается от истинного решения стационарной задачи, производная по времени не равна нулю, и возникает динамика. Но поскольку граничные условия нестационарной задачи берутся из исходной стационарной задачи и не зависят от времени, при $t\to\infty$ решение «устанавливается»:

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

Процесс пошагового приближения к стационарному решению называют итерационным процессом, переход от $n$-го слоя к $(n{+}1)$-му — итерацией, $\Delta t$ — шагом итерации, $n$ — номером итерации.

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

Для нестационарного уравнения записывают неявную схему (с соблюдением правила выбора разности для $\partial u/\partial x$) — 10.4:

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

Эта схема (в отличие от явной) абсолютно устойчива, поэтому шаг итерации $\Delta t$ можно выбирать произвольно (грубым). Это ускоряет сходимость: число итераций $n\sim\dfrac{1}{h}$ против $n\sim\dfrac{1}{h^{2}}$ у явной схемы (метода простой итерации). Порядок аппроксимации — $O(\Delta t,h)$.

Метод решения — прогонка

Неявная схема приводится к трёхдиагональному виду; на каждой итерации решается методом прогонки. Прогоночные коэффициенты:

$$a_{j}=-\sigma\frac{\Delta t}{h^{2}},\quad b_{j}=1+v\frac{\Delta t}{h}+2\sigma\frac{\Delta t}{h^{2}},\quad c_{j}=-v\frac{\Delta t}{h}-\sigma\frac{\Delta t}{h^{2}},\quad \xi_{j}^{n}=u_{j}^{n}+\Delta t\,f(x_j).$$

Достаточное условие сходимости прогонки выполняется за счёт фиктивной производной (единица в $b_j$):

$$|a_{j}|+|c_{j}|=v\frac{\Delta t}{h}+2\sigma\frac{\Delta t}{h^{2}}<1+v\frac{\Delta t}{h}+2\sigma\frac{\Delta t}{h^{2}}=|b_{j}|.$$

Итерационным выражением служит прогоночное соотношение $u_j^{n}=\alpha_j u_{j+1}^{n}+\beta_j$ (как в общем случае, гл. 4.2.2).

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

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

$$u_{j}^{0}=f(x_j).$$

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

$$\left\|u^{n+1}-u^{n}\right\|=\sqrt{h\sum_{j=1}^{N}\left(u_{j}^{n+1}-u_{j}^{n}\right)^{2}}\le\varepsilon.$$

Важные правила метода установления

1. Исходное стационарное уравнение должно быть приведено к виду $v\,du/dx=\sigma\,d^2u/dx^2+\dots$, т.е. вторая производная — в правой части с положительным знаком, первая — в левой части.
2. Фиктивная производная по времени вводится в левую часть с положительным знаком (см. обобщение 10.6).
3. При $v<0$ для $\partial u/\partial x$ берётся правая конечная разность.

Метод установления с неявной схемой — основа лабораторной работы 4; блок-схемы — в разделе «Блок-схемы».

Суть метода установления

Метод применяется к ОДУ 2-го порядка $v\dfrac{du}{dx}=\sigma\dfrac{d^{2}u}{dx^{2}}-ku+f(x)$, $\sigma>0,v>0$ (гл. 10.1) в случае $k=0$, когда метод прогонки неприменим (нарушено достаточное условие сходимости). Стационарную задачу превращают в нестационарную, добавляя фиктивную производную по времени (10.2):

$$v\frac{du}{dx}=\sigma\frac{d^{2}u}{dx^{2}}+f(x)\ \rightarrow\ \frac{\partial\widetilde u}{\partial t}+v\frac{\partial\widetilde u}{\partial x}=\sigma\frac{\partial^{2}\widetilde u}{\partial x^{2}}+f(x),\qquad u(x)\to\widetilde u(x,t).$$

Так как граничные условия берутся из стационарной задачи и не зависят от времени, при $t\to\infty$ производная по времени стремится к нулю, а решение «устанавливается»: $\widetilde u(x,t)\to[u(x)]$, $\dfrac{\partial\widetilde u}{\partial t}\to0$. Пошаговое приближение к стационару — итерационный процесс; переход на $(n{+}1)$-й слой — итерация, $\Delta t$ — шаг итерации.

Схема Кранка-Николсона

Для нестационарного уравнения записывают схему Кранка-Николсона (полусумма явной и неявной аппроксимаций пространственных операторов) — 10.5:

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+\frac{v}{2}\frac{u_{j}^{n+1}-u_{j-1}^{n+1}}{h}+\frac{v}{2}\frac{u_{j}^{n}-u_{j-1}^{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).$$

Как и неявная схема, схема Кранка-Николсона абсолютно устойчива — шаг итерации можно выбирать произвольно. Её преимущество — второй порядок аппроксимации по времени: $O(\Delta t^{2},h)$. Это позволяет брать ещё более грубый шаг итерации, чем у неявной схемы, для той же точности, что дополнительно ускоряет сходимость:

$$n<\frac{1}{h}\quad(\text{против }n\sim\tfrac{1}{h}\text{ у неявной и }n\sim\tfrac{1}{h^{2}}\text{ у явной}).$$

Метод решения — прогонка

Схема приводится к трёхдиагональному виду и на каждой итерации решается методом прогонки (гл. 6.4). Прогоночные коэффициенты:

$$a_{j}=-\sigma\frac{\Delta t}{2h^{2}},\quad b_{j}=1+\frac{v}{2}\frac{\Delta t}{h}+\sigma\frac{\Delta t}{h^{2}},\quad c_{j}=-\frac{v}{2}\frac{\Delta t}{h}-\frac{\sigma}{2}\frac{\Delta t}{h^{2}},$$

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

Достаточное условие сходимости прогонки выполняется (единица в $b_j$ обеспечивает преобладание диагонали):

$$|a_{j}|+|c_{j}|=\frac{v}{2}\frac{\Delta t}{h}+\frac{\sigma\Delta t}{h^{2}}<1+\frac{v}{2}\frac{\Delta t}{h}+\frac{\sigma\Delta t}{h^{2}}=|b_{j}|.$$

Итерационное выражение — прогоночное соотношение $u_j^{n}=\alpha_j u_{j+1}^{n}+\beta_j$ (как в общем случае, 4.2.2).

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

Нулевая итерация (начальное условие из-за фиктивной производной) задаётся свободным членом: $u_{j}^{0}=f(x_j)$.

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

$$\left\|u^{n+1}-u^{n}\right\|=\sqrt{h\sum_{j=1}^{N}\left(u_{j}^{n+1}-u_{j}^{n}\right)^{2}}\le\varepsilon.$$

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

Схема	Порядок	Устойчивость	Решение	Число итераций
Явная (простая итерация)	$O(\Delta t,h)$	условно устойчива, $\Delta t\le\frac{h^2}{vh+2\sigma}$	рекуррентная формула	$n\sim 1/h^{2}$
Неявная	$O(\Delta t,h)$	абсолютно устойчива	прогонка	$n\sim 1/h$
Кранка-Николсона	$O(\Delta t^{2},h)$	абсолютно устойчива	прогонка	$n<1/h$

(см. обобщение, гл. 10.6). Алгоритмы — в разделе «Блок-схемы»; практика — лабораторная работа 4.

Постановка задачи

Дифференциальное уравнение эллиптического типа в общем виде (гл. 11.1):

$$\sigma\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)+f(x,y)=0;\quad x\in[a,b],\ y\in[c,d],\ \sigma>0,$$

с двумя ГУ 1-го рода по каждой координате: $u(a,y)=\varphi_1(y)$, $u(b,y)=\varphi_2(y)$, $u(x,c)=\psi_1(x)$, $u(x,d)=\psi_2(x)$. Примеры — стационарный реактор с продольным и поперечным перемешиванием, стационарное двумерное температурное поле $\lambda(\partial^2 T/\partial x^2+\partial^2 T/\partial y^2)+q=0$.

Суть метода установления

Прямая разностная аппроксимация даёт неразрешимую систему (одно уравнение связывает значения соседних узлов без явного выделения неизвестного). Поэтому, как и для ОДУ 2-го порядка, применяют метод установления — превращение стационарной задачи в нестационарную добавлением фиктивной производной по времени (11.2):

$$\sigma\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)+f(x,y)=0\ \rightarrow\ \frac{\partial\widetilde u}{\partial t}=\sigma\left(\frac{\partial^{2}\widetilde u}{\partial x^{2}}+\frac{\partial^{2}\widetilde u}{\partial y^{2}}\right)+f(x,y),\qquad u(x,y)\to\widetilde u(x,y,t).$$

Получается двумерное параболическое уравнение (гл. 7). Так как и граничные условия, и свободный член $f(x,y)$ не зависят от времени, при $t\to\infty$ решение «устанавливается»:

$$t\to\infty,\quad \widetilde u(x,y,t)\to u(x,y),\quad \frac{\partial\widetilde u}{\partial t}\to0.$$

Пошаговое приближение к стационару — итерационный процесс, $\Delta t$ — шаг итерации.

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

Переход с $n$-го слоя на $(n{+}1)$-й разбивается на два дробных шага через промежуточный слой $n+1/2$: на первом шаге неявно обрабатывается оператор по $x$, на втором — по $y$ (11.4):

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}=\sigma\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+f(x_j,y_k),$$

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

Каждая подсхема — аналог одномерной неявной параболической схемы, поэтому обе абсолютно устойчивы. Шаг итерации можно выбирать произвольно (грубым) — это ускоряет сходимость по сравнению с методом простой итерации: число итераций $n\sim\dfrac{1}{h}$ (против $n\sim\dfrac{1}{h^{2}}$ у явной схемы). Порядок аппроксимации $O(\Delta t,h_x^2,h_y^2)$.

Метод решения — прогонки по направлениям

Каждая из подсхем решается методом прогонки. Прогоночные коэффициенты:

1-я подсхема (прогонка по $x$): $a_j=c_j=-\sigma\dfrac{\Delta t}{h_x^2}$, $b_j=1+2\sigma\dfrac{\Delta t}{h_x^2}$, $\xi_{j,k}^{n}=u_{j,k}^{n}+\Delta t\,f(x_j,y_k)$.

2-я подсхема (прогонка по $y$): $\tilde a_k=\tilde c_k=-\sigma\dfrac{\Delta t}{h_y^2}$, $\tilde b_k=1+2\sigma\dfrac{\Delta t}{h_y^2}$, $\tilde\xi_{j,k}^{n+1/2}=u_{j,k}^{n+1/2}$.

Для обеих подсхем достаточное условие сходимости прогонки выполняется. Итерационные (прогоночные) выражения: для 1-й подсхемы $u_{j,k}^{n+1/2}=\alpha_j u_{j+1,k}^{n+1/2}+\beta_j$, для 2-й — $u_{j,k}^{n+1}=\widetilde\alpha_k u_{j,k+1}^{n+1}+\widetilde\beta_k$. Алгоритм аналогичен схеме расщепления для параболического 2D-уравнения (7.6.4).

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

Нулевая итерация (начальное условие из-за фиктивной производной) задаётся свободным членом: $u_{j,k}^{0}=f(x_j,y_k)$.

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

$$\left\|u^{n+1}-u^{n}\right\|=\sqrt{h^{2}\sum_{k=1}^{N}\sum_{j=1}^{N}\left(u_{j,k}^{n+1}-u_{j,k}^{n}\right)^{2}}\le\varepsilon.$$

Правила метода установления (11.7): вторые производные приводят к правой части с положительным знаком; фиктивная производная вводится в левую часть с положительным знаком. Алгоритмы — «Блок-схемы», практика — лабораторная работа 4.

Постановка задачи и суть метода установления

Эллиптическое уравнение (гл. 11.1):

$$\sigma\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)+f(x,y)=0;\quad \sigma>0,$$

с ГУ 1-го рода по обеим координатам. Прямая разностная схема неразрешима, поэтому применяют метод установления: стационарную задачу превращают в нестационарную добавлением фиктивной производной по времени (11.2):

$$\sigma\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)+f(x,y)=0\ \rightarrow\ \frac{\partial\widetilde u}{\partial t}=\sigma\left(\frac{\partial^{2}\widetilde u}{\partial x^{2}}+\frac{\partial^{2}\widetilde u}{\partial y^{2}}\right)+f(x,y),\qquad u(x,y)\to\widetilde u(x,y,t).$$

Получается двумерное параболическое уравнение. Поскольку и ГУ, и свободный член не зависят от времени, при $t\to\infty$ решение устанавливается: $\widetilde u(x,y,t)\to u(x,y)$, $\dfrac{\partial\widetilde u}{\partial t}\to0$. Пошаговое приближение к стационару — итерационный процесс, $\Delta t$ — шаг итерации.

Схема переменных направлений

Переход с $n$-го на $(n{+}1)$-й слой выполняется через промежуточный слой $n+1/2$ двумя подсхемами. Отличие от схемы расщепления: на каждом дробном шаге присутствуют оба пространственных оператора (с коэффициентом $\sigma/2$), но неявным попеременно становится то оператор по $x$, то по $y$ (отсюда «переменные направления») — 11.5:

$$\frac{u_{j,k}^{n+1/2}-u_{j,k}^{n}}{\Delta t}=\frac{\sigma}{2}\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+\frac{\sigma}{2}\frac{u_{j,k+1}^{n}-2u_{j,k}^{n}+u_{j,k-1}^{n}}{h_{y}^{2}}+f(x_j,y_k),$$

$$\frac{u_{j,k}^{n+1}-u_{j,k}^{n+1/2}}{\Delta t}=\frac{\sigma}{2}\frac{u_{j+1,k}^{n+1/2}-2u_{j,k}^{n+1/2}+u_{j-1,k}^{n+1/2}}{h_{x}^{2}}+\frac{\sigma}{2}\frac{u_{j,k+1}^{n+1}-2u_{j,k}^{n+1}+u_{j,k-1}^{n+1}}{h_{y}^{2}}.$$

Схема абсолютно устойчива (как и схема расщепления) — шаг итерации произвольный. Её преимущество — второй порядок аппроксимации по времени, что позволяет брать более грубый шаг итерации, чем у схемы расщепления, для той же точности и тем самым ускорить сходимость:

$$n<\frac{1}{h}.$$

Метод решения — прогонки по переменным направлениям

Каждая подсхема решается методом прогонки (обозначения $\lambda_{xx},\lambda_{yy}$ — разностные операторы вторых производных по $x$ и $y$):

1-я подсхема (прогонка по $x$): $a_j=c_j=-\dfrac{\sigma}{2}\dfrac{\Delta t}{h_x^2}$, $b_j=1+\sigma\dfrac{\Delta t}{h_x^2}$, $\xi_{j,k}^{n}=u_{j,k}^{n}+\dfrac{\sigma}{2}\Delta t\,\lambda_{yy}u_{j,k}^{n}+\Delta t\,f(x_j,y_k)$.

2-я подсхема (прогонка по $y$): $\tilde a_k=\tilde c_k=-\dfrac{\sigma}{2}\dfrac{\Delta t}{h_y^2}$, $\tilde b_k=1+\sigma\dfrac{\Delta t}{h_y^2}$, $\tilde\xi_{j,k}^{n+1/2}=u_{j,k}^{n+1/2}+\dfrac{\sigma}{2}\Delta t\,\lambda_{xx}u_{j,k}^{n+1/2}$.

Для обеих подсхем достаточное условие сходимости прогонки выполняется. Итерационные (прогоночные) выражения: для 1-й подсхемы $u_{j,k}^{n+1/2}=\alpha_j u_{j+1,k}^{n+1/2}+\beta_j$, для 2-й — $u_{j,k}^{n+1}=\widetilde\alpha_k u_{j,k+1}^{n+1}+\widetilde\beta_k$. Алгоритм аналогичен схеме расщепления для параболического 2D-уравнения (7.6.4), методики прогоночных коэффициентов — 4.2.2.

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

Нулевая итерация (начальное условие из-за фиктивной производной): $u_{j,k}^{0}=f(x_j,y_k)$.

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

$$\left\|u^{n+1}-u^{n}\right\|=\sqrt{h^{2}\sum_{k=1}^{N}\sum_{j=1}^{N}\left(u_{j,k}^{n+1}-u_{j,k}^{n}\right)^{2}}\le\varepsilon.$$

Сравнение схем метода установления для эллиптических уравнений

Схема	Порядок по времени	Устойчивость	Число итераций
Явная (простая итерация)	$O(\Delta t)$	условно устойчива, $\Delta t\le\frac{h^2}{4\sigma}$	$n\sim 1/h^{2}$
Расщепление (дробные шаги)	$O(\Delta t)$	абсолютно устойчива	$n\sim 1/h$
Переменные направления	$O(\Delta t^{2})$	абсолютно устойчива	$n<1/h$

(см. обобщение, гл. 11.7). Все подсхемы решаются прогонкой; нулевая итерация — свободный член. Блок-схемы — «Блок-схемы», практика — лабораторная работа 4.

1. Сложные системы уравнений и основной подход к их решению

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

Методы численного решения уравнений (явные и неявные схемы) применимы и к системам. Однако чтобы решить какое-либо уравнение системы, нужно знать значения всех входящих в него функций, определяемых из других уравнений. Поэтому при записи разностных схем для систем применяют принцип замороженных коэффициентов: все «чужие» функции берут с предыдущего расчётного шага $n$ (на первом шаге — из начальных условий). Это делает разностную схему разрешимой.

Иллюстрация на простой системе двух одномерных параболических уравнений (учебник, гл. 15.1):

$$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}-k u v ; \qquad \frac{\partial v}{\partial t}=\frac{\partial^{2} v}{\partial x^{2}}-k u v .$$

Свободный член первого уравнения содержит $v$, а второго — $u$. Неявная схема с замороженными коэффициентами:

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}=\frac{u_{j+1}^{n+1}-2 u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}-k\,u_{j}^{n+1} v_{j}^{n};$$ $$\frac{v_{j}^{n+1}-v_{j}^{n}}{\Delta t}=\frac{v_{j+1}^{n+1}-2 v_{j}^{n+1}+v_{j-1}^{n+1}}{h^{2}}-k\,u_{j}^{n} v_{j}^{n+1}.$$

«Чужая» функция ($v_j^n$ в первом, $u_j^n$ во втором) взята с $n$-го шага — схема разрешима.

2. Математическая модель процесса массовой кристаллизации

Рассмотрим ёмкостной кристаллизатор периодического действия, в котором массовая кристаллизация идёт за счёт охлаждения раствора. Аппарат идеального смешения, поэтому все градиенты (концентрации, температуры) отсутствуют. Модель базируется на законах сохранения массы и энергии и имеет вид (учебник, гл. 15.2):

— уравнение изменения концентрации раствора (ОДУ + интеграл, материальный баланс):

$$\frac{d c}{d t}=-\int_{0}^{R} \rho_{2}^{0} f \eta\, d r \tag{15.2}$$

— уравнение изменения температуры в реакторе (ОДУ + интеграл, тепловой баланс):

$$\left[\rho_{1} C_{1 T}+\rho_{2}^{0} C_{2 T} \int_{0}^{R} f r\, d r\right] \frac{d T}{d t}=\Delta H \int_{0}^{R} \rho_{2}^{0} f \eta\, d r-K F\left(T-T_{x}\right) \tag{15.3}$$

— уравнение баланса числа частиц (дифференциальное уравнение в частных производных 1-го порядка по координате-размеру $r$):

$$\frac{\partial f}{\partial t}+\frac{\partial(f \eta)}{\partial r}=0 \tag{15.4}$$

— выражение для скорости роста кристалла (алгебраическое):

$$\eta=\frac{d r}{d t}=k_{2} S_{r}\left(c-c_{S}\right)^{m} \tag{15.5}$$

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

$$c(t=0)=c_{0}, \quad T(t=0)=T_{0}, \quad f(t=0, r)=0 \tag{15.6}$$

— граничное условие к уравнению (15.4):

$$f\left(t, r=r_{0}\right) \cdot \eta\left(t, r=r_{0}\right)=I=k_{1}\left(c-c_{S}\right)^{p} \tag{15.7}$$

— равновесная концентрация (алгебраическое, линеаризация в узком диапазоне температур):

$$c_{S}=a+b T \tag{15.8}$$

Здесь $c$ — концентрация кристаллизующегося компонента; $\rho_1$, $\rho_2^0$ — плотности раствора и кристалла; $T$ — температура; $C_{1T},C_{2T}$ — теплоёмкости раствора и кристалла; $\eta$ — скорость роста кристаллов; $f(r)\,dr$ — число кристаллов в единице объёма с размером от $r$ до $r+dr$; $r_0$ — размер зародыша; $I$ — скорость зародышеобразования; $\Delta H$ — тепловой эффект; $K$ — коэффициент теплопередачи; $F$ — поверхность кристаллизатора; $T_x$ — температура хладагента; $c_S$ — равновесная концентрация; $(c-c_S)$ — пересыщение; $S_r$ — поверхность кристалла размером $r$; $k_1,k_2$ — кинетические константы; $m,p$ — показатели степени при пересыщении.

Состав системы (15.2)–(15.8): два интегро-дифференциальных уравнения (15.2), (15.3), одно дифференциальное уравнение в частных производных 1-го порядка (15.4), два алгебраических уравнения (15.5), (15.8), а также начальные (15.6) и граничные (15.7) условия. Это и есть искомая система, состоящая из уравнений в частных производных и обыкновенных (а также интегро-дифференциальных) дифференциальных уравнений.

3. Метод численного решения системы

Интеграл в правых частях (15.2), (15.3) аппроксимируется методом прямоугольников (суммой по разностной сетке размеров $r_j$, $j=1\dots N_r$ с шагом $\Delta r$), см. учебник, гл. 14.3.1. Поскольку искомой функцией для (15.2) является $c$, а для (15.3) — $T$, то «чужие» функции под интегралом ($f$, $\eta$) берутся с $n$-го шага (замороженные коэффициенты).

Разностные схемы для (15.2) и (15.3):

$$\frac{c^{n+1}-c^{n}}{\Delta t}=-\sum_{j=1}^{N_{r}} \rho_{2}^{0} f_{j}^{n} \eta_{j}^{n} \Delta r,$$ $$\left[\rho_{1}^{n} C_{1 T}^{n}+\rho_{2}^{0} C_{2 T} \sum_{j=1}^{N_{r}} f_{j}^{n} r_{j} \Delta r\right] \frac{T^{n+1}-T^{n}}{\Delta t}=\Delta H \sum_{j=1}^{N_{r}} \rho_{2}^{0} f_{j}^{n} \eta_{j}^{n} \Delta r+K F\left(T^{n+1}-T_{x}\right) \tag{15.9}$$

Эти схемы разрешаются явно — рекуррентными соотношениями:

$$c^{n+1}=c^{n}-\Delta t \sum_{j=1}^{N_{r}} \rho_{2}^{0} f_{j}^{n} \eta_{j}^{n} \Delta r,$$ $$T^{n+1}=\frac{\left[\rho_{1}^{n} C_{1 T}^{n}+\rho_{2}^{0} C_{2 T} \sum_{j} f_{j}^{n} r_{j} \Delta r\right] T^{n}+\Delta t\left(\Delta H \sum_{j} \rho_{2}^{0} f_{j}^{n} \eta_{j}^{n} \Delta r-K F T_{x}\right)}{\rho_{1}^{n} C_{1 T}^{n}+\rho_{2}^{0} C_{2 T} \sum_{j} f_{j}^{n} r_{j} \Delta r-\Delta t\, K F}.$$

Уравнение в частных производных 1-го порядка (15.4) аппроксимируется с учётом правила выбора конечной разности по $r$ (для растущих кристаллов $\eta>0$ — левая разность) и принципа замороженных коэффициентов для $\eta$ (берётся с $n$-го шага):

$$\frac{f_{j}^{n+1}-f_{j}^{n}}{\Delta t}+\frac{f_{j}^{n+1} \eta_{j}^{n}-f_{j-1}^{n+1} \eta_{j-1}^{n}}{\Delta r}=0 \tag{15.10}$$

откуда рекуррентное соотношение для плотности распределения:

$$f_{j}^{n+1}=\frac{f_{j}^{n}+\dfrac{\Delta t}{\Delta r} f_{j-1}^{n+1} \eta_{j-1}^{n}}{1+\dfrac{\Delta t}{\Delta r} \eta_{j}^{n}}.$$

Алгебраические уравнения и условия в разностном виде:

$$c^{0}=c_{0}, \quad T^{0}=T_{0}, \quad f_{j}^{0}=0;$$ $$c_{S}^{n+1}=a+b T^{n+1}; \quad \eta_{j}^{n+1}=k_{2} S_{j}\left(c^{n+1}-c_{S}^{n+1}\right)^{m};$$ $$I^{n+1}=k_{1}\left(c^{n+1}-c_{S}^{n+1}\right)^{p}; \quad f_{1}^{n+1}=I^{n+1} / \eta_{1}^{n+1}.$$

Порядок вычислений на каждом шаге по времени

1. По известным с $n$-го шага значениям $f_j^n$, $\eta_j^n$ вычисляем суммы (интегралы) и по явным рекуррентным формулам находим $c^{n+1}$ и $T^{n+1}$.
2. По $T^{n+1}$ из (15.8) находим равновесную концентрацию $c_S^{n+1}$.
3. По $c^{n+1}$ и $c_S^{n+1}$ из (15.5) находим скорости роста $\eta_j^{n+1}$, а из (15.7) — скорость зародышеобразования $I^{n+1}$ и граничное значение $f_1^{n+1}$.
4. По рекуррентному соотношению (15.10) последовательно (по $j$) рассчитываем всю функцию распределения $f_j^{n+1}$.
5. Переходим к следующему шагу по времени.

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

Источники: гл. 1.8 (модель кристаллизации — пример системы интегро-дифференциальных уравнений), гл. 14.3.1 (решение интегро-дифференциальных уравнений модели), гл. 15.1 (принцип замороженных коэффициентов), гл. 15.2.1 и 15.2.2 (уравнения модели и рекуррентные соотношения).

1. Зачем нужно обезразмеривание

Значения, получаемые численными методами, отличаются от истинных из-за ошибки аппроксимации. Если уравнения модели содержат переменные, значения которых отличаются по порядкам, алгоритм может оказаться непригодным: погрешности при определении величин больших порядков, не значимые для них самих, будут сильно искажать значения величин меньших порядков. Например, в модели кристаллизации функция $f$ имеет порядок $\sim 10^{20}$, а скорость роста $\eta\sim 10^{-10}$. Поэтому перед построением алгоритма уравнения приводят к безразмерному виду (обезразмеривают переменные), чтобы все переменные модели имели одинаковый порядок (учебник, гл. 2.1.1).

2. Общий алгоритм приведения к безразмерному виду

Покажем алгоритм на уравнении диффузии с реакцией (семинар 1):

$$\frac{\partial C}{\partial t}=D \frac{\partial^{2} C}{\partial x^{2}}-k C \tag{1}$$

Шаг 1. Вводим характерные параметры процесса (масштабы) — постоянные с индексом 0, имеющие размерность соответствующих величин:

$$C_{0},\ t_{0},\ x_{0},\ D_{0},\ k_{0}=\mathrm{const}.$$

Часть масштабов обычно известна из постановки задачи (например, характерное время $t_0$, характерный размер $x_0$, характерная концентрация $C_0$), а часть подлежит определению.

Шаг 2. Вводим безразмерные переменные как отношение величины к её масштабу:

$$C'=\frac{C}{C_{0}}, \quad t'=\frac{t}{t_{0}}, \quad x'=\frac{x}{x_{0}}, \quad D'=\frac{D}{D_{0}}, \quad k'=\frac{k}{k_{0}}.$$

Отсюда выражаем размерные переменные через масштабы и безразмерные:

$$C=C' C_{0}, \quad t=t' t_{0}, \quad x=x' x_{0}, \quad D=D' D_{0}, \quad k=k' k_{0}.$$

Шаг 3. Подставляем в исходное уравнение и выносим масштабы за знаки производных:

$$\frac{C_{0}}{t_{0}} \frac{\partial C'}{\partial t'}=\frac{D_{0} C_{0}}{x_{0}^{2}} D' \frac{\partial^{2} C'}{\partial x'^{2}}-k_{0} C_{0}\, k' C' \tag{2}$$

Шаг 4. Делим уравнение на коэффициент при производной по времени $\left(\dfrac{C_0}{t_0}\right)$, чтобы при $\partial C'/\partial t'$ стояла единица. При остальных слагаемых образуются безразмерные комплексы характерных параметров:

$$\frac{\partial C'}{\partial t'}=\left[\frac{D_{0} t_{0}}{x_{0}^{2}}\right] D' \frac{\partial^{2} C'}{\partial x'^{2}}-\left[k_{0} t_{0}\right] k' C' \tag{3}$$

Шаг 5. Находим неизвестные масштабы из условия совпадения с исходным уравнением. Чтобы безразмерное уравнение совпало по форме с исходным, безразмерные комплексы приравнивают единице:

$$\left[\frac{D_{0} t_{0}}{x_{0}^{2}}\right]=1 \;\Rightarrow\; D_{0}=\frac{x_{0}^{2}}{t_{0}}\ \left[\tfrac{\text{м}^2}{\text{с}}\right], \qquad \left[k_{0} t_{0}\right]=1 \;\Rightarrow\; k_{0}=\frac{1}{t_{0}}\ \left[\tfrac{1}{\text{с}}\right].$$

Шаг 6. Итоговое безразмерное уравнение по форме совпадает с исходным, но все переменные в нём имеют одинаковый порядок:

$$\frac{\partial C'}{\partial t'}=D' \frac{\partial^{2} C'}{\partial x'^{2}}-k' C'. \tag{4}$$

3. Применение к системе уравнений кристаллизации

Для системы алгоритм тот же, но выполняется для каждого уравнения, причём масштабы общие для всей системы. Рассмотрим модель кристаллизации в ёмкостном реакторе идеального смешения (учебник, гл. 2.1.2, семинар 1, пример 5):

$$\frac{d c}{d t}=-\int_{0}^{R} \rho_{2}^{0} f \eta\, d r, \qquad \frac{\partial f}{\partial t}+\eta \frac{\partial f}{\partial r}=0. \tag{2.1}$$

Безразмерные переменные:

$$c'=\frac{c}{c_{0}}, \; t'=\frac{t}{t_{0}}, \; r'=\frac{r}{r_{0}}, \; \eta'=\frac{\eta}{\eta_{0}}, \; f'=\frac{f}{f_{0}}, \; \rho_{2}^{\prime 0}=\frac{\rho_{2}^{0}}{\rho_{0}}.$$

Известны характерные время $t_0$, размер кристалла $r_0$ и концентрация $c_0$. Поскольку плотность и концентрация имеют одинаковую размерность, принимают $\rho_0=c_0$. А характерную скорость роста $\eta_0$ и характерную плотность функции распределения $f_0$ измерить непосредственно нельзя — их находят из безразмерных комплексов.

Выражаем размерные переменные через масштабы ($\rho_2^0=\rho_2^{\prime 0}\cdot c_0$) и подставляем в (2.1):

$$\frac{c_{0}}{t_{0}} \frac{d c'}{d t'}=-c_{0} f_{0} \eta_{0} r_{0} \int_{0}^{R} \rho_{2}^{\prime 0} f' \eta'\, d r', \qquad \frac{f_{0}}{t_{0}} \frac{\partial f'}{\partial t'}+\frac{f_{0} \eta_{0}}{r_{0}} \eta' \frac{\partial f'}{\partial r'}=0. \tag{2.2}$$

Второе уравнение делим на $f_0/t_0$:

$$\frac{\partial f'}{\partial t'}+\frac{t_{0} \eta_{0}}{r_{0}} \eta' \frac{\partial f'}{\partial r'}=0.$$

Приравниваем комплекс перед вторым слагаемым к единице — это даёт характерную скорость роста (типичный пример комплекса вида $v_0 t_0 / x_0 = 1$):

$$\frac{t_{0} \eta_{0}}{r_{0}}=1 \;\Rightarrow\; \eta_{0}=\frac{r_{0}}{t_{0}}. \tag{2.3}$$

Первое уравнение делим на $c_0/t_0$:

$$\frac{d c'}{d t'}=-t_{0} f_{0} \eta_{0} r_{0} \int_{0}^{R} \rho_{2}^{\prime 0} f' \eta'\, d r'.$$

Приравниваем комплекс перед интегралом к единице и с учётом (2.3) находим характерную плотность распределения:

$$t_{0} f_{0} \eta_{0} r_{0}=1 \;\Rightarrow\; f_{0}=\frac{1}{t_{0} r_{0} \eta_{0}}=\frac{1}{t_{0} r_{0} \cdot \frac{r_{0}}{t_{0}}}=\frac{1}{r_{0}^{2}}. \tag{2.4}$$

Вывод: если масштабы $\eta_0$, $f_0$ выбраны по формулам (2.3), (2.4), то безразмерная система (2.2) по форме полностью совпадает с исходной (2.1), но теперь все переменные имеют одинаковый порядок, и расчётные ошибки не приводят к искажению функций малого порядка.

4. Замечания и критерии подобия

• Степень свободы. Если число неизвестных масштабов больше числа комплексов (как в кристаллизации: три неизвестных $f_0,\eta_0,r_0$ при двух комплексах), один масштаб (например $r_0$) задают произвольно, остальные выражают через него (семинар 1).
• Уравнение без времени. Если в уравнении нет переменной времени (например реактор идеального вытеснения $v\,dC/dx=-kC$), вводят $t_0$ и умножают уравнение на $t_0$, чтобы комплексы имели физический смысл; тогда $v_0 t_0/x_0=1\Rightarrow v_0=x_0/t_0$ (семинар 0, семинар 1, пример 4).
• Критерии подобия. Если масштабов меньше, чем комплексов, оставшиеся комплексы не приравнивают единице, а оставляют в виде безразмерных критериев подобия $\Pi_1,\beta,\dots$ (семинар 1, пример с движением тела в атмосфере).

Источники: гл. 2.1.1 (необходимость безразмерных переменных), гл. 2.1.2 (методика определения неизвестных характерных параметров), семинар 0 и семинар 1 (алгоритм и 5 разобранных примеров).

1. Когда применяют метод тестовых задач

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

2. Идея метода — два опорных факта

Факт 1. Устойчивость разностной схемы не зависит от вида свободного члена дифференциального уравнения, если он содержит только независимые переменные и не содержит саму искомую функцию. Иными словами, разностные схемы, отличающиеся только свободным членом вида

$$f\left(t^{n}, x_{j}, y_{k}, \ldots\right) \tag{15.11}$$

обладают одинаковым типом устойчивости.

Факт 2. Если известно истинное решение дифференциального уравнения, то его можно сравнить с численным решением, полученным по разностной схеме, и определить, устойчива она или нет:

$$\left\|u_{j, k, \ldots}^{n}-\left[u\left(t^{n}, x_{j}, y_{k}, \ldots\right)\right]\right\| \leq \varepsilon \;\Rightarrow\; \text{схема устойчива};$$ $$\left\|u_{j, k, \ldots}^{n}-\left[u\left(t^{n}, x_{j}, y_{k}, \ldots\right)\right]\right\| > \varepsilon \;\Rightarrow\; \text{схема неустойчива}. \tag{15.12}$$

Здесь $[u(t^n,x_j,\dots)]$ — истинное решение в узле сетки, $u_{j,\dots}^n$ — решение разностной схемы в том же узле, $\varepsilon$ — заданная допустимая погрешность.

3. Суть метода (последовательность действий)

1. Задают тест — некоторую функцию от независимых переменных $\tilde u$ (называемую тестом).
2. Строят тестовую задачу — дифференциальное уравнение, для которого выбранный тест является точным истинным решением; при этом новое уравнение должно отличаться от исходного только свободным членом вида (15.11). (Тогда по Факту 1 тип устойчивости схемы для тестовой задачи и для исходного уравнения один и тот же.)
3. Решают тестовую задачу исследуемой разностной схемой и сравнивают численный результат с истинным решением (тестом) в узлах сетки по критерию (15.12).
4. Делают вывод: если решение тестовой задачи подтверждает устойчивость, то эту схему можно использовать и для исходного уравнения, точное решение которого неизвестно; если выявлена неустойчивость — выбирают другую схему.

Правило выбора теста: целесообразно брать функцию того же типа, что и свободный член (15.11): если он алгебраический — тест берут алгебраическим; если тригонометрический — тригонометрическим; если экспоненциальный — экспонентой, и т. д.

4. Пример построения тестовой задачи

(учебник, гл. 15.3.2) Дано одномерное параболическое уравнение, для которого нужно подобрать устойчивую схему:

$$\frac{\partial u}{\partial t}=10 \frac{\partial^{2} u}{\partial x^{2}}+x^{2}, \quad x \in[0,1], \; t \in[0,1]. \tag{15.13}$$

Свободный член $x^2$ — алгебраический, поэтому тест задаём в виде алгебраической функции:

$$\tilde{u}=t x^{2}. \tag{15.14}$$

Шаг 1. Находим производные теста, входящие в уравнение:

$$\frac{\partial \tilde{u}}{\partial t}=x^{2}, \quad \frac{\partial \tilde{u}}{\partial x}=2 t x, \quad \frac{\partial^{2} \tilde{u}}{\partial x^{2}}=2 t. \tag{15.15}$$

Шаг 2. Записываем уравнение, истинным решением которого должен быть тест, добавляя в свободный член неизвестную функцию $\varphi(t,x)$:

$$\frac{\partial \tilde{u}}{\partial t}=10 \frac{\partial^{2} \tilde{u}}{\partial x^{2}}+x^{2}+\varphi(t, x).$$

Шаг 3. Подставляем производные (15.15) и определяем $\varphi$:

$$x^{2}=10\cdot 2t+x^{2}+\varphi(t, x) \;\Rightarrow\; \varphi(t, x)=-20 t.$$

Шаг 4. Получаем тестовую задачу, отличающуюся от исходной (15.13) только свободным членом, и для которой тест (15.14) является точным решением:

$$\frac{\partial \tilde{u}}{\partial t}=10 \frac{\partial^{2} \tilde{u}}{\partial x^{2}}+x^{2}-20 t. \tag{15.16}$$

Шаг 5. Начальное и граничные условия задаём с помощью самого теста (15.14), чтобы получить замкнутую корректную задачу:

$$\tilde{u}(t=0, x)=0; \quad \tilde{u}(t, x=0)=0; \quad \tilde{u}(t, x=1)=t.$$

Теперь задачу (15.16) решают исследуемой разностной схемой и сравнивают численное решение с точным $\tilde u = t x^2$ по критерию (15.12). Так как тестовая задача отличается от исходной только свободным членом (15.11), вывод об устойчивости автоматически переносится на исходное уравнение (15.13).

5. Применение к модели кристаллизации

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

Источники: гл. 15.3.1 (метод тестовых задач), гл. 15.3.2 (пример построения тестовой задачи), гл. 15.3.3 (устойчивость схем модели кристаллизации). О независимости устойчивости от свободного члена см. также гл. 3.6.

1. Постановка задачи

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

$$\frac{\partial c}{\partial t}+v \frac{\partial c}{\partial x}=D_{L} \frac{\partial^{2} c}{\partial x^{2}}-k c.$$

В общем виде рассматривается уравнение:

$$\frac{\partial u}{\partial t}+v \frac{\partial u}{\partial x}=\sigma \frac{\partial^{2} u}{\partial x^{2}}+f(t, x); \quad v>0,\ \sigma>0. \tag{6.1}$$

с начальным и граничными условиями (для определённости — 1-го рода):

$$u(t=0, x)=\xi(x), \quad u(t, x=0)=\varphi_{1}(t), \quad u(t, x=l)=\varphi_{2}(t).$$

Правило выбора конечной разности для конвективного члена $v\,\partial u/\partial x$ (по направлению потока): при $v>0$ используют левую конечную разность, при $v<0$ — правую. Ниже рассмотрен случай $v>0$ (левая разность); случай $v<0$ аналогичен. Отдельно рассматривается центральная разность.

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

(учебник, гл. 6.2.1) Конвективный член — левой разностью, всё на $n$-м шаге:

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+v \frac{u_{j}^{n}-u_{j-1}^{n}}{h}=\sigma \frac{u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}}{h^{2}}+f\left(t^{n}, x_{j}\right). \tag{6.2}$$

• Порядок аппроксимации $O(\Delta t, h)$ — первый и по времени, и по координате (левая разность даёт $O(h)$, что грубее, чем $O(h^2)$ от второй производной).
• Устойчивость условная (спектральный метод: собственные числа лежат на окружности с центром $(1-q,0)$ радиусом $r=v\frac{\Delta t}{h}$):

$$v \frac{\Delta t}{h}+2 \sigma \frac{\Delta t}{h^{2}} \leq 1. \tag{6.3}$$

Метод решения — явный, по рекуррентному соотношению (одна неизвестная на $(n+1)$-м шаге):

$$u_{j}^{n+1}=u_{j}^{n}+v \frac{\Delta t}{h}\left(u_{j-1}^{n}-u_{j}^{n}\right)+\sigma \frac{\Delta t}{h^{2}}\left(u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}\right)+\Delta t\, f\left(t^{n}, x_{j}\right). \tag{6.4}$$

Граничные точки $u_1^{n+1}, u_N^{n+1}$ определяются из граничных условий. Алгоритм аналогичен явной схеме для обычного параболического уравнения (гл. 6.2.2).

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

(учебник, гл. 6.3.1) Конвективный член и вторую производную берут на $(n+1)$-м шаге:

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+v \frac{u_{j}^{n+1}-u_{j-1}^{n+1}}{h}=\sigma \frac{u_{j+1}^{n+1}-2 u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}+f\left(t^{n}, x_{j}\right). \tag{6.5}$$

• Порядок аппроксимации $O(\Delta t, h)$.
• Абсолютно устойчива (спектральный метод: $1/\lambda=1+q-r e^{-i\alpha}$ лежит на окружности с центром $(1+q,0)$ радиусом $r$; так как $1+q>r$, при любых шагах $|1/\lambda|\ge 1\Rightarrow|\lambda|\le 1$).

Метод решения — метод прогонки (три неизвестных на $(n+1)$-м шаге). Приводим к трёхдиагональному виду:

$$-\sigma \frac{\Delta t}{h^{2}} u_{j+1}^{n+1}+\left(1+v \frac{\Delta t}{h}+2 \sigma \frac{\Delta t}{h^{2}}\right) u_{j}^{n+1}-\left(v \frac{\Delta t}{h}+\sigma \frac{\Delta t}{h^{2}}\right) u_{j-1}^{n+1}=u_{j}^{n}+\Delta t\, f\left(t^{n}, x_{j}\right).$$

Коэффициенты прогонки: $a_j=-\sigma\frac{\Delta t}{h^2}$, $b_j=1+v\frac{\Delta t}{h}+2\sigma\frac{\Delta t}{h^2}$, $c_j=-v\frac{\Delta t}{h}-\sigma\frac{\Delta t}{h^2}$. Достаточное условие сходимости прогонки $|a_j|+|c_j|<|b_j|$ выполняется автоматически (гл. 6.3.2).

4. Разностная схема Кранка-Николсона

(учебник, гл. 6.4) Чтобы сохранить второй порядок по времени, и конвективный член, и вторую производную представляют как полусумму слагаемых на $(n+1)$-м и $n$-м шагах:

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+\frac{v}{2} \frac{u_{j}^{n+1}-u_{j-1}^{n+1}}{h}+\frac{v}{2} \frac{u_{j}^{n}-u_{j-1}^{n}}{h}=\frac{\sigma}{2} \frac{u_{j+1}^{n+1}-2 u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}+\frac{\sigma}{2} \frac{u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}}{h^{2}}+f\left(t^{n}, x_{j}\right). \tag{6.7}$$

• Порядок аппроксимации $O(\Delta t^{2}, h)$ — повышенный по времени.
• Абсолютно устойчива.
• Решается методом прогонки (три неизвестных). Коэффициенты: $a_j=-\sigma\frac{\Delta t}{2h^2}$, $b_j=1+v\frac{\Delta t}{2h}+\sigma\frac{\Delta t}{h^2}$, $c_j=-v\frac{\Delta t}{2h}-\sigma\frac{\Delta t}{2h^2}$; достаточное условие сходимости выполняется.

5. Неявная схема с центральной разностью конвективного члена

(учебник, гл. 6.5.1) Производную $\partial u/\partial x$ аппроксимируют центральной разностью:

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+v \frac{u_{j+1}^{n+1}-u_{j-1}^{n+1}}{2 h}=\sigma \frac{u_{j+1}^{n+1}-2 u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}+f\left(t^{n}, x_{j}\right). \tag{6.8}$$

• Порядок аппроксимации $O(\Delta t, h^{2})$ — повышенный по координате (центральная разность точнее).
• Абсолютно устойчива при любом знаке $v$ ($|1/\lambda|=\sqrt{(1+4\sigma\frac{\Delta t}{h^2}\sin^2\frac{\alpha}{2})^2+(v\frac{\Delta t}{h}\sin\alpha)^2}\ge 1$).
• Решается прогонкой, но достаточное условие сходимости накладывает ограничение: либо на шаг по координате $h\le \dfrac{2\sigma}{|v|}$ (тогда ограничений на $\Delta t$ нет), либо при бо́льших $h$ — условие $|v|\frac{\Delta t}{h}<1+2\sigma\frac{\Delta t}{h^2}$ (гл. 6.5.2).

6. Сравнительная характеристика и выбор схемы

Схема	Аппроксимация конв. члена	Порядок	Устойчивость	Решение
Явная	левая разность (на $n$)	$O(\Delta t, h)$	условная: $v\frac{\Delta t}{h}+2\sigma\frac{\Delta t}{h^2}\le 1$	рекуррентно (6.4)
Неявная	левая разность (на $n{+}1$)	$O(\Delta t, h)$	абсолютная	прогонка
Кранка-Николсона	полусумма (левая)	$O(\Delta t^2, h)$	абсолютная	прогонка
Неявная (центр. разность)	центральная (на $n{+}1$)	$O(\Delta t, h^2)$	абсолютная (любой знак $v$)	прогонка (условие на $h$/$\Delta t$)

При $v<0$ для конвективного члена берут правую разность. Если $v$ знакопеременна или знак заранее неизвестен — применяют неявную схему с центральной разностью. Практические примеры записи всех трёх схем (явной, неявной с прогонкой и Кранка-Николсона) с аппроксимацией граничных условий разобраны на семинаре 7.

Источники: гл. 6.1 (постановка), 6.2 (явная), 6.3 (неявная), 6.4 (Кранка-Николсона), 6.5 (центральная разность), 6.6 (сравнение), семинар 7.