Контрольная работа № 1
подробные решения примерного варианта · темы семинаров 0–5

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

Определить тип уравнения и обосновать выбор. Уравнение:

$$0{,}2\,\frac{\partial u}{\partial t}=5\,\frac{\partial^2 u}{\partial x^2}+2x^2,\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$ — гиперболический (как у волнового уравнения).

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

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

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

Производная по времени $0{,}2\,\partial u/\partial t$ — первого порядка, в дискриминант старших производных не входит. Тогда

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

Вывод

$D=0$ — уравнение параболического типа. Это и видно структурно: первая производная по времени плюс вторая по координате — классическое уравнение теплопроводности (диффузии) с источником $2x^2$. Свободный член $2x^2$ и коэффициенты на тип не влияют — он определяется только старшими производными.

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

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

Определить порядок аппроксимации разностной схемы в точке $(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},$$ $$\text{аппроксимирующей}\qquad \frac{\partial u}{\partial t}+8\,\frac{\partial u}{\partial x}=e^{tx}.$$

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

Подставляем в схему точное решение $u(t,x)$ и раскладываем сеточные значения в ряд Тейлора относительно узла $(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+\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 8\,\frac{u_j^n-u_{j-1}^n}{h}=8\left.\frac{\partial u}{\partial x}\right|_j^n-\underbrace{8\cdot\frac{h}{2}\left.\frac{\partial^2 u}{\partial x^2}\right|_j^n}_{O(h)}+\dots$$

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

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

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

Схема имеет первый порядок аппроксимации по времени и первый по координате: $\;\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.$$

Метод гармоник (вопрос 7, семинар 3)

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

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

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

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

Поделим на $0{,}2$: $\;\dfrac{u_j^{n+1}-u_j^n}{\Delta t}=\sigma\,\dfrac{u_{j+1}^n-2u_j^n+u_{j-1}^n}{h^2}$, где эффективный коэффициент

$$\sigma=\frac{5}{0{,}2}=25.$$

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

Подставляем $u_j^n=\lambda^n e^{i\alpha j}$ и делим на $\lambda^n e^{i\alpha j}$. Используем $e^{i\alpha}-2+e^{-i\alpha}=2\cos\alpha-2=-4\sin^2\dfrac{\alpha}{2}$:

$$\frac{\lambda-1}{\Delta t}=\sigma\,\frac{e^{i\alpha}-2+e^{-i\alpha}}{h^2}=-\,\frac{4\sigma}{h^2}\sin^2\frac{\alpha}{2}.$$ $$\boxed{\;\lambda=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\le 1$ всегда. Жёсткое условие — снизу: $\lambda\ge -1$. Наихудший случай $\sin^2\frac{\alpha}{2}=1$:

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

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

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

Явная схема условно устойчива: шаг по времени жёстко ограничен квадратом шага по координате. (Проверить условие на числах можно на странице «Инструменты».) Для снятия ограничения используют неявные схемы (см. КР1, задача 4).

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

$$\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=-5t$; оба граничных условия — 2-го рода (Неймана). Сетка: $t^n=n\Delta t$, $x_j=(j-1)h$, $j=1\dots N_x$, $u_j^n=u(t^n,x_j)$.

1. Схема Кранка-Николсона (вопрос 10, семинар 5)

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

$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=\frac{7}{2}\,\frac{u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{h^2}+\frac{7}{2}\,\frac{u_{j+1}^{n}-2u_j^{n}+u_{j-1}^{n}}{h^2}-5\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{7\Delta t}{2h^2},\qquad b_j=1+\frac{7\Delta t}{h^2},$$ $$\xi_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\,\Delta t.$$

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

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

$$|a_j|+|c_j|=\frac{7\Delta t}{h^2}<1+\frac{7\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}$, получаем прямой ход:

$$\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)=0$ аппроксимируем правой разностью на левой границе: $\dfrac{u_1^{n+1}-u_0^{n+1}}{h}=0\Rightarrow u_1^{n+1}=u_2^{n+1}$ (нумерация от $j{=}1$; левый узел выражается через соседний). В форме прогонки это даёт стартовые

$$\alpha_1=1,\qquad \beta_1=0.$$

(Эквивалентно: фиктивный узел даёт $\alpha_0=1,\ \beta_0=0$, откуда $\alpha_1=\dfrac{-a_1}{b_1+c_1},\ \beta_1=\dfrac{\xi_1}{b_1+c_1}$ — обе записи дают $u_1=u_2$.)

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

Условие $\dfrac{\partial u}{\partial x}(t,1)=0$ аппроксимируем левой разностью: $\dfrac{u_{N_x}^{n+1}-u_{N_x-1}^{n+1}}{h}=0\Rightarrow u_{N_x}^{n+1}=u_{N_x-1}^{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}=\frac{\beta_{N_x-1}}{1-\alpha_{N_x-1}}.$$

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

7. Алгоритм

Цикл по слоям $n$: вычислить $\xi_j$ → прямой ход ($\alpha_j,\beta_j$ от $j{=}1$) → $u_{N_x}^{n+1}$ из правого ГУ → обратный ход. Блок-схема — ниже и на странице «Блок-схемы» (неявная схема, метод прогонки).

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

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

Конвективный коэффициент (скорость) $v=0{,}8>0$, диффузия $\sigma=3$, реакция $-2u$. Граничные условия — 1-го рода (Дирихле). Сетка: $t^n=n\Delta t$, $x_j=(j-1)h$.

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

Все пространственные операторы берём на известном слое $n$. Производная по времени — правая разность. Вторая производная — центральная разность. Конвективный член при $v>0$ — левая («против потока») разность; реакционный член $-2u$ — на слое $n$:

$$\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.$$

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

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

$$u_j^{n+1}=\frac{3\Delta t}{h^2}\,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}+\Big(\frac{3\Delta t}{h^2}+\frac{0{,}8\,\Delta t}{h}\Big)u_{j-1}^{n}.$$

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

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

1. Инициализация: $u_j^0=0$ для всех внутренних узлов (начальное условие).
2. Цикл по слоям $n=0,1,2,\dots$:
   • границы из ГУ: $u_1^{n+1}=t^{n+1}=(n{+}1)\Delta t$, $\;u_{N_x}^{n+1}=(t^{n+1})^2=\big((n{+}1)\Delta t\big)^2$;
   • внутренние узлы $j=2\dots N_x-1$ — по рекуррентной формуле выше.
3. Контроль устойчивости (явная схема условно устойчива): шаг $\Delta t$ ограничен сверху $\sim h^2/(2\sigma)$ — см. метод гармоник («Инструменты»).

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