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

Что требуется

Определить тип уравнения (эллиптическое / параболическое / гиперболическое) и обосновать выбор. Уравнение:

$$0{,}2\,\frac{\partial u}{\partial t}=5\,\frac{\partial^2 u}{\partial x^2}+2x^2,\qquad u=u(t,x).$$

Как классифицируют (метод курса)

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

$$a_{11}\frac{\partial^2 u}{\partial \xi_1^2}+2a_{12}\frac{\partial^2 u}{\partial \xi_1\,\partial \xi_2}+a_{22}\frac{\partial^2 u}{\partial \xi_2^2}+\dots=0,\qquad D=a_{12}^2-a_{11}a_{22}.$$

• $D<0$ — эллиптический тип (обе вторые производные одного знака, как у уравнения Лапласа);
• $D=0$ — параболический тип (одна из старших производных отсутствует, как у уравнения теплопроводности);
• $D>0$ — гиперболический тип (как у волнового уравнения).

Применяем к нашему уравнению

Независимые переменные — $t$ и $x$. Выпишем все вторые производные. Присутствует только $\dfrac{\partial^2 u}{\partial x^2}$ с коэффициентом $5$; вторых производных по $t$ и смешанной нет:

$$a_{tt}=0,\qquad a_{tx}=0,\qquad a_{xx}=5.$$

Член $0{,}2\,\dfrac{\partial u}{\partial t}$ — производная первого порядка, в дискриминант старших производных она не входит. Аналогично свободный член $2x^2$ (правая часть) на тип не влияет. Тогда

$$D=a_{tx}^2-a_{tt}\,a_{xx}=0^2-0\cdot 5=0.$$

Вывод

Поскольку $D=0$ — уравнение параболического типа. Это видно и структурно: первая производная по времени плюс одна вторая производная по координате — это классическое уравнение теплопроводности (диффузии)

$$\frac{\partial u}{\partial t}=a^2\frac{\partial^2 u}{\partial x^2}+f(x),\qquad a^2=\frac{5}{0{,}2}=25,$$

с источником $f(x)=\dfrac{2x^2}{0{,}2}=10x^2$. Члены первого порядка и правая часть тип не меняют — он определяется исключительно старшими (вторыми) производными.

Что требуется

Определить порядок аппроксимации в точке $(t^n,x_j)$ схемы

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}+8\,\frac{u_j^n-u_{j-1}^n}{h}=e^{\,n\Delta t\,(j-1)h},$$

аппроксимирующей $\dfrac{\partial u}{\partial t}+8\dfrac{\partial u}{\partial x}=e^{tx}$.

Метод

Подставляем точное гладкое решение $u(t,x)$ и раскладываем сеточные значения $u_j^{n+1}$, $u_{j-1}^n$ в ряд Тейлора у узла $(t^n,x_j)$. Разность разностного и дифференциального операторов — невязка $\psi$; её главный член (наименьшие степени $\Delta t$, $h$) задаёт порядок, отдельно по каждой переменной. Сетка равномерная: $t^n=n\Delta t$, $x_j=(j-1)h$, поэтому источник $e^{\,n\Delta t(j-1)h}=e^{t^n x_j}$ совпадает с $e^{tx}$ в узле — аппроксимирован точно, в $\psi$ не входит.

Время — правая (вперёд) разность

$$u_j^{n+1}=u_j^n+\left.\frac{\partial u}{\partial t}\right|_j^n\Delta t+\frac{1}{2}\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t^2+\dots$$$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=\left.\frac{\partial u}{\partial t}\right|_j^n+\underbrace{\frac{1}{2}\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t}_{O(\Delta t)}+\dots$$

Первый порядок по $\Delta t$.

Координата — левая (назад) разность

$u_{j-1}^n=u(t^n,x_j-h)$, знаки чередуются:

$$u_{j-1}^n=u_j^n-\left.\frac{\partial u}{\partial x}\right|_j^n h+\frac{1}{2}\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n h^2-\dots$$$$8\,\frac{u_j^n-u_{j-1}^n}{h}=8\left.\frac{\partial u}{\partial x}\right|_j^n-\underbrace{4h\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n}_{O(h)}+\dots$$

Первый порядок по $h$.

Невязка и вывод

Подстановка даёт $u_t+8u_x+\psi=e^{t^n x_j}$, где

$$\psi=\frac{\Delta t}{2}\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n-4h\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n+\dots=O(\Delta t)+O(h).$$

При $\Delta t\to0$, $h\to0$ имеем $\psi\to0$ — схема согласована. Итог: первый порядок по времени и первый по координате, $\psi=O(\Delta t+h)$.

Что требуется

Методом гармоник исследовать устойчивость явной разностной схемы

$$0{,}2\,\frac{u_j^{n+1}-u_j^n}{\Delta t}=5\,\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}+2h^2(j-1)^2,$$

аппроксимирующей уравнение $0{,}2\,\dfrac{\partial u}{\partial t}=5\,\dfrac{\partial^2 u}{\partial x^2}+2x^2$ в точке $(t^n,x_j)$.

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

Погрешность округления удовлетворяет той же однородной разностной схеме, что и сама сеточная функция. Поэтому свободный член $2h^2(j-1)^2$ на устойчивость не влияет (он не содержит искомую функцию $u$) и отбрасывается. В оставшуюся однородную схему подставляем гармонику

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

где $\lambda$ — множитель перехода за один шаг по времени, $\alpha=kh$ — фаза. Необходимое условие устойчивости (равномерной ограниченности гармоник при $n\to\infty$):

$$|\lambda|\le 1\quad\text{при всех }\alpha.$$

Схема явная: пространственный оператор взят на нижнем слое $n$.

Приведение к стандартному виду

Делим обе части (уже без источника) на коэффициент $0{,}2$ при производной по времени:

$$\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},\qquad \sigma=\frac{5}{0{,}2}=25.$$

Здесь $\sigma$ — эффективный коэффициент диффузии.

Подстановка гармоники

Подставляем $u_j^n=\lambda^n e^{i\alpha j}$. Тогда $u_j^{n+1}=\lambda^{n+1}e^{i\alpha j}$, $u_{j\pm1}^n=\lambda^n e^{i\alpha(j\pm1)}=\lambda^n e^{i\alpha j}e^{\pm i\alpha}$. Делим всё уравнение на $\lambda^n e^{i\alpha j}$:

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

Используем тождество $e^{i\alpha}+e^{-i\alpha}=2\cos\alpha$ и формулу понижения $2\cos\alpha-2=-4\sin^2\dfrac{\alpha}{2}$:

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

Отсюда множитель перехода

$$\frac{\lambda-1}{\Delta t}=-\frac{4\sigma}{h^2}\sin^2\frac{\alpha}{2}\quad\Longrightarrow\quad \boxed{\;\lambda(\alpha)=1-\frac{4\sigma\,\Delta t}{h^2}\sin^2\frac{\alpha}{2}=1-\frac{100\,\Delta t}{h^2}\sin^2\frac{\alpha}{2}.\;}$$

Анализ условия $|\lambda|\le 1$

Величина $\lambda$ вещественна, поэтому условие $|\lambda|\le1$ распадается на два неравенства $-1\le\lambda\le1$. Обозначим

$$\mu=\frac{4\sigma\,\Delta t}{h^2}\sin^2\frac{\alpha}{2}\ge 0,\qquad \lambda=1-\mu.$$

• Сверху: поскольку $\mu\ge0$, всегда $\lambda=1-\mu\le1$ — выполнено автоматически при любом $\alpha$.
• Снизу: требуется $\lambda\ge-1$, то есть $1-\mu\ge-1$, или $\mu\le2$. Это самое жёсткое условие; его проверяют в наихудшем случае, когда $\sin^2\dfrac{\alpha}{2}$ максимален.

Максимум $\sin^2\dfrac{\alpha}{2}=1$ достигается при $\alpha=\pi$ (самая короткая гармоника, «пилообразная» ошибка). Тогда

$$1-\frac{4\sigma\,\Delta t}{h^2}\ge-1\quad\Longrightarrow\quad \frac{4\sigma\,\Delta t}{h^2}\le2\quad\Longrightarrow\quad \frac{\sigma\,\Delta t}{h^2}\le\frac12.$$

Вывод — условие устойчивости

Подставляя $\sigma=25$:

$$\frac{25\,\Delta t}{h^2}\le\frac12\quad\Longleftrightarrow\quad \boxed{\;\Delta t\le\frac{h^2}{50}.\;}$$

Явная схема условно устойчива: шаг по времени жёстко ограничен квадратом шага по координате (для параболической задачи $\Delta t\sim h^2$). Свободный член $2x^2$ (его сеточный аналог $2h^2(j-1)^2$) на устойчивость не влияет. Для снятия ограничения на шаг применяют неявные схемы, которые абсолютно (безусловно) устойчивы при любом $\Delta t$.

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

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

$$\frac{\partial u}{\partial t}=7\,\frac{\partial^2 u}{\partial x^2}-5t,\qquad u(t{=}0,x)=0,\qquad \frac{\partial u}{\partial x}(t,0)=0,\quad \frac{\partial u}{\partial x}(t,1)=0.$$

Коэффициент диффузии $\sigma=7$, свободный член $f(t)=-5t$ (от $x$ не зависит). Оба граничных условия — 2-го рода (Неймана), задают поток на концах. Вводим сетку: $t^n=n\,\Delta t$ по времени и $x_j=(j-1)h$, $j=1,2,\dots,N$ по пространству, где $h=\dfrac{1}{N-1}$; узловые значения $u_j^n\approx u(t^n,x_j)$. Узлы $j=1$ и $j=N$ — границы $x=0$ и $x=1$.

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

Производная по времени аппроксимируется разностью между слоями $n$ и $n{+}1$, а вторая производная по $x$ берётся как полусумма (центр $t^{n+1/2}$) на обоих слоях; свободный член вычисляется в полуцелой точке $t^{n+1/2}=\big(n+\tfrac12\big)\Delta t$. Схема имеет порядок аппроксимации $O(\Delta t^2,\;h^2)$ и абсолютно устойчива:

$$\frac{u_j^{n+1}-u_j^{n}}{\Delta t}=\frac{7}{2}\cdot\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}+\frac{7}{2}\cdot\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}-5\Big(n+\tfrac12\Big)\Delta t.$$

Обозначим для краткости $r=\dfrac{7\,\Delta t}{2h^2}$.

2. Приведение к трёхдиагональному виду

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

$$A_j\,u_{j-1}^{n+1}+C_j\,u_j^{n+1}+B_j\,u_{j+1}^{n+1}=F_j,$$

с коэффициентами (конвективного члена нет, поэтому $A_j=B_j$):

$$A_j=B_j=-\frac{7\,\Delta t}{2h^2}=-r,\qquad C_j=1+\frac{7\,\Delta t}{h^2}=1+2r,$$ $$F_j=u_j^{n}+\frac{7\,\Delta t}{2h^2}\big(u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}\big)-5\Big(n+\tfrac12\Big)\Delta t\cdot\Delta t.$$

Система записана для внутренних узлов $j=2,\dots,N-1$; замыкается двумя граничными соотношениями (пп. 5–6).

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

Достаточное условие устойчивости прогонки — диагональное преобладание $|C_j|\ge|A_j|+|B_j|$:

$$|A_j|+|B_j|=\frac{7\,\Delta t}{2h^2}+\frac{7\,\Delta t}{2h^2}=\frac{7\,\Delta t}{h^2}=2r,$$ $$|C_j|=1+\frac{7\,\Delta t}{h^2}=1+2r>2r=|A_j|+|B_j|.$$

Неравенство строгое и выполняется при любых $\Delta t>0,\;h>0$. Значит, знаменатели прогоночных коэффициентов не обращаются в ноль и метод прогонки устойчив.

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

Ищем решение в форме прямого хода $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$. Подставляя $u_{j-1}^{n+1}=\alpha_{j-1}u_j^{n+1}+\beta_{j-1}$ в трёхдиагональное уравнение и выражая $u_j^{n+1}$ через $u_{j+1}^{n+1}$, получаем рекуррентные формулы для коэффициентов:

$$\alpha_j=\frac{-B_j}{C_j+A_j\,\alpha_{j-1}},\qquad \beta_j=\frac{F_j-A_j\,\beta_{j-1}}{C_j+A_j\,\alpha_{j-1}},\qquad j=2,\dots,N-1.$$

С учётом $A_j=B_j=-r$, $C_j=1+2r$:

$$\alpha_j=\frac{r}{1+2r-r\,\alpha_{j-1}},\qquad \beta_j=\frac{F_j+r\,\beta_{j-1}}{1+2r-r\,\alpha_{j-1}}.$$

5. Левое граничное условие — $\alpha_1,\ \beta_1$

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

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

Сравнивая с формой прямого хода $u_1^{n+1}=\alpha_1\,u_2^{n+1}+\beta_1$, получаем стартовые значения

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

(Эквивалентная запись через фиктивный узел $\alpha_0=1,\ \beta_0=0$ даёт $\alpha_1=\dfrac{-B_1}{C_1+A_1},\ \beta_1=\dfrac{F_1}{C_1+A_1}$; обе формы выражают одно и то же равенство $u_1^{n+1}=u_2^{n+1}$.)

6. Правое граничное условие — $u_N^{n+1}$

Условие $\dfrac{\partial u}{\partial x}(t,1)=0$ аппроксимируем левой (односторонней) разностью на правой границе $j=N$:

$$\frac{u_N^{n+1}-u_{N-1}^{n+1}}{h}=0\;\Longrightarrow\; u_N^{n+1}=u_{N-1}^{n+1}.$$

Подставляем сюда прямой ход $u_{N-1}^{n+1}=\alpha_{N-1}\,u_N^{n+1}+\beta_{N-1}$:

$$u_N^{n+1}=\alpha_{N-1}\,u_N^{n+1}+\beta_{N-1}\;\Longrightarrow\;\boxed{\;u_N^{n+1}=\frac{\beta_{N-1}}{1-\alpha_{N-1}}.\;}$$

Знаменатель $1-\alpha_{N-1}\neq 0$ в силу диагонального преобладания ($0<\alpha_j<1$ для $j\ge 2$).

7. Алгоритм счёта по слоям

1. Начальный слой: $u_j^{0}=0$ при всех $j$ (из НУ).
2. Цикл по $n=0,1,2,\dots$: по известному слою $u_j^{n}$ вычисляем правые части $F_j$ ($j=2,\dots,N-1$).
3. Прямой ход: задаём $\alpha_1=1,\ \beta_1=0$ и по рекуррентным формулам считаем $\alpha_j,\beta_j$ до $j=N-1$.
4. Из правого ГУ находим $u_N^{n+1}=\dfrac{\beta_{N-1}}{1-\alpha_{N-1}}$.
5. Обратный ход: $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$ для $j=N-1,\dots,1$.

Шаг $\Delta t$ ограничен только требуемой точностью ($O(\Delta t^2)$), но не устойчивостью — схема Кранка–Николсона абсолютно устойчива.

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

Дано уравнение параболического типа с конвекцией и реакцией:

$$\frac{\partial u}{\partial t}+0{,}8\,\frac{\partial u}{\partial x}=3\,\frac{\partial^2 u}{\partial x^2}-2u,\qquad u=u(t,x),\quad x\in[0,1].$$

Начальное условие $u(t{=}0,x)=0$; граничные условия 1-го рода (Дирихле):

$$\begin{cases} u(t,x{=}0)=t,\\[2pt] u(t,x{=}1)=t^2. \end{cases}$$

Физический смысл коэффициентов: скорость переноса $v=0{,}8>0$, коэффициент диффузии $\sigma=3$, реакционный (поглощающий) член $-2u$.

Вводим равномерную разностную сетку: $n$ — номер слоя по времени, $j$ — номер узла по координате; шаги $\Delta t$ и $h$; $t^{n}=n\Delta t$, $x_j=(j-1)h$, $u_j^{n}=u(t^n,x_j)$, $j=1,\dots,N_x$.

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

В явной схеме все пространственные операторы берутся на известном слое $n$. Аппроксимации:

• производная по времени — правая (вперёд) разность: $\dfrac{\partial u}{\partial t}\approx\dfrac{u_j^{n+1}-u_j^{n}}{\Delta t}$;
• вторая производная — центральная разность: $\dfrac{\partial^2 u}{\partial x^2}\approx\dfrac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}$;
• конвективный член при $v=0{,}8>0$ — левая разность («против потока», для устойчивости): $\dfrac{\partial u}{\partial x}\approx\dfrac{u_j^{n}-u_{j-1}^{n}}{h}$;
• реакционный член $-2u$ — на старом слое: $-2u_j^{n}$.

Подставляя, получаем разностное уравнение во внутренних узлах $j=2,\dots,N_x-1$:

$$\frac{u_j^{n+1}-u_j^{n}}{\Delta t}+0{,}8\,\frac{u_j^{n}-u_{j-1}^{n}}{h}=3\,\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}-2u_j^{n}.$$

Начальное условие на сетке: $u_j^{0}=0$ для всех $j$.

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

$$u_1^{n+1}=t^{n+1}=(n{+}1)\Delta t,\qquad u_{N_x}^{n+1}=\big(t^{n+1}\big)^2=\big((n{+}1)\Delta t\big)^2.$$

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

В разностном уравнении единственное неизвестное слоя $(n{+}1)$ — это $u_j^{n+1}$. Умножаем на $\Delta t$ и выражаем его явно, группируя коэффициенты при $u_{j-1}^{n},u_j^{n},u_{j+1}^{n}$:

$$\boxed{\,u_j^{n+1}=\Big(\frac{3\Delta t}{h^2}+\frac{0{,}8\,\Delta t}{h}\Big)u_{j-1}^{n}+\Big(1-\frac{6\Delta t}{h^2}-\frac{0{,}8\,\Delta t}{h}-2\Delta t\Big)u_j^{n}+\frac{3\Delta t}{h^2}\,u_{j+1}^{n}\,}$$

Введя сеточные параметры $r=\dfrac{3\Delta t}{h^2}$ (диффузионное число) и $c=\dfrac{0{,}8\,\Delta t}{h}$ (число Куранта), запишем компактно:

$$u_j^{n+1}=(r+c)\,u_{j-1}^{n}+\big(1-2r-c-2\Delta t\big)\,u_j^{n}+r\,u_{j+1}^{n}.$$

Соотношение явное: правая часть содержит только значения известного слоя $n$, поэтому новое значение в каждом узле вычисляется напрямую, без решения системы уравнений.

Порядок аппроксимации: схема имеет первый порядок по времени и по пространству, $O(\Delta t+h)$. Первый порядок по $h$ обусловлен противопотоковой (левой) разностью для конвективного члена; диффузионный оператор сам по себе аппроксимирован со вторым порядком.

3. Условие устойчивости (без вывода)

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

$$\frac{\Delta t}{h^2}\le\frac{1}{2\sigma}=\frac{1}{2\cdot 3}=\frac{1}{6}\quad\Longleftrightarrow\quad \boxed{\,\Delta t\le\frac{h^2}{6}\,}$$

Дополнительно для устойчивости конвекции желательно $c=\dfrac{0{,}8\,\Delta t}{h}\le 1$. Шаги $\Delta t$ и $h$ выбираются так, чтобы оба условия выполнялись.

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

1. Задать сетку: $h$, число узлов $N_x$ ($x_j=(j-1)h$), шаг $\Delta t$ с учётом условия устойчивости $\Delta t\le h^2/6$; вычислить $r=3\Delta t/h^2$, $c=0{,}8\,\Delta t/h$.
2. Инициализация (НУ): положить $u_j^{0}=0$ для всех $j=1,\dots,N_x$.
3. Цикл по слоям $n=0,1,2,\dots$ до нужного времени:
   • задать границы из ГУ: $u_1^{n+1}=(n{+}1)\Delta t$, $\;u_{N_x}^{n+1}=\big((n{+}1)\Delta t\big)^2$;
   • для внутренних узлов $j=2,\dots,N_x-1$ пересчитать значения по рекуррентной формуле: $$u_j^{n+1}=(r+c)\,u_{j-1}^{n}+\big(1-2r-c-2\Delta t\big)\,u_j^{n}+r\,u_{j+1}^{n};$$
   • перейти на следующий слой (массив слоя $n{+}1$ становится текущим).
4. Повторять шаг 3, пока не достигнуто требуемое время счёта; при пересчёте контролировать выполнение условия устойчивости.

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

Что требуется

Определить тип уравнения (эллиптическое / параболическое / гиперболическое) и обосновать выбор. Уравнение:

$$3\,\frac{\partial^2 u}{\partial x^2}+7\,\frac{\partial^2 u}{\partial y^2}=e^{x+y},\qquad u=u(x,y).$$

Как классифицируют (метод курса, вопрос 1, семинар 0)

Тип уравнения 2-го порядка определяется только по старшим (вторым) производным. Общий вид для функции двух переменных:

$$a_{11}\frac{\partial^2 u}{\partial \xi_1^2}+2a_{12}\frac{\partial^2 u}{\partial \xi_1\partial \xi_2}+a_{22}\frac{\partial^2 u}{\partial \xi_2^2}+\dots=0,\qquad D=a_{12}^2-a_{11}a_{22}.$$

• $D<0$ — эллиптический тип (обе вторые производные одного знака, как у уравнения Лапласа);
• $D=0$ — параболический (одна из старших производных отсутствует, как у уравнения теплопроводности);
• $D>0$ — гиперболический (как у волнового уравнения).

Применяем к нашему уравнению

Независимые переменные — $x$ и $y$. Выпишем все вторые производные: присутствуют $\dfrac{\partial^2 u}{\partial x^2}$ (коэффициент $a_{11}=3$) и $\dfrac{\partial^2 u}{\partial y^2}$ (коэффициент $a_{22}=7$); смешанной производной $\dfrac{\partial^2 u}{\partial x\,\partial y}$ нет:

$$a_{11}=3,\qquad a_{12}=0,\qquad a_{22}=7.$$

Первых производных в уравнении нет, а правая часть $e^{x+y}$ — свободный член; ни то, ни другое в дискриминант старших производных не входит. Тогда

$$D=a_{12}^2-a_{11}\,a_{22}=0^2-3\cdot 7=-21<0.$$

Вывод

$D=-21<0$ — уравнение эллиптического типа. Это видно и структурно: обе вторые производные присутствуют и стоят с коэффициентами одного знака ($3>0$ и $7>0$), как у уравнения Пуассона $\dfrac{\partial^2 u}{\partial x^2}+\dfrac{\partial^2 u}{\partial y^2}=f$. Свободный член $e^{x+y}$ и численные значения коэффициентов $3,7$ на тип не влияют — тип определяется только знаком дискриминанта старших производных (важно лишь, что оба коэффициента ненулевые и одного знака).

Для эллиптических уравнений курса далее строят разностную аппроксимацию оператора Лапласа на сетке и решают полученную систему (см. соответствующий семинар и страницу «Инструменты» для проверки порядка аппроксимации).

Что требуется

Определить порядок аппроксимации разностной схемы в точке $(t^n,x_j)$:

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=6\,\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}+\ln\bigl[(j-1)h\bigr],$$ $$\text{аппроксимирующей}\qquad \frac{\partial u}{\partial t}=6\,\frac{\partial^2 u}{\partial x^2}+\ln x.$$

Это уравнение теплопроводности с источником; схема явная (по координате значения берутся с известного слоя $n$).

Метод (вопрос 5, семинар 2)

Подставляем в схему точное (гладкое) решение $u(t,x)$ и раскладываем сеточные значения в ряд Тейлора относительно узла $(t^n,x_j)$. Разность между разностным и дифференциальным операторами — это ошибка аппроксимации $\psi$; её главный (с наименьшими степенями $\Delta t$ и $h$) член задаёт порядок. На сетке $t^n=n\Delta t$, $x_j=(j-1)h$, поэтому источник $\ln[(j-1)h]=\ln x_j$ совпадает с $\ln x$ в узле — он аппроксимирован точно и в $\psi$ не входит.

Производная по времени — правая (вперёд) разность

Раскладываем значение на верхнем слое:

$$u_j^{n+1}=u_j^n+\left.\frac{\partial u}{\partial t}\right|_j^n\Delta t+\frac12\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t^2+\dots$$ $$\Rightarrow\quad \frac{u_j^{n+1}-u_j^n}{\Delta t}=\left.\frac{\partial u}{\partial t}\right|_j^n+\underbrace{\frac12\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t}_{O(\Delta t)}+\dots$$

Разность по времени вперёд даёт первый порядок по $\Delta t$.

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

Раскладываем оба соседних узла (члены при нечётных степенях $h$ взаимно сокращаются):

$$u_{j\pm1}^n=u_j^n\pm\left.\frac{\partial u}{\partial x}\right|_j^n h+\frac12\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n h^2\pm\frac16\left.\frac{\partial^3 u}{\partial x^3}\right|_j^n h^3+\frac1{24}\left.\frac{\partial^4 u}{\partial x^4}\right|_j^n h^4\pm\dots$$

Сложив $u_{j+1}^n+u_{j-1}^n$ и вычтя $2u_j^n$, получаем:

$$u_{j+1}^n-2u_j^n+u_{j-1}^n=\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n h^2+\frac1{12}\left.\frac{\partial^4 u}{\partial x^4}\right|_j^n h^4+\dots$$ $$\Rightarrow\quad 6\,\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}=6\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n+\underbrace{6\cdot\frac{h^2}{12}\left.\frac{\partial^4 u}{\partial x^4}\right|_j^n}_{O(h^2)}+\dots$$

Центральная разность второй производной даёт второй порядок по $h$ (первый ненулевой остаток — при $h^2$).

Источник

$\ln[(j-1)h]=\ln x_j=\left.\ln x\right|_j^n$ — берётся в узле без погрешности, вклада в $\psi$ нет.

Ошибка аппроксимации и вывод

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

$$\psi=\frac{\Delta t}{2}\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n-\frac{h^2}{2}\left.\frac{\partial^4 u}{\partial x^4}\right|_j^n+\dots=O(\Delta t)+O(h^2).$$

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

$$\boxed{\;\psi=O(\Delta t)+O(h^2).\;}$$

(Разложения можно проверить интерактивно на странице «Инструменты».)

Что требуется

Методом гармоник исследовать устойчивость явной разностной схемы

$$\dfrac{u_j^{n+1}-u_j^n}{\Delta t} - 25\,\dfrac{u_j^n - u_{j-1}^n}{h} = e^{5n\Delta t\,(j-1)h},$$

аппроксимирующей уравнение переноса $\dfrac{\partial u}{\partial t} - 25\,\dfrac{\partial u}{\partial x} = e^{5tx}$ в точке $(t^n, x_j)$.

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

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

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

где $\lambda$ — множитель перехода за один шаг по времени, $\alpha$ — фаза. Необходимое условие устойчивости — $|\lambda|\le 1$ при всех $\alpha$. Схема явная: пространственная разность берётся на нижнем слое $n$.

Стандартный вид

Перепишем схему как явную схему уравнения переноса $\dfrac{\partial u}{\partial t}+v\,\dfrac{\partial u}{\partial x}=f$. Член $-25\,(u_j^n-u_{j-1}^n)/h$ соответствует $v\,(u_j^n-u_{j-1}^n)/h$ с

$$v=-25<0,$$

то есть характеристическая скорость отрицательна, а пространственная производная аппроксимирована левой (назад) разностью. Свободный член $e^{5n\Delta t\,(j-1)h}$ искомой функции не содержит — отбрасываем.

Подстановка гармоники

Подставляем $u_j^n=\lambda^n e^{i\alpha j}$ в однородную схему и делим на $\lambda^n e^{i\alpha j}$. Используем $u_{j-1}^n\to \lambda^n e^{i\alpha(j-1)}$, то есть множитель $e^{-i\alpha}$:

$$\frac{\lambda-1}{\Delta t} - 25\,\frac{1-e^{-i\alpha}}{h}=0 \quad\Rightarrow\quad \lambda=1+25\,\frac{\Delta t}{h}\bigl(1-e^{-i\alpha}\bigr).$$

Анализ $|\lambda|\le 1$

$\lambda$ комплексно, поэтому условие $|\lambda|\le 1$ означает, что точки $\lambda$ должны лежать в круге радиуса $1$ с центром в нуле. Обозначим

$$r=25\,\frac{\Delta t}{h}>0,\qquad \lambda=1+r - r\,e^{-i\alpha}.$$

Это окружность с центром в точке $(1+r,0)$ и радиусом $r$ (так как $|{-r}\,e^{-i\alpha}|=r$). Центр смещён вправо от нуля на $1+r$, поэтому при любом $r>0$ окружность целиком лежит вне единичного круга: ближайшая к нулю точка $(1+r)-r=1$, самая дальняя $(1+r)+r=1+2r>1$. Значит

$$\max_{\alpha}|\lambda|=1+2r=1+\frac{50\,\Delta t}{h}>1\qquad\text{при любом }\Delta t>0.$$

То же видно из модуля напрямую: $|\lambda|^2=(1+r-r\cos\alpha)^2+(r\sin\alpha)^2=1+2r(1+r)(1-\cos\alpha)\ge 1$, с равенством лишь при $\alpha=0$.

Вывод

Условие $|\lambda|\le 1$ не выполняется ни при каком соотношении шагов $\Delta t$ и $h$.

$$\boxed{\;\text{Схема абсолютно (безусловно) неустойчива при любых }\Delta t,\,h.\;}$$

Причина — неверная сторона разности относительно направления переноса: при $v=-25<0$ характеристики идут влево, и устойчивую явную схему даёт правая (вперёд) разность $(u_{j+1}^n-u_j^n)/h$ с условием Куранта $0<-v\,\Delta t/h\le 1$, то есть $\Delta t\le h/25$. Использованная же левая разность при $v<0$ устойчивости не имеет ни при каком шаге. (Проверить можно на странице «Инструменты».)

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

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

$$\frac{\partial u}{\partial t}+3\frac{\partial u}{\partial x}=0{,}1\frac{\partial^2 u}{\partial x^2}+2t\big(x^2-t\big),$$ $$u(t{=}0,x)=0,\qquad \frac{\partial u}{\partial x}(t,0)=2u(t,0)+3,\qquad u(t,1)=5t.$$

Перепишем в каноническом виде $u_t=\sigma\,u_{xx}+C\,u_x+f$, перенеся конвекцию вправо:

$$\frac{\partial u}{\partial t}=0{,}1\,\frac{\partial^2 u}{\partial x^2}-3\,\frac{\partial u}{\partial x}+2t\big(x^2-t\big).$$

Здесь $\sigma=0{,}1$, конвективный коэффициент $C=-3$ (при $+C\,u_x$), источник $f(t,x)=2t(x^2-t)$. Левое ГУ — 3-го рода, правое — 1-го рода. Сетка: $t^n=n\Delta t$, $x_j=(j-1)h$, $j=1\dots N$ (так $x_1=0$, $x_N=1$), $u_j^n=u(t^n,x_j)$.

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

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

$$\frac{u_j^{n+1}-u_j^{n}}{\Delta t}+3\,\frac{u_{j+1}^{n+1}-u_{j-1}^{n+1}}{2h}=0{,}1\,\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^{2}}+2t^{n+1}\big(x_j^{2}-t^{n+1}\big).$$

2. Приведение к трёхдиагональному виду

Умножаем на $\Delta t$, все неизвестные слоя $(n{+}1)$ собираем слева в виде $A_j u_{j-1}^{n+1}+C_j u_j^{n+1}+B_j u_{j+1}^{n+1}=F_j$. Член $+C\,u_x$ с $C=-3$ даёт добавки: при $u_{j+1}$ слагаемое $-C\dfrac{\Delta t}{2h}=+\dfrac{3\Delta t}{2h}$, при $u_{j-1}$ слагаемое $+C\dfrac{\Delta t}{2h}=-\dfrac{3\Delta t}{2h}$. Получаем:

$$A_j=-\frac{0{,}1\,\Delta t}{h^{2}}-\frac{3\Delta t}{2h},\qquad C_j=1+\frac{0{,}2\,\Delta t}{h^{2}},\qquad B_j=-\frac{0{,}1\,\Delta t}{h^{2}}+\frac{3\Delta t}{2h},$$ $$F_j=u_j^{n}+2t^{n+1}\big(x_j^{2}-t^{n+1}\big)\,\Delta t.$$

(Знаки проверены символически: при $+C\,\partial u/\partial x$ неявно у $u_{j+1}$ коэффициент $-C\Delta t/2h$, у $u_{j-1}$ коэффициент $+C\Delta t/2h$.)

3. Сходимость прогонки

Достаточное условие — диагональное преобладание $|C_j|\ge|A_j|+|B_j|$. Член $A_j$ всегда отрицателен, $|A_j|=\dfrac{0{,}1\Delta t}{h^2}+\dfrac{3\Delta t}{2h}$. Знак $B_j=-\dfrac{0{,}1\Delta t}{h^2}+\dfrac{3\Delta t}{2h}=\dfrac{(15h-1)\Delta t}{10h^2}$ зависит от шага: $B_j\ge0$ при $h\ge 1/15$ и $B_j\le0$ при $h\le 1/15$.

• При мелкой пространственной сетке ($h\le 1/15$, когда $\dfrac{3}{2h}\le\dfrac{0{,}1}{h^2}$) имеем $B_j\le0$, тогда $|A_j|+|B_j|=\Big(\dfrac{0{,}1\Delta t}{h^2}+\dfrac{3\Delta t}{2h}\Big)+\Big(\dfrac{0{,}1\Delta t}{h^2}-\dfrac{3\Delta t}{2h}\Big)=\dfrac{0{,}2\Delta t}{h^{2}}<1+\dfrac{0{,}2\Delta t}{h^{2}}=|C_j|$ — преобладание строгое и безусловное при любых $\Delta t,h$. Это и есть условие отсутствия осцилляций от центральной аппроксимации конвекции (число Пекле $|C|h\le2\sigma$, т.е. $3h\le0{,}2$, $h\le1/15$).
• При крупном шаге ($h\ge 1/15$, $B_j\ge0$) имеем $|A_j|+|B_j|=\dfrac{3\Delta t}{h}$, и условие $1+\dfrac{0{,}2\Delta t}{h^2}\ge\dfrac{3\Delta t}{h}$ выполняется лишь дополнительно, при ограничении на $\Delta t$; центральная аппроксимация конвекции при этом может давать осцилляции.

Таким образом, при сеточном условии $h\le0{,}2/3=1/15$ прогонка устойчива (диагональное преобладание строгое) при любом $\Delta t$.

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

Решение ищем в форме $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$. Подставив $u_{j-1}^{n+1}=\alpha_{j-1}u_j^{n+1}+\beta_{j-1}$ в трёхдиагональное уравнение, получаем прямой ход (от $j=1$ к $j=N-1$):

$$\alpha_j=\frac{-B_j}{C_j+A_j\,\alpha_{j-1}},\qquad \beta_j=\frac{F_j-A_j\,\beta_{j-1}}{C_j+A_j\,\alpha_{j-1}}.$$

5. Левое граничное условие (3-го рода) — $\alpha_1,\beta_1$

Условие $\dfrac{\partial u}{\partial x}(t,0)=2u(t,0)+3$ аппроксимируем правой разностью на левой границе ($x_1=0$):

$$\frac{u_2^{n+1}-u_1^{n+1}}{h}=2u_1^{n+1}+3.$$

Выразим левый узел через соседний $u_1^{n+1}=\alpha_1 u_2^{n+1}+\beta_1$:

$$u_2^{n+1}-u_1^{n+1}=2h\,u_1^{n+1}+3h\;\Rightarrow\;u_1^{n+1}=\frac{u_2^{n+1}-3h}{1+2h}.$$

Отсюда стартовые коэффициенты:

$$\boxed{\;\alpha_1=\frac{1}{1+2h},\qquad \beta_1=-\frac{3h}{1+2h}.\;}$$

6. Правое граничное условие (1-го рода) — $u_N^{n+1}$

Условие $u(t,1)=5t$ задано прямо в узле $x_N=1$, поэтому

$$\boxed{\;u_N^{n+1}=5\,t^{n+1}=5(n+1)\Delta t.\;}$$

Это известное значение служит точкой старта обратного хода $u_j^{n+1}=\alpha_j u_{j+1}^{n+1}+\beta_j$ от $j=N-1$ к $j=1$.

7. Алгоритм

1. $u^0=0$ (начальное условие).
2. Цикл по слоям $n$: вычислить $F_j$ при всех внутренних $j$.
3. Прямой ход: $\alpha_1,\beta_1$ из левого ГУ 3-го рода, далее $\alpha_j,\beta_j$ по рекуррентным формулам до $j=N-1$.
4. Правая граница: $u_N^{n+1}=5(n+1)\Delta t$ (ГУ 1-го рода).
5. Обратный ход $u_j^{n+1}=\alpha_j u_{j+1}^{n+1}+\beta_j$ от $j=N-1$ до $j=1$.
6. Повторять по $n$.

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

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

Дано параболическое уравнение теплопроводности с источником:

$$\frac{\partial u}{\partial t}=0{,}02\,\frac{\partial^2 u}{\partial x^2}+e^{-t},\qquad u=u(t,x),$$ $$u(t{=}0,x)=0,\qquad u(t,x{=}0)=t,\qquad u(t,x{=}1)=2t.$$

Коэффициент диффузии $\sigma=0{,}02>0$, конвективных (первых) производных нет, свободный член $f(t)=e^{-t}$. Оба граничных условия — 1-го рода (Дирихле). Вводим равномерную сетку: по времени $t^n=n\Delta t$, по пространству $x_j=(j-1)h$, $j=1,\dots,N_x$, $h=\dfrac{1}{N_x-1}$.

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

В неявной схеме все пространственные операторы берутся на искомом слое $n{+}1$. Производная по времени — правая (односторонняя) разность, вторая производная — центральная разность на слое $n{+}1$, свободный член — на слое $n{+}1$:

$$\frac{u_j^{n+1}-u_j^{n}}{\Delta t}=0{,}02\,\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}+e^{-t^{n+1}},\qquad t^{n+1}=(n{+}1)\Delta t.$$

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

2. Рекуррентное соотношение (трёхдиагональный вид)

Умножаем на $\Delta t$ и переносим все значения слоя $(n{+}1)$ в левую часть, известные ($n$-й слой и источник) — в правую. Введём безразмерный параметр

$$r=\frac{0{,}02\,\Delta t}{h^2}.$$

Тогда

$$-r\,u_{j-1}^{n+1}+\bigl(1+2r\bigr)u_j^{n+1}-r\,u_{j+1}^{n+1}=u_j^{n}+\Delta t\,e^{-t^{n+1}}.$$

Это трёхдиагональная СЛАУ вида $c_j u_{j-1}^{n+1}+b_j u_j^{n+1}+a_j u_{j+1}^{n+1}=\xi_j$ с коэффициентами

$$a_j=c_j=-r=-\frac{0{,}02\,\Delta t}{h^2},\qquad b_j=1+2r=1+\frac{0{,}04\,\Delta t}{h^2},\qquad \xi_j=u_j^{n}+\Delta t\,e^{-t^{n+1}}.$$

Сходимость прогонки (диагональное преобладание): $|a_j|+|c_j|=2r<1+2r=|b_j|$ выполнено при любых $\Delta t,h$ — схема абсолютно устойчива (безусловно устойчива), ограничения на шаг по времени нет. Анализ фон Неймана подтверждает: множитель перехода $\lambda=\dfrac{1}{1+4r\sin^2(\varphi/2)}$, так что $0<\lambda\le 1$ при любом $r\ge 0$.

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

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

3. Учёт граничных условий (прогоночные коэффициенты на краях)

Оба ГУ — 1-го рода, значения на границах заданы явно, поэтому прогоночные коэффициенты на левом крае тривиальны:

• Левая граница $j=1$: $u_1^{n+1}=t^{n+1}=(n{+}1)\Delta t\;\Rightarrow\;\alpha_1=0,\ \beta_1=(n{+}1)\Delta t.$
• Правая граница $j=N_x$: $u_{N_x}^{n+1}=2t^{n+1}=2(n{+}1)\Delta t$ — стартовое значение для обратного хода прогонки.

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

1. Инициализация (НУ): $u_j^{0}=0$ для всех узлов $j=1,\dots,N_x$ (начальное условие $u(0,x)=0$).
2. Цикл по слоям $n=0,1,2,\dots$ (до достижения конечного времени):
   1. вычислить источник $\Delta t\,e^{-t^{n+1}}$ и правые части $\xi_j=u_j^{n}+\Delta t\,e^{-t^{n+1}}$ для внутренних узлов $j=2,\dots,N_x-1$;
   2. прямой ход прогонки: положить $\alpha_1=0,\ \beta_1=(n{+}1)\Delta t$ (ЛГУ) и для $j=2,\dots,N_x-1$ пересчитать $\alpha_j,\beta_j$ по формулам п. 2 с $a_j=c_j=-r$, $b_j=1+2r$;
   3. задать значение на правой границе из ГУ: $u_{N_x}^{n+1}=2(n{+}1)\Delta t$;
   4. обратный ход прогонки: для $j=N_x-1,\dots,1$ вычислить $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$;
   5. зафиксировать границы $u_1^{n+1}=(n{+}1)\Delta t$, $u_{N_x}^{n+1}=2(n{+}1)\Delta t$ и перейти к следующему слою ($u^{n}\leftarrow u^{n+1}$).

Поскольку коэффициенты $a_j,b_j,c_j$ постоянны (не зависят от $n$ и $j$), на каждом слое решается одна трёхдиагональная система методом прогонки; пересчёта по устойчивости не требуется — схема безусловно устойчива.

Порядок аппроксимации схемы: $O\!\left(\Delta t,\,h^2\right)$ (первый по времени, второй по пространству).

Что требуется

Определить тип уравнения (эллиптическое / параболическое / гиперболическое) и обосновать выбор. Уравнение:

$$3\,\frac{\partial^2 u}{\partial x^2}-2\,\frac{\partial^2 u}{\partial y^2}=xy,\qquad u=u(x,y).$$

Как классифицируют (метод курса, вопрос 1, семинар 0)

Тип линейного уравнения 2-го порядка определяется по старшим (вторым) производным. Общий вид для функции двух переменных $\xi_1,\xi_2$:

$$a_{11}\frac{\partial^2 u}{\partial \xi_1^2}+2a_{12}\frac{\partial^2 u}{\partial \xi_1\partial \xi_2}+a_{22}\frac{\partial^2 u}{\partial \xi_2^2}+\dots=0,\qquad D=a_{12}^2-a_{11}a_{22}.$$

• $D<0$ — эллиптический тип (обе вторые производные одного знака, как у уравнения Лапласа);
• $D=0$ — параболический (одна из старших производных отсутствует, как у уравнения теплопроводности);
• $D>0$ — гиперболический (производные разных знаков, как у волнового уравнения).

Применяем к нашему уравнению

Независимые переменные — $x$ и $y$. Выпишем все вторые производные. Присутствуют $\dfrac{\partial^2 u}{\partial x^2}$ с коэффициентом $3$ и $\dfrac{\partial^2 u}{\partial y^2}$ с коэффициентом $-2$; смешанной производной $\dfrac{\partial^2 u}{\partial x\,\partial y}$ нет:

$$a_{xx}=3,\qquad a_{xy}=0,\qquad a_{yy}=-2.$$

Считаем дискриминант старших производных:

$$D=a_{xy}^2-a_{xx}\,a_{yy}=0^2-3\cdot(-2)=6>0.$$

Та же проверка в стандартной записи $A\,u_{xx}+B\,u_{xy}+C\,u_{yy}$ с $A=3,\ B=0,\ C=-2$ даёт $B^2-4AC=0-4\cdot3\cdot(-2)=24>0$ — знак тот же, вывод тот же.

Вывод

$D=6>0$ — уравнение гиперболического типа. Это видно и структурно: две вторые производные входят с коэффициентами разных знаков ($+3$ при $\partial^2 u/\partial x^2$ и $-2$ при $\partial^2 u/\partial y^2$), то есть $3\,u_{xx}-2\,u_{yy}=\dots$ — это «волновая» структура $u_{\tau\tau}-c^2 u_{\eta\eta}$ (после нормировки $\tau=x/\sqrt3$, $\eta=y/\sqrt2$ левая часть принимает вид $u_{\tau\tau}-u_{\eta\eta}$, а правая переходит в $xy=\sqrt6\,\tau\eta$).

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

Для гиперболических уравнений курса далее применяют соответствующие разностные схемы; проверить порядок аппроксимации и устойчивость можно интерактивно на странице «Инструменты».

Что требуется

Определить порядок аппроксимации разностной схемы в точке $(t^n,x_j)$:

$$\dfrac{u_j^{n+1}-u_j^n}{\Delta t}+9\,\dfrac{u_{j+1}^n-u_j^n}{h}=2\,\dfrac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}+(j-1)^2h^2,$$

аппроксимирующей уравнение

$$\dfrac{\partial u}{\partial t}+9\,\dfrac{\partial u}{\partial x}=2\,\dfrac{\partial^2 u}{\partial x^2}+x^2.$$

Метод (вопрос 5, семинар 2)

Подставляем в схему точное гладкое решение $u(t,x)$ дифференциального уравнения и раскладываем сеточные значения в узлах в ряд Тейлора относительно точки $(t^n,x_j)$, где $t^n=n\Delta t$, $x_j=(j-1)h$. Разность между разностным и дифференциальным операторами — это ошибка аппроксимации $\psi$; её главный член (с наименьшими степенями $\Delta t$ и $h$) задаёт порядок схемы отдельно по каждой переменной.

Свободный член

Так как $x_j=(j-1)h$, дискретный источник совпадает с непрерывным точно в узле:

$$(j-1)^2h^2=\big((j-1)h\big)^2=x_j^2.$$

Поэтому источник аппроксимирован точно и в ошибку $\psi$ не входит.

Производная по времени — правая (вперёд) разность

$$u_j^{n+1}=u_j^n+\left.\dfrac{\partial u}{\partial t}\right|_j^n\Delta t+\dfrac{1}{2!}\left.\dfrac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t^2+\dots$$ $$\Rightarrow\quad \dfrac{u_j^{n+1}-u_j^n}{\Delta t}=\left.\dfrac{\partial u}{\partial t}\right|_j^n+\underbrace{\dfrac{1}{2}\left.\dfrac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t}_{O(\Delta t)}+\dots$$

Правая разность по времени даёт первый порядок по $\Delta t$.

Конвективный член — правая (вперёд) разность по $x$

Из разложения $u_{j+1}^n$ относительно узла $x_j$:

$$u_{j+1}^n=u_j^n+\left.\dfrac{\partial u}{\partial x}\right|_j^n h+\dfrac{1}{2!}\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n h^2+\dots$$ $$\Rightarrow\quad 9\,\dfrac{u_{j+1}^n-u_j^n}{h}=9\left.\dfrac{\partial u}{\partial x}\right|_j^n+\underbrace{9\cdot\dfrac{h}{2}\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n}_{O(h)}+\dots$$

Правая разность первой производной по координате даёт первый порядок по $h$.

Диффузионный член — центральная разность второй производной

Складывая разложения $u_{j+1}^n$ и $u_{j-1}^n$, нечётные степени $h$ сокращаются:

$$u_{j+1}^n=u_j^n+\left.\dfrac{\partial u}{\partial x}\right|_j^n h+\dfrac{1}{2!}\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n h^2+\dfrac{1}{3!}\left.\dfrac{\partial^3 u}{\partial x^3}\right|_j^n h^3+\dfrac{1}{4!}\left.\dfrac{\partial^4 u}{\partial x^4}\right|_j^n h^4+\dots$$ $$u_{j-1}^n=u_j^n-\left.\dfrac{\partial u}{\partial x}\right|_j^n h+\dfrac{1}{2!}\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n h^2-\dfrac{1}{3!}\left.\dfrac{\partial^3 u}{\partial x^3}\right|_j^n h^3+\dfrac{1}{4!}\left.\dfrac{\partial^4 u}{\partial x^4}\right|_j^n h^4-\dots$$ $$u_{j+1}^n-2u_j^n+u_{j-1}^n=\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n h^2+\dfrac{2}{4!}\left.\dfrac{\partial^4 u}{\partial x^4}\right|_j^n h^4+\dots$$ $$\Rightarrow\quad 2\,\dfrac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}=2\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n+\underbrace{\dfrac{1}{6}\left.\dfrac{\partial^4 u}{\partial x^4}\right|_j^n h^2}_{O(h^2)}+\dots$$

Центральная аппроксимация второй производной даёт второй порядок по $h$.

Ошибка аппроксимации и вывод

Подставляя разложения в схему и вычитая дифференциальное уравнение, получаем (перенося диффузионный член влево); выписываем младшие по $\Delta t$ и $h$ члены:

$$\psi=\underbrace{\dfrac{\Delta t}{2}\,\dfrac{\partial^2 u}{\partial t^2}}_{O(\Delta t)}+\underbrace{\dfrac{9h}{2}\,\dfrac{\partial^2 u}{\partial x^2}}_{O(h)}+\;O(h^2)+\dots$$

(При $h^2$ конкурируют поправки от диффузии $-\tfrac{1}{6}\partial^4 u/\partial x^4$ и следующий порядок разности вперёд конвекции $+\tfrac{3}{2}\partial^3 u/\partial x^3$ — оба члены $O(h^2)$ и на порядок не влияют.) По координате конкурируют члены $O(h)$ (от правой разности конвекции) и $O(h^2)$ (от центральной разности диффузии); главным является младший — $O(h)$. Итог:

$$\psi=O(\Delta t)+O(h).$$

Схема имеет первый порядок аппроксимации по времени и первый по координате: $\;\psi=O(\Delta t+h)$. (Проверить разложения можно на странице «Инструменты».)

Что требуется

Методом гармоник исследовать устойчивость неявной разностной схемы

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

аппроксимирующей уравнение $\dfrac{\partial u}{\partial t}=9\dfrac{\partial^2 u}{\partial x^2}+x^2$. Здесь $\sigma=9$, а свободный член $(j-1)^2h^2=x_j^2$ — источник, не содержащий искомую функцию.

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

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

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

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

Однородная схема

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=9\,\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}.$$

Подстановка гармоники

Подставляем $u_j^n=\lambda^n e^{i\alpha j}$ и делим на $\lambda^n e^{i\alpha j}$. Производная по времени даёт множитель $\dfrac{\lambda-1}{\Delta t}$, а все пространственные члены — на слое $(n+1)$, поэтому дают общий множитель $\lambda$:

$$\frac{\lambda-1}{\Delta t}=9\,\frac{\lambda\bigl(e^{i\alpha}-2+e^{-i\alpha}\bigr)}{h^2}.$$

Тождество для второй разности

Используем $e^{i\alpha}-2+e^{-i\alpha}=2\cos\alpha-2=-4\sin^2\dfrac{\alpha}{2}$:

$$\frac{\lambda-1}{\Delta t}=-\frac{36\,\lambda}{h^2}\sin^2\frac{\alpha}{2}.$$

Выражаем множитель перехода

Сгруппируем члены с $\lambda$:

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

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

$$\boxed{\;\frac{1}{\lambda}=1+\frac{36\,\Delta t}{h^2}\sin^2\frac{\alpha}{2}.\;}$$

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

Для неявных схем условие $|\lambda|\le 1$ удобно переписать как $\left|\dfrac{1}{\lambda}\right|\ge 1$. Введём обозначение

$$q=\frac{36\,\Delta t}{h^2}\sin^2\frac{\alpha}{2}\ge 0.$$

Тогда $\dfrac{1}{\lambda}=1+q$ — величина вещественная и не меньше единицы при любых $\alpha$, $\Delta t$, $h$ (так как $\sigma=9>0$ и $\sin^2\frac{\alpha}{2}\ge 0$):

$$\frac{1}{\lambda}=1+q\ge 1\quad\Longrightarrow\quad 0<\lambda=\frac{1}{1+q}\le 1.$$

Значит $|\lambda|\le 1$ для всех фаз $\alpha$ — без каких-либо ограничений на шаги. Наихудший случай $\sin^2\frac{\alpha}{2}=1$ лишь уменьшает $\lambda$ (приближает к нулю), но границу $|\lambda|\le 1$ не нарушает.

Вывод

Неявная (чисто неявная, с разностями на слое $n+1$) схема абсолютно (безусловно) устойчива: ограничений на $\Delta t$ и $h$ нет. Это коренное отличие от явной схемы того же уравнения, где появляется жёсткое условие $\sigma\dfrac{\Delta t}{h^2}\le\dfrac12$ (ср. КР1, задача 3). Свободный член $x^2$ на устойчивость не влияет. Проверить условие на числах можно на странице «Инструменты: устойчивость».

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

$$\frac{\partial u}{\partial t}=4\,\frac{\partial^2 u}{\partial x^2}+2x^2-t,\qquad u(t{=}0,x)=0,\qquad \frac{\partial u}{\partial x}(t,0)=3,\quad \frac{\partial u}{\partial x}(t,1)=5\,u(t,1).$$

Здесь $\sigma=4$, свободный член $f(t,x)=2x^2-t$. Левое граничное условие — 2-го рода (Нейман, неоднородный), правое — 3-го рода (Робена). Конвективного члена $\partial u/\partial x$ в уравнении нет, поэтому система симметрична. Сетка: $t^n=n\Delta t$, $x_j=(j-1)h$, $j=1\dots N_x$ (так что $x_1=0$, $x_{N_x}=1$), $u_j^n=u(t^n,x_j)$.

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

Вторую производную берём «пополам» на слоях $n$ и $n{+}1$, свободный член — в точке $t^{n+1/2}$. Схема имеет порядок $O(\Delta t^2,h^2)$ и абсолютно устойчива:

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=\frac{4}{2}\,\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}+\frac{4}{2}\,\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}+2x_j^2-\Big(n+\tfrac12\Big)\Delta t.$$

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

Умножаем на $\Delta t$, неизвестные слоя $(n{+}1)$ переносим налево, известные — направо. Получаем трёхдиагональную систему $a_j u_{j+1}^{n+1}+b_j u_j^{n+1}+c_j u_{j-1}^{n+1}=\xi_j$:

$$a_j=c_j=-\frac{2\Delta t}{h^2},\qquad b_j=1+\frac{4\Delta t}{h^2},$$ $$\xi_j=u_j^n+\frac{2\Delta t}{h^2}\big(u_{j+1}^n-2u_j^n+u_{j-1}^n\big)+\Big(2x_j^2-\big(n+\tfrac12\big)\Delta t\Big)\Delta t.$$

3. Сходимость прогонки

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

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

Выполнено при любых $\Delta t,h$ — прогонка устойчива.

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

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

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

5. Левое граничное условие — $\alpha_1,\beta_1$

Условие $\dfrac{\partial u}{\partial x}(t,0)=3$ (2-го рода, неоднородное) аппроксимируем правой разностью на левой границе ($x_1=0$):

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

Сравнивая с формой прогонки $u_1^{n+1}=\alpha_1 u_2^{n+1}+\beta_1$, получаем стартовые коэффициенты

$$\boxed{\;\alpha_1=1,\qquad \beta_1=-3h.\;}$$

(Неоднородность ГУ ушла именно в $\beta_1$; для однородного условия было бы $\beta_1=0$.)

6. Правое граничное условие — $u_{N_x}^{n+1}$

Условие $\dfrac{\partial u}{\partial x}(t,1)=5\,u(t,1)$ (3-го рода) аппроксимируем левой разностью в узле $x_{N_x}=1$:

$$\frac{u_{N_x}^{n+1}-u_{N_x-1}^{n+1}}{h}=5\,u_{N_x}^{n+1}.$$

Подставляем обратный ход $u_{N_x-1}^{n+1}=\alpha_{N_x-1}u_{N_x}^{n+1}+\beta_{N_x-1}$:

$$u_{N_x}^{n+1}-\alpha_{N_x-1}u_{N_x}^{n+1}-\beta_{N_x-1}=5h\,u_{N_x}^{n+1} \;\Longrightarrow\; u_{N_x}^{n+1}\big(1-\alpha_{N_x-1}-5h\big)=\beta_{N_x-1},$$ $$\boxed{\;u_{N_x}^{n+1}=\frac{\beta_{N_x-1}}{\,1-\alpha_{N_x-1}-5h\,}.\;}$$

Далее обратный ход $u_j^{n+1}=\alpha_j u_{j+1}^{n+1}+\beta_j$ от $j=N_x-1$ до $1$.

7. Алгоритм

Цикл по слоям $n$: вычислить $\xi_j$ $\to$ прямой ход ($\alpha_j,\beta_j$ от $j{=}1$ со стартом $\alpha_1=1,\ \beta_1=-3h$) $\to$ найти $u_{N_x}^{n+1}$ из правого ГУ 3-го рода $\to$ обратный ход. Схема неявная (Кранка–Николсона), решается методом прогонки.

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

Дано уравнение переноса (гиперболическое, 1-го порядка):

$$\frac{\partial u}{\partial t}-0{,}7\,\frac{\partial u}{\partial x}=0,\qquad u(t{=}0,x)=0,\qquad u(t,x{=}1)=t.$$

Запишем его в стандартном виде $\dfrac{\partial u}{\partial t}+a\,\dfrac{\partial u}{\partial x}=0$ со скоростью переноса $a=-0{,}7<0$. Знак скорости задаёт направление «потока»: при $a<0$ возмущения распространяются справа налево (от $x=1$ к $x=0$), поэтому граничное условие задано именно на правом конце $x=1$. Сетка: $t^n=n\Delta t$, $x_j=(j-1)h$, $j=1\dots N_x$ (так что $x_1=0$, $x_{N_x}=1$), $u_j^n=u(t^n,x_j)$.

1. Выбор разностей и неявная схема (вопрос 6, семинар 4)

Производную по времени берём правой (вперёд) разностью. Для конвективного члена применяем разность «против потока»: так как $a<0$, устойчивой является правая разность $\dfrac{u_{j+1}-u_j}{h}$ (узел берётся со стороны, откуда приходит возмущение). В неявной схеме все пространственные операторы относим на новый слой $n{+}1$:

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}-0{,}7\,\frac{u_{j+1}^{n+1}-u_j^{n+1}}{h}=0.$$

Свободного члена и второй производной в уравнении нет, поэтому схема двухточечная по пространству. Неявная схема для уравнения переноса абсолютно устойчива (ограничения типа Куранта на $\Delta t$ нет), порядок аппроксимации $O(\Delta t,h)$.

2. Приведение к прогоночному виду

Умножаем на $\Delta t$ и собираем неизвестные слоя $(n{+}1)$ слева. Введём число Куранта $r=\dfrac{0{,}7\,\Delta t}{h}>0$:

$$u_j^{n+1}-r\big(u_{j+1}^{n+1}-u_j^{n+1}\big)=u_j^n\quad\Longrightarrow\quad (1+r)\,u_j^{n+1}-r\,u_{j+1}^{n+1}=u_j^n.$$

Это система с двухдиагональной матрицей: каждое уравнение связывает $u_j^{n+1}$ только с соседом справа $u_{j+1}^{n+1}$ (член $u_{j-1}^{n+1}$ отсутствует, $c_j=0$). В обозначениях прогонки $a_j u_{j+1}^{n+1}+b_j u_j^{n+1}+c_j u_{j-1}^{n+1}=\xi_j$:

$$a_j=-r,\qquad b_j=1+r,\qquad c_j=0,\qquad \xi_j=u_j^n.$$

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

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

$$|a_j|+|c_j|=r+0=r<1+r=|b_j|.$$

Матрица устойчива, обусловленность гарантирована. Анализ Фурье даёт множитель перехода $\lambda=\dfrac{1}{1+r-r\,e^{i\varphi}}$, причём $\max|\lambda|=1$ при любом $r>0$ — схема абсолютно устойчива.

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

Поскольку $c_j=0$, полный прямой ход прогонки не нужен — неизвестное $u_j^{n+1}$ выражается явно через уже найденного правого соседа $u_{j+1}^{n+1}$ (обратный ход справа налево):

$$\boxed{\,u_j^{n+1}=\frac{u_j^n+r\,u_{j+1}^{n+1}}{1+r}=\frac{u_j^n+\dfrac{0{,}7\,\Delta t}{h}\,u_{j+1}^{n+1}}{1+\dfrac{0{,}7\,\Delta t}{h}}\,}$$

В стандартной форме прогонки $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$ коэффициенты здесь постоянны (не зависят от $j$):

$$\alpha_j=\frac{r}{1+r}=\frac{0{,}7\,\Delta t}{0{,}7\,\Delta t+h},\qquad \beta_j=\frac{u_j^n}{1+r}=\frac{h\,u_j^n}{0{,}7\,\Delta t+h}.$$

Так как $0<\alpha_j<1$, ошибки при обратном ходе не нарастают.

5. Граничное условие (старт обратного хода)

Прогонка стартует с правой границы, где задано ГУ 1-го рода:

$$u_{N_x}^{n+1}=t^{n+1}=(n+1)\Delta t.$$

Левая граница $x=0$ — выходная (возмущение «вытекает»), отдельного условия там нет: значение $u_1^{n+1}$ получается из той же рекуррентной формулы при $j=1$ через $u_2^{n+1}$.

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

1. Инициализация (НУ): $u_j^0=0$ для всех $j=1,\dots,N_x$.
2. Цикл по слоям $n=0,1,2,\dots,N_t-1$:
   • применить правое ГУ: $u_{N_x}^{n+1}=(n+1)\Delta t$;
   • обратный ход (справа налево) для $j=N_x-1,\,N_x-2,\,\dots,1$: $$u_j^{n+1}=\frac{u_j^n+r\,u_{j+1}^{n+1}}{1+r},\qquad r=\frac{0{,}7\,\Delta t}{h};$$
   • перейти к следующему слою (заменить $u_j^n\leftarrow u_j^{n+1}$).
3. Выход по достижении нужного момента времени $t=N_t\,\Delta t$.

Поскольку схема неявная и абсолютно устойчивая, шаг $\Delta t$ выбирается из соображений точности, а не устойчивости. Блок-схема — на странице «Блок-схемы» (неявная схема, обратный ход вместо полной прогонки из-за двухдиагональной матрицы).

Что требуется

Определить тип уравнения (эллиптическое / параболическое / гиперболическое) и обосновать выбор. Уравнение:

$$\dfrac{\partial u}{\partial t} - 3\,\dfrac{\partial u}{\partial x} = 0{,}5\,\dfrac{\partial^2 u}{\partial x^2} - 10xt,\qquad u=u(t,x).$$

Как классифицируют (метод курса, вопрос 1, семинар 0)

Тип линейного уравнения 2-го порядка определяется только по старшим (вторым) производным. Общий вид для функции двух переменных:

$$a_{11}\frac{\partial^2 u}{\partial \xi_1^2}+2a_{12}\frac{\partial^2 u}{\partial \xi_1\partial \xi_2}+a_{22}\frac{\partial^2 u}{\partial \xi_2^2}+\dots=0,\qquad D=a_{12}^2-a_{11}a_{22}.$$

• $D<0$ — эллиптический тип (обе вторые производные одного знака, как у Лапласа);
• $D=0$ — параболический (одна из старших производных отсутствует, как у уравнения теплопроводности);
• $D>0$ — гиперболический (как у волнового уравнения).

Приводим к стандартному виду

Перенесём все слагаемые в одну часть, чтобы выделить коэффициенты при старших производных:

$$\dfrac{\partial u}{\partial t} - 0{,}5\,\dfrac{\partial^2 u}{\partial x^2} - 3\,\dfrac{\partial u}{\partial x} + 10xt = 0.$$

Применяем к нашему уравнению

Независимые переменные — $t$ и $x$. Выпишем все вторые производные: присутствует только $\dfrac{\partial^2 u}{\partial x^2}$ (её коэффициент равен $-0{,}5$); вторых производных по $t$ и смешанной нет:

$$a_{tt}=0,\qquad a_{tx}=0,\qquad a_{xx}=-0{,}5.$$

Член $\dfrac{\partial u}{\partial t}$ — первого порядка по времени, а $-3\,\dfrac{\partial u}{\partial x}$ — снос (конвективный член) первого порядка; в дискриминант старших производных они не входят. Тогда

$$D=a_{tx}^2-a_{tt}\,a_{xx}=0^2-0\cdot(-0{,}5)=0.$$

Вывод

$D=0$ — уравнение параболического типа. Это видно и структурно: первая производная по времени плюс единственная вторая производная по координате — это уравнение конвекции–диффузии (теплопроводности со сносом) с источником $-10xt$. Конвективный член $-3\,\partial u/\partial x$, свободный член $-10xt$ и числовые коэффициенты на тип не влияют — тип определяется исключительно старшими (вторыми) производными.

Для параболических уравнений курса далее применяют явные/неявные разностные схемы (см. семинар 4, страницу «Инструменты» для проверки порядка аппроксимации).

Что требуется

Определить порядок аппроксимации разностной схемы в точке $(t^n,x_j)$:

$$\dfrac{u_j^{n+1}-u_j^n}{\Delta t}+8\,\dfrac{u_{j+1}^n-u_{j-1}^n}{2h}=2n\cdot\Delta t,$$ $$\text{аппроксимирующей}\qquad \dfrac{\partial u}{\partial t}+8\,\dfrac{\partial u}{\partial x}=2t.$$

Метод (вопрос 5, семинар 2)

Подставляем в разностную схему точное гладкое решение $u(t,x)$ и раскладываем сеточные значения в ряд Тейлора относительно узла $(t^n,x_j)$ по $\Delta t$ и $h$. Разность между разностным и дифференциальным операторами — это ошибка (невязка) аппроксимации $\psi$; её главный член (с наименьшими степенями $\Delta t$ и $h$) задаёт порядок. На сетке $t^n=n\Delta t$, поэтому правая часть $2n\cdot\Delta t=2t^n=2t$ в узле совпадает с правой частью уравнения точно — источник аппроксимирован без погрешности и в $\psi$ не входит.

Производная по времени — правая (вперёд) разность

Раскладываем верхний слой $u_j^{n+1}=u(t^n+\Delta t,\,x_j)$:

$$u_j^{n+1}=u_j^n+\left.\dfrac{\partial u}{\partial t}\right|_j^n\Delta t+\dfrac12\left.\dfrac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t^2+\dots$$ $$\Rightarrow\quad \dfrac{u_j^{n+1}-u_j^n}{\Delta t}=\left.\dfrac{\partial u}{\partial t}\right|_j^n+\underbrace{\dfrac12\left.\dfrac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t}_{O(\Delta t)}+\dots$$

Разность по времени вперёд даёт первый порядок по $\Delta t$.

Производная по координате — центральная разность

Раскладываем оба соседних узла по $x$:

$$u_{j+1}^n=u_j^n+\left.\dfrac{\partial u}{\partial x}\right|_j^n h+\dfrac12\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n h^2+\dfrac16\left.\dfrac{\partial^3 u}{\partial x^3}\right|_j^n h^3+\dots$$ $$u_{j-1}^n=u_j^n-\left.\dfrac{\partial u}{\partial x}\right|_j^n h+\dfrac12\left.\dfrac{\partial^2 u}{\partial x^2}\right|_j^n h^2-\dfrac16\left.\dfrac{\partial^3 u}{\partial x^3}\right|_j^n h^3+\dots$$

При вычитании члены чётного порядка (в том числе $u_{xx}$) сокращаются:

$$u_{j+1}^n-u_{j-1}^n=2\left.\dfrac{\partial u}{\partial x}\right|_j^n h+\dfrac13\left.\dfrac{\partial^3 u}{\partial x^3}\right|_j^n h^3+\dots$$ $$\Rightarrow\quad 8\,\dfrac{u_{j+1}^n-u_{j-1}^n}{2h}=8\left.\dfrac{\partial u}{\partial x}\right|_j^n+\underbrace{\dfrac{4h^2}{3}\left.\dfrac{\partial^3 u}{\partial x^3}\right|_j^n}_{O(h^2)}+\dots$$

Симметрия центральной разности убивает член $O(h)$, поэтому она даёт второй порядок по $h$.

Ошибка аппроксимации и вывод

Складываем разложения и вычитаем дифференциальное уравнение (точная часть $\;u_t+8u_x-2t\;$ обращается в ноль на решении):

$$\psi=\dfrac{\Delta t}{2}\,\dfrac{\partial^2 u}{\partial t^2}+\dfrac{4h^2}{3}\,\dfrac{\partial^3 u}{\partial x^3}+\dots=O(\Delta t)+O(h^2).$$

Схема имеет первый порядок аппроксимации по времени и второй по координате: $\;\psi=O(\Delta t+h^2)$. (Проверить разложения можно интерактивно на странице «Инструменты», /taylor.html.)

Что требуется

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

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}-3\,\frac{u_{j+1}^{n+1}-u_j^{n+1}}{h}=\frac{(j-1)h}{n\,\Delta t+1},$$

аппроксимирующей уравнение переноса $\dfrac{\partial u}{\partial t}-3\dfrac{\partial u}{\partial x}=\dfrac{x}{t+1}$ в точке $(t^n,x_j)$.

Метод гармоник

Погрешность решения удовлетворяет той же (однородной) разностной схеме, что и сама сеточная функция. Свободный член $\dfrac{(j-1)h}{n\Delta t+1}$ — это дискретизация источника $\dfrac{x}{t+1}$; он не содержит искомую функцию $u$ и на устойчивость не влияет — отбрасываем его. Подставляем гармонику

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

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

Схема неявная: пространственная разность $\dfrac{u_{j+1}^{n+1}-u_j^{n+1}}{h}$ берётся на верхнем слое $n+1$, поэтому после подстановки она даст множитель $\lambda^{\,n+1}$.

Тип схемы

Запишем уравнение в виде $u_t+v\,u_x=f$. Здесь $\dfrac{\partial u}{\partial t}-3\dfrac{\partial u}{\partial x}=f$, то есть скорость переноса $v=-3<0$. Использована правая (вперёд) разность $\dfrac{u_{j+1}^{n+1}-u_j^{n+1}}{h}$ — это корректное направление «против потока» для $v<0$.

Подстановка гармоники

Подставляем $u_j^n=\lambda^{\,n}e^{i\alpha j}$ в однородную схему и делим всё на $\lambda^{\,n}e^{i\alpha j}$. Учтём, что верхний слой даёт множитель $\lambda$:

$$\frac{\lambda-1}{\Delta t}-3\,\frac{\lambda e^{i\alpha}-\lambda}{h}=0.$$

Вынесем $\lambda$ из второго слагаемого и умножим на $\Delta t$:

$$\lambda-1-3\frac{\Delta t}{h}\,\lambda\bigl(e^{i\alpha}-1\bigr)=0 \quad\Longrightarrow\quad \lambda\!\left(1-3\frac{\Delta t}{h}\bigl(e^{i\alpha}-1\bigr)\right)=1.$$

Выразим удобную для анализа неявных схем обратную величину:

$$\frac{1}{\lambda}=1-3\frac{\Delta t}{h}\bigl(e^{i\alpha}-1\bigr)=1+3\frac{\Delta t}{h}-3\frac{\Delta t}{h}\,e^{i\alpha}.$$

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

Для неявных схем условие $|\lambda|\le 1$ удобно переписать как $\left|\dfrac1\lambda\right|\ge 1$: обратные множителю величины должны лежать вне единичного круга (или на его границе). Обозначим

$$r=|v|\frac{\Delta t}{h}=3\frac{\Delta t}{h}>0,\qquad \frac{1}{\lambda}=(1+r)-r\,e^{i\alpha}.$$

При изменении $\alpha$ точка $\dfrac1\lambda$ описывает на комплексной плоскости окружность с центром в $(1+r,\,0)$ и радиусом $r$. Её ближайшая к началу координат точка имеет абсциссу

$$(1+r)-r=1.$$

Эквивалентно, прямой подсчёт модуля даёт

$$\left|\frac1\lambda\right|^2=1+2r(1+r)\bigl(1-\cos\alpha\bigr)\ge 1\qquad(r>0,\ 1-\cos\alpha\ge 0).$$

Значит, при любом $r>0$ вся окружность лежит вне (на границе — лишь при $\alpha=0$) единичного круга:

$$\left|\frac1\lambda\right|\ge 1\quad\Longrightarrow\quad |\lambda|\le 1\quad\text{при всех }\alpha.$$

(Знаменатель $\tfrac1\lambda=(1+r)-r e^{i\alpha}$ в нуль не обращается, так как $|(1+r)/r|>1$, поэтому $\lambda$ всегда конечен. Численная проверка: для $r=0{,}1;\ 1;\ 5;\ 100$ получается $\max_\alpha|\lambda|=1$.)

Вывод

Неравенство $|\lambda|\le 1$ выполнено при любых $\Delta t$ и $h$ — никакого ограничения на соотношение шагов не возникает. Неявная схема переноса с правой разностью (при $v=-3<0$, то есть в направлении против потока) абсолютно (безусловно) устойчива. Источник $\dfrac{(j-1)h}{n\Delta t+1}$ на устойчивость не влияет.

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

Для уравнения теплопроводности с источником

$$\dfrac{\partial u}{\partial t}=\dfrac{\partial^2 u}{\partial x^2}+x\,t,\qquad u=u(t,x),$$

с начальным условием $u(t{=}0,x)=0$ и краевыми условиями

$$u(t,x{=}0)=t \quad\text{(левое, 1-го рода — Дирихле)},$$ $$2\,\dfrac{\partial u}{\partial x}(t,x{=}1)=3\big[u(t,x{=}1)-10\big]\quad\text{(правое, 3-го рода — Робина)}.$$

Здесь коэффициент диффузии $\sigma=1$, конвективного члена нет, свободный член $f(t,x)=x\,t$ (зависит и от $t$, и от $x$). Сетка: по времени $t^{n}=n\,\Delta t$; по координате $x_j=(j-1)h$, $j=1,\dots,N$, $h=\dfrac{1}{N-1}$, $u_j^{n}=u(t^{n},x_j)$. Узел $j{=}1$ отвечает $x{=}0$, узел $j{=}N$ — $x{=}1$.

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

В неявной (чисто имплицитной) схеме производную по времени берут левой разностью, а все пространственные операторы — на новом слое $(n{+}1)$; свободный член также на слое $(n{+}1)$:

$$\dfrac{u_j^{n+1}-u_j^{n}}{\Delta t}=\dfrac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^{2}}+x_j\,t^{n+1}.$$

Схема абсолютно устойчива (безусловно по $\Delta t$ и $h$); порядок аппроксимации внутренней схемы $O(\Delta t,\,h^{2})$.

2. Приведение к трёхдиагональному виду

Домножаем на $\Delta t$ и переносим все неизвестные слоя $(n{+}1)$ налево, а известные (слой $n$ и источник) — направо. Записываем в форме

$$A_j\,u_{j-1}^{n+1}+C_j\,u_j^{n+1}+B_j\,u_{j+1}^{n+1}=F_j,$$

где (поскольку $\sigma=1$, конвекции нет)

$$A_j=-\dfrac{\Delta t}{h^{2}},\qquad C_j=1+\dfrac{2\Delta t}{h^{2}},\qquad B_j=-\dfrac{\Delta t}{h^{2}},$$ $$F_j=u_j^{n}+\Delta t\,x_j\,t^{n+1}.$$

(Конвективного члена $C\,\partial u/\partial x$ в уравнении нет, поэтому добавок $\mp C\Delta t/(2h)$ к коэффициентам при $u_{j\pm1}$ не возникает, и $A_j=B_j$.)

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

Достаточное условие — диагональное преобладание $|C_j|\ge|A_j|+|B_j|$:

$$|A_j|+|B_j|=\dfrac{2\Delta t}{h^{2}}\;<\;1+\dfrac{2\Delta t}{h^{2}}=|C_j|.$$

Неравенство строгое и выполняется при любых $\Delta t,h>0$ — преобладание обеспечивает единица в $C_j$, пришедшая из производной по времени. Прогонка устойчива (корректна) безусловно.

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

Решение ищем в виде $u_j^{n+1}=\alpha_j\,u_{j+1}^{n+1}+\beta_j$. Подставляя $u_{j-1}^{n+1}=\alpha_{j-1}u_j^{n+1}+\beta_{j-1}$ в разностное уравнение, получаем прямой ход прогонки ($j=2,\dots,N-1$):

$$\alpha_j=\dfrac{-B_j}{C_j+A_j\,\alpha_{j-1}},\qquad \beta_j=\dfrac{F_j-A_j\,\beta_{j-1}}{C_j+A_j\,\alpha_{j-1}}.$$

С учётом $A_j=B_j=-\dfrac{\Delta t}{h^2}$, $C_j=1+\dfrac{2\Delta t}{h^2}$:

$$\alpha_j=\dfrac{\dfrac{\Delta t}{h^2}}{\,1+\dfrac{2\Delta t}{h^2}-\dfrac{\Delta t}{h^2}\alpha_{j-1}\,},\qquad \beta_j=\dfrac{F_j+\dfrac{\Delta t}{h^2}\beta_{j-1}}{\,1+\dfrac{2\Delta t}{h^2}-\dfrac{\Delta t}{h^2}\alpha_{j-1}\,}.$$

5. Левое граничное условие — $\alpha_1,\beta_1$

Левое ГУ 1-го рода $u(t,x{=}0)=t$ задаёт значение в узле $j{=}1$ напрямую:

$$u_1^{n+1}=t^{n+1}.$$

Сравнивая с прогоночной формой $u_1^{n+1}=\alpha_1\,u_2^{n+1}+\beta_1$ (значение не зависит от $u_2^{n+1}$), получаем стартовые коэффициенты:

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

6. Правое граничное условие — $u_N^{n+1}$

Правое ГУ 3-го рода $2\,\dfrac{\partial u}{\partial x}(t,1)=3\big[u(t,1)-10\big]$ аппроксимируем левой (крайней) односторонней разностью на правой границе (точность по $x$ на границе — $O(h)$):

$$2\,\dfrac{u_N^{n+1}-u_{N-1}^{n+1}}{h}=3\big(u_N^{n+1}-10\big).$$

Домножая на $h$ и раскрывая: $2u_N^{n+1}-2u_{N-1}^{n+1}=3h\,u_N^{n+1}-30h$, откуда $2u_{N-1}^{n+1}=(2-3h)\,u_N^{n+1}+30h$. Подставляя обратный ход $u_{N-1}^{n+1}=\alpha_{N-1}u_N^{n+1}+\beta_{N-1}$:

$$2\big(\alpha_{N-1}u_N^{n+1}+\beta_{N-1}\big)=(2-3h)\,u_N^{n+1}+30h \;\Rightarrow\;\big(2\alpha_{N-1}-2+3h\big)u_N^{n+1}=30h-2\beta_{N-1},$$

и окончательно значение на правой границе:

$$\boxed{u_N^{n+1}=\dfrac{2\beta_{N-1}-30h}{\,2-2\alpha_{N-1}-3h\,}.}$$

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

1. Начальный слой $u_j^{0}=0$ (из начального условия).
2. Цикл по слоям $n=0,1,2,\dots$ — на каждом слое прогонка:
   • из левого ГУ: $\alpha_1=0,\ \beta_1=t^{n+1}$;
   • прямой ход $j=2,\dots,N-1$: вычислить $F_j=u_j^{n}+\Delta t\,x_j\,t^{n+1}$ и $\alpha_j,\beta_j$;
   • из правого ГУ 3-го рода: $u_N^{n+1}=\dfrac{2\beta_{N-1}-30h}{2-2\alpha_{N-1}-3h}$;
   • обратный ход $j=N-1,\dots,1$: $u_j^{n+1}=\alpha_j u_{j+1}^{n+1}+\beta_j$.
3. Перейти к следующему слою до нужного момента времени.

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

$$\dfrac{\partial u}{\partial t}+6\,\dfrac{\partial u}{\partial x}=0,\qquad u(t{=}0,x)=0,\qquad u(t,x{=}0)=e^{t}.$$

Это уравнение переноса (1-го порядка): вторых производных нет, диффузии и реакции нет, свободный член $f=0$. Конвективная скорость $v=6>0$. Сетка: $t^{n}=n\Delta t$, $x_j=(j-1)h$, $j=1,\dots,N_x$; узел $j=1$ отвечает $x=0$.

Выбор разности по знаку скорости. При $v>0$ поток входит слева, поэтому для $u_x$ берётся левая («против потока», upwind) разность $\dfrac{u_j-u_{j-1}}{h}$ и используется левое граничное условие $u(t,0)=e^{t}$. Правое условие для такой задачи не требуется (правая граница — выходная по потоку).

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

Производную по времени берём правой разностью, пространственный оператор $u_x$ — на новом $(n{+}1)$-м слое (отсюда «неявная»):

$$\frac{u_j^{n+1}-u_j^{n}}{\Delta t}+6\,\frac{u_j^{n+1}-u_{j-1}^{n+1}}{h}=0,\qquad j=2,\dots,N_x.$$

Начальное условие на сетке: $u_j^{0}=0$ для всех $j$.

2. Приведение к рекуррентному виду

В отличие от уравнений 2-го порядка (где неявная схема даёт трёхдиагональную систему и решается прогонкой), здесь на новом слое присутствуют лишь два неизвестных — $u_j^{n+1}$ и $u_{j-1}^{n+1}$. Значит, систему решать не нужно: схема приводится к прямому рекуррентному соотношению. Умножаем на $\Delta t$ и вводим $r=\dfrac{6\Delta t}{h}$:

$$u_j^{n+1}-u_j^{n}+r\big(u_j^{n+1}-u_{j-1}^{n+1}\big)=0\;\Rightarrow\;(1+r)\,u_j^{n+1}=u_j^{n}+r\,u_{j-1}^{n+1}.$$

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

Неизвестное узла выражается через значение на старом слое $u_j^{n}$ и через уже найденного «по потоку» левого соседа $u_{j-1}^{n+1}$:

$$\boxed{\,u_j^{n+1}=\frac{u_j^{n}+\dfrac{6\Delta t}{h}\,u_{j-1}^{n+1}}{1+\dfrac{6\Delta t}{h}}\,},\qquad j=2,\dots,N_x.$$

Соотношение замкнуто: при счёте слева направо величина $u_{j-1}^{n+1}$ к моменту вычисления $u_j^{n+1}$ уже известна (начиная с границы $u_1^{n+1}=e^{t^{n+1}}$). Поэтому полноценная прогонка (прямой/обратный ход с коэффициентами $\alpha_j,\beta_j$) не нужна — достаточно одного прохода по $j$.

4. Устойчивость

Неявная upwind-схема для уравнения переноса абсолютно устойчива. Подстановка гармоники $u_j^n=\lambda^n e^{\,i j\varphi}$ даёт множитель перехода $\lambda=\dfrac{1}{1+r-r e^{-i\varphi}}$, откуда $|\lambda|^2=\dfrac{1}{1+2r(r+1)(1-\cos\varphi)}\le1$ при любых $\varphi$ и любом $r=\dfrac{6\Delta t}{h}>0$. Знаменатель $1+\dfrac{6\Delta t}{h}>1$ при любых $\Delta t,h>0$, ограничения на шаг $\Delta t$ нет (в отличие от явной схемы, где требуется условие Куранта $\dfrac{6\Delta t}{h}\le1$).

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

1. Инициализация (НУ): $u_j^{0}=0$ для всех $j=1,\dots,N_x$.
2. Цикл по слоям $n=0,1,\dots,N_t-1$ (текущее время $t^{n+1}=(n{+}1)\Delta t$):
   • Левая граница (ГУ): $u_1^{n+1}=e^{t^{n+1}}=e^{(n+1)\Delta t}$.
   • Прямой проход слева направо $j=2,\dots,N_x$: по рекуррентной формуле $$u_j^{n+1}=\frac{u_j^{n}+\frac{6\Delta t}{h}\,u_{j-1}^{n+1}}{1+\frac{6\Delta t}{h}}$$ (на каждом шаге $u_{j-1}^{n+1}$ уже вычислено).
3. Переход к следующему слою: принять $u_j^{n}\leftarrow u_j^{n+1}$ и повторить, пока не достигнуто $t_k$.

Схема имеет порядок аппроксимации $O(\Delta t,h)$, абсолютно устойчива; решение находится прямым рекуррентным пересчётом за один проход по координате на каждом временном слое.

Что требуется

Определить тип уравнения (эллиптическое / параболическое / гиперболическое) и обосновать выбор. Уравнение:

$$\dfrac{\partial u}{\partial t} + 2\,\dfrac{\partial u}{\partial x} - 3\,\dfrac{\partial u}{\partial y} = 6\,\dfrac{\partial^2 u}{\partial x^2} + xy,\qquad u=u(t,x,y).$$

Как классифицируют (метод курса, вопрос 1, семинар 0)

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

$$a_{11}\dfrac{\partial^2 u}{\partial \xi_1^2}+2a_{12}\dfrac{\partial^2 u}{\partial \xi_1\partial \xi_2}+a_{22}\dfrac{\partial^2 u}{\partial \xi_2^2}+\dots=f,\qquad D=a_{12}^2-a_{11}a_{22}.$$

• $D<0$ — эллиптический тип (обе вторые производные одного знака, как у Лапласа);
• $D=0$ — параболический (одна из старших производных отсутствует, как у уравнения теплопроводности);
• $D>0$ — гиперболический (как у волнового уравнения).

Применяем к нашему уравнению

Независимых переменных три: $t,\,x,\,y$. Перенесём всё в одну часть и выпишем главную (старшую) часть:

$$6\,\dfrac{\partial^2 u}{\partial x^2}\;-\;\dfrac{\partial u}{\partial t}-2\,\dfrac{\partial u}{\partial x}+3\,\dfrac{\partial u}{\partial y}\;+\;xy=0.$$

Среди вторых производных присутствует только одна — $\dfrac{\partial^2 u}{\partial x^2}$ с коэффициентом $6$. Вторых производных по $t$, по $y$ и смешанных нет. Значит, коэффициенты старшей части:

$$a_{xx}=6,\qquad a_{tt}=0,\qquad a_{yy}=0,\qquad a_{tx}=a_{ty}=a_{xy}=0.$$

Производные $\dfrac{\partial u}{\partial t}$, $2\dfrac{\partial u}{\partial x}$, $-3\dfrac{\partial u}{\partial y}$ — первого порядка и в дискриминант не входят; слагаемое $xy$ — свободный член (источник), на тип также не влияет.

Матрица коэффициентов старших производных по $(t,x,y)$ диагональна: $A=\operatorname{diag}(0,6,0)$, поэтому $\det A=0$ и её собственные значения равны $0,\,0,\,6$ — одно положительное, два нулевых. Эквивалентно, по любой паре переменных, где могла бы появиться смешанная/вторая производная вместе с $x$, дискриминант равен нулю:

$$D=a_{tx}^2-a_{tt}\,a_{xx}=0^2-0\cdot 6=0,\qquad D=a_{xy}^2-a_{xx}\,a_{yy}=0^2-6\cdot 0=0.$$

Вывод

$D=0$ — уравнение параболического типа. Структурно это видно сразу: есть первая производная по времени $\dfrac{\partial u}{\partial t}$ и ровно одна вторая производная по пространственной координате $6\dfrac{\partial^2 u}{\partial x^2}$ — это классическая схема уравнения теплопроводности (диффузии). Дополнительные конвективные члены $2\dfrac{\partial u}{\partial x}-3\dfrac{\partial u}{\partial y}$ (перенос) и источник $xy$ остаются членами младшего порядка и тип не меняют: уравнение остаётся параболическим.

Для параболических уравнений курса далее применяют явные/неявные разностные схемы (см. семинар 4; страницу «Инструменты» — для проверки порядка аппроксимации).

Что требуется

Определить порядок аппроксимации разностной схемы в точке $(t^n,x_j)$:

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}-3\,\frac{u_{j+1}^n-u_j^n}{h}=n\,\Delta t-3(j-1)h,$$ $$\text{аппроксимирующей}\qquad \frac{\partial u}{\partial t}-3\,\frac{\partial u}{\partial x}=t-3x.$$

Метод

Подставляем в схему точное гладкое решение $u(t,x)$ и раскладываем сеточные значения в ряд Тейлора относительно узла $(t^n,x_j)$. Разность между разностным и дифференциальным операторами — это ошибка аппроксимации $\psi$; её главный член (с наименьшими степенями $\Delta t$ и $h$) задаёт порядок. Здесь $t^n=n\,\Delta t$, $x_j=(j-1)h$, поэтому правая часть $n\,\Delta t-3(j-1)h=t^n-3x_j$ совпадает с $t-3x$ в узле — источник аппроксимирован точно и в $\psi$ не входит.

Производная по времени — правая (вперёд) разность

$$u_j^{n+1}=u_j^n+\left.\frac{\partial u}{\partial t}\right|_j^n\Delta t+\frac12\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t^2+\dots$$ $$\Rightarrow\quad \frac{u_j^{n+1}-u_j^n}{\Delta t}=\left.\frac{\partial u}{\partial t}\right|_j^n+\underbrace{\frac12\left.\frac{\partial^2 u}{\partial t^2}\right|_j^n\Delta t}_{O(\Delta t)}+\dots$$

Правая разность по времени даёт первый порядок по $\Delta t$.

Производная по координате — правая (вперёд) разность

$$u_{j+1}^n=u_j^n+\left.\frac{\partial u}{\partial x}\right|_j^n h+\frac12\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n h^2+\dots$$ $$\Rightarrow\quad -3\,\frac{u_{j+1}^n-u_j^n}{h}=-3\left.\frac{\partial u}{\partial x}\right|_j^n-\underbrace{3\cdot\frac{h}{2}\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n}_{O(h)}-\dots$$

Правая (односторонняя) разность по координате даёт первый порядок по $h$. (Если бы стояла центральная разность $\dfrac{u_{j+1}^n-u_{j-1}^n}{2h}$, был бы $O(h^2)$, но здесь разность односторонняя.)

Ошибка аппроксимации

Складываем разложения временного и координатного слагаемых и вычитаем точно совпавшую правую часть. Дифференциальное уравнение $\left.\dfrac{\partial u}{\partial t}\right|_j^n-3\left.\dfrac{\partial u}{\partial x}\right|_j^n-(t^n-3x_j)=0$ обнуляет главные члены, и остаётся:

$$\psi=\frac{\Delta t}{2}\,\frac{\partial^2 u}{\partial t^2}\bigg|_j^n-\frac{3h}{2}\,\frac{\partial^2 u}{\partial x^2}\bigg|_j^n+\dots=O(\Delta t)+O(h).$$

Вывод

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

$$\psi=O(\Delta t)+O(h)=O(\Delta t+h).$$

Что требуется

Методом гармоник исследовать устойчивость явной разностной схемы

$$\dfrac{u_j^{n+1}-u_j^n}{\Delta t} - 3\,\dfrac{u_{j+1}^n - u_j^n}{h} = n\,\Delta t - 3(j-1)h,$$

аппроксимирующей уравнение переноса $\dfrac{\partial u}{\partial t} - 3\dfrac{\partial u}{\partial x} = t - 3x$ в точке $(t^n, x_j)$.

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

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

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

где $\lambda$ — множитель перехода за один шаг по времени, $\alpha$ — фаза. Необходимое условие устойчивости: $|\lambda|\le 1$ при всех $\alpha$. Схема явная — пространственная разность берётся на нижнем слое $n$.

Стандартный вид

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

$$\dfrac{u_j^{n+1}-u_j^n}{\Delta t} - 3\,\dfrac{u_{j+1}^n - u_j^n}{h} = 0\quad(\text{однородная схема}).$$

Это явная схема уравнения переноса $\dfrac{\partial u}{\partial t}+v\dfrac{\partial u}{\partial x}=f$ с аппроксимацией $\dfrac{\partial u}{\partial x}$ правой разностью: $-3\,(u_{j+1}^n-u_j^n)/h=+v\,(u_{j+1}^n-u_j^n)/h$ при $v=-3<0$. Свободный член $n\Delta t-3(j-1)h$ искомой функции не содержит — отброшен.

Подстановка гармоники

Подставляем $u_j^n=\lambda^n e^{i\alpha j}$ в однородную схему и делим на $\lambda^n e^{i\alpha j}$ (используем $u_{j+1}^n=\lambda^n e^{i\alpha(j+1)}=u_j^n\,e^{i\alpha}$):

$$\dfrac{\lambda-1}{\Delta t}-3\,\dfrac{e^{i\alpha}-1}{h}=0\quad\Rightarrow\quad \boxed{\;\lambda=1+3\,\dfrac{\Delta t}{h}\bigl(e^{i\alpha}-1\bigr).\;}$$

Анализ $|\lambda|\le 1$

Множитель $\lambda$ комплексный, поэтому $|\lambda|\le 1$ означает, что точки $\lambda(\alpha)$ лежат в единичном круге с центром в нуле. Обозначим

$$r=-v\,\dfrac{\Delta t}{h}=3\,\dfrac{\Delta t}{h}>0,\qquad \lambda=1-r+r\,e^{i\alpha}.$$

При изменении $\alpha$ от $0$ до $2\pi$ точка $\lambda$ описывает окружность с центром в $(1-r,\,0)$ и радиусом $r$ (так как $|r\,e^{i\alpha}|=r$). Эта окружность всегда проходит через точку $\lambda(0)=1$ (она лежит на единичной окружности). Граничной — определяющей выход за единичный круг — является диаметрально противоположная точка на вещественной оси, $\lambda(\pi)=1-2r$. Окружность целиком лежит в замкнутом единичном круге тогда и только тогда, когда эта точка не выходит наружу, то есть

$$|1-2r|\le 1\quad\Longleftrightarrow\quad -1\le 1-2r\le 1\quad\Longleftrightarrow\quad 0\le r\le 1.$$

(Это удобно подтвердить и точной оценкой модуля: $|\lambda|^2=(1-r+r\cos\alpha)^2+(r\sin\alpha)^2=1-4r(1-r)\sin^2\dfrac{\alpha}{2}$. При $r\le 1$ коэффициент $4r(1-r)\ge 0$, поэтому $|\lambda|^2\le 1$ при всех $\alpha$ — максимум $|\lambda|^2=1$ достигается на $\alpha=0$; при $r>1$ в худшем случае $\sin^2\frac{\alpha}{2}=1$ (то есть $\alpha=\pi$) получаем $|\lambda|^2=(1-2r)^2>1$, и схема неустойчива.)

Вывод — условие устойчивости

Условие $r\le 1$ даёт

$$3\,\dfrac{\Delta t}{h}\le 1\quad\Longleftrightarrow\quad \boxed{\;\Delta t\le \dfrac{h}{3}.\;}$$

Схема условно устойчива: шаг по времени ограничен (число Куранта $r=3\,\Delta t/h\le 1$). Это согласуется с правилом курса: при $v<0$ для уравнения переноса именно правая разность даёт условно устойчивую явную схему (условие $-1\le v\,\Delta t/h<0$, то есть $0<3\,\Delta t/h\le 1$). Если бы при $v=-3<0$ была взята левая разность, схема была бы абсолютно неустойчивой. Источник $t-3x$ на устойчивость не влияет.

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

$$\dfrac{d^2u}{dx^2}-11u=3x^3,\qquad u=u(x),$$ $$\begin{cases} 3\,\dfrac{du}{dx}(0)=u(0)-5,\\[2mm] 3\,\dfrac{du}{dx}(1)=2.\end{cases}$$

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

Это краевая задача для линейного ОДУ 2-го порядка (временного слоя нет, метод установления не требуется). Общий вид (гл. 10.1): $v\,\dfrac{du}{dx}=\sigma\dfrac{d^2u}{dx^2}-k\,u+f(x)$, $\sigma>0$. Заданное уравнение уже в этом виде:

$$0\cdot u' = 1\cdot u'' - 11\,u + f(x)\;\Rightarrow\; v=0,\quad \sigma=1,\quad k=+11,\quad f(x)=-3x^3.$$

Конвективный член отсутствует ($v=0$): первой производной в уравнении нет, выбор «левой/правой» разности для $u'$ не требуется. $k=11>0$ — слагаемое $-11u$ усиливает главную диагональ, поэтому достаточное условие сходимости прогонки выполнено и краевую задачу решаем методом прогонки напрямую.

2. Разностная схема ($O(h^2)$ внутри области)

Сетка $x_j=(j-1)h$, $j=1,\dots,N$, $h=\dfrac{1}{N-1}$. Вторую производную берём центральной разностью:

$$\dfrac{u_{j+1}-2u_j+u_{j-1}}{h^2}-11\,u_j=3x_j^3,\qquad j=2,\dots,N-1.$$

Порядок аппроксимации во внутренних узлах $O(h^2)$. Граничные условия ниже аппроксимируются односторонними разностями $O(h)$, поэтому глобальный порядок схемы — $O(h)$.

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

В форме $A_j u_{j-1}+C_j u_j+B_j u_{j+1}=F_j$:

$$\underbrace{\dfrac{1}{h^2}}_{A_j}u_{j-1}\;+\;\underbrace{\Big(-\dfrac{2}{h^2}-11\Big)}_{C_j}u_j\;+\;\underbrace{\dfrac{1}{h^2}}_{B_j}u_{j+1}\;=\;\underbrace{3x_j^3}_{F_j}.$$ $$A_j=B_j=\dfrac{1}{h^2},\qquad C_j=-\dfrac{2}{h^2}-11,\qquad F_j=3x_j^3.$$

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

Достаточное условие $|C_j|\ge|A_j|+|B_j|$:

$$|A_j|+|B_j|=\dfrac{2}{h^2}\;<\;\dfrac{2}{h^2}+11=\Big|-\dfrac{2}{h^2}-11\Big|=|C_j|.$$

Неравенство строгое при любом $h$ (запас $|C_j|-(|A_j|+|B_j|)=11>0$ даёт член $k=11>0$), поэтому прогонка устойчива и применима напрямую.

5. Прогоночное (рекуррентное) соотношение

Ищем решение в виде $u_j=\alpha_j\,u_{j+1}+\beta_j$. Подставляя $u_{j-1}=\alpha_{j-1}u_j+\beta_{j-1}$ в схему, получаем прямой ход ($j=2,\dots,N-1$):

$$\alpha_j=\dfrac{-B_j}{C_j+A_j\alpha_{j-1}},\qquad \beta_j=\dfrac{F_j-A_j\beta_{j-1}}{C_j+A_j\alpha_{j-1}}.$$

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

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

$$3\,\dfrac{u_2-u_1}{h}=u_1-5\;\Rightarrow\;\dfrac{3}{h}(u_2-u_1)=u_1-5\;\Rightarrow\;(h+3)\,u_1=3\,u_2+5h.$$

Откуда $u_1=\dfrac{3}{h+3}\,u_2+\dfrac{5h}{h+3}$, и, сравнивая с $u_1=\alpha_1 u_2+\beta_1$:

$$\boxed{\;\alpha_1=\dfrac{3}{h+3},\qquad \beta_1=\dfrac{5h}{h+3}\;}.$$

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

Правое условие $3\,\dfrac{du}{dx}(1)=2$, т.е. $\dfrac{du}{dx}(1)=\dfrac{2}{3}$, аппроксимируем левой разностью и подставляем обратный ход $u_{N-1}=\alpha_{N-1}u_N+\beta_{N-1}$:

$$3\,\dfrac{u_N-u_{N-1}}{h}=2\;\Rightarrow\;u_N-\big(\alpha_{N-1}u_N+\beta_{N-1}\big)=\dfrac{2h}{3}.$$ $$\boxed{\;u_N=\dfrac{\beta_{N-1}+\dfrac{2h}{3}}{1-\alpha_{N-1}}\;}.$$

8. Алгоритм

1. Задать сетку $x_j=(j-1)h$, $j=1,\dots,N$.
2. Из левого ГУ: $\alpha_1=\dfrac{3}{h+3}$, $\beta_1=\dfrac{5h}{h+3}$.
3. Прямой ход $j=2,\dots,N-1$: вычислить $A_j=B_j=\dfrac{1}{h^2}$, $C_j=-\dfrac{2}{h^2}-11$, $F_j=3x_j^3$ и $\alpha_j=\dfrac{-B_j}{C_j+A_j\alpha_{j-1}}$, $\beta_j=\dfrac{F_j-A_j\beta_{j-1}}{C_j+A_j\alpha_{j-1}}$.
4. Из правого ГУ: $u_N=\dfrac{\beta_{N-1}+\tfrac{2h}{3}}{1-\alpha_{N-1}}$.
5. Обратный ход $j=N-1,\dots,1$: $u_j=\alpha_j u_{j+1}+\beta_j$ — решение задачи.

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

$$\frac{\partial u}{\partial t}-0{,}9\,\frac{\partial u}{\partial x}=5\,\frac{\partial^2 u}{\partial x^2}-3u,\qquad u(t{=}0,x)=0,$$ $$u(t,x{=}0)=t,\qquad u(t,x{=}1)=2t-t^2,\qquad x\in[0,1].$$

Приведём к стандартному виду $\dfrac{\partial u}{\partial t}+v\dfrac{\partial u}{\partial x}=\sigma\dfrac{\partial^2 u}{\partial x^2}+f$. Здесь конвективный коэффициент (скорость переноса) $v=-0{,}9<0$, коэффициент диффузии $\sigma=5>0$, реакционный член $-3u$ (свободный член $f=0$). Граничные условия — 1-го рода (Дирихле). Сетка: $t^n=n\Delta t$, $\;x_j=(j-1)h$, $\;j=1,\dots,N_x$.

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

Все пространственные операторы берём на известном слое $n$. Производную по времени аппроксимируем правой разностью; вторую производную — центральной разностью на слое $n$. Конвективный член, так как $v=-0{,}9<0$, аппроксимируем правой («против потока») разностью — именно этот выбор даёт устойчивую (монотонную) явную схему при $v<0$. Реакционный член $-3u$ берём на слое $n$:

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}-0{,}9\,\frac{u_{j+1}^n-u_j^n}{h}=5\,\frac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}-3u_j^n.$$

Начальное условие: $u_j^0=0$ для всех узлов.

Граничные условия: $\;u_1^{n+1}=t^{n+1}=(n{+}1)\Delta t,\qquad u_{N_x}^{n+1}=2t^{n+1}-(t^{n+1})^2=2(n{+}1)\Delta t-\big((n{+}1)\Delta t\big)^2.$

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

В уравнении единственное неизвестное слоя $(n{+}1)$ — это $u_j^{n+1}$. Выражаем его:

$$u_j^{n+1}=u_j^n+0{,}9\,\frac{\Delta t}{h}\bigl(u_{j+1}^n-u_j^n\bigr)+5\,\frac{\Delta t}{h^2}\bigl(u_{j+1}^n-2u_j^n+u_{j-1}^n\bigr)-3\Delta t\,u_j^n.$$

Сгруппируем по узлам:

$$\boxed{\,u_j^{n+1}=\Big(\frac{5\Delta t}{h^2}+\frac{0{,}9\,\Delta t}{h}\Big)u_{j+1}^{n}+\Big(1-\frac{10\Delta t}{h^2}-\frac{0{,}9\,\Delta t}{h}-3\Delta t\Big)u_j^{n}+\frac{5\Delta t}{h^2}\,u_{j-1}^{n}\,}$$

Формула явная: правая часть содержит только значения известного слоя $n$ — новое значение в каждом внутреннем узле считается напрямую, без решения системы уравнений. (Отметим, что в отличие от случая $v>0$ здесь «лишний» множитель $\tfrac{0{,}9\Delta t}{h}$ добавляется к коэффициенту при $u_{j+1}^n$, а не при $u_{j-1}^n$, — следствие правой разностной аппроксимации сноса. При этом все три пространственных коэффициента неотрицательны при выполнении условия устойчивости, что обеспечивает монотонность схемы.)

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

Подстановка гармоники $u_j^n=\lambda^n e^{\,\mathrm{i}\,j\varphi}$ в рекуррентное соотношение даёт множитель перехода $\lambda=1-2\sigma\tfrac{\Delta t}{h^2}(1-\cos\varphi)-3\Delta t+\tfrac{0{,}9\Delta t}{h}(e^{\mathrm{i}\varphi}-1)$. Требование $|\lambda|\le1$ для всех $\varphi$ приводит (метод гармоник, гл. 6.2) к достаточному условию вида $|v|\dfrac{\Delta t}{h}+2\sigma\dfrac{\Delta t}{h^2}\le1$. Подставляя $|v|=0{,}9$, $\sigma=5$:

$$\boxed{\,0{,}9\,\frac{\Delta t}{h}+10\,\frac{\Delta t}{h^2}\le1\,}$$

откуда шаг по времени ограничен сверху; в грубой диффузионной оценке $\Delta t\lesssim h^2/(2\sigma)=h^2/10$.

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

1. Подготовка. Задать $h$ и $\Delta t$ (с учётом условия устойчивости $0{,}9\frac{\Delta t}{h}+10\frac{\Delta t}{h^2}\le1$), число узлов $N_x$ по координате и число слоёв $N_t$ по времени; сетка $x_j=(j-1)h$, $t^n=n\Delta t$.
2. Инициализация (начальный слой). $u_j^0=0$ для всех $j=1,\dots,N_x$ (начальное условие).
3. Цикл по слоям $n=0,1,\dots,N_t-1$:
   • внутренние узлы $j=2,\dots,N_x-1$ — по рекуррентной формуле (используются $u_{j-1}^n,u_j^n,u_{j+1}^n$);
   • левая граница из ГУ: $u_1^{n+1}=(n{+}1)\Delta t$;
   • правая граница из ГУ: $u_{N_x}^{n+1}=2(n{+}1)\Delta t-\big((n{+}1)\Delta t\big)^2$.
4. Контроль устойчивости. Шаг $\Delta t$ выбирается из условия п.3; результат — массив $u_j^n$ (сеточная функция приближённого решения).

Блок-схема алгоритма (явная разностная схема, параболическое уравнение) — на странице «Блок-схемы».