Неоднородные задачи: метод Фурье для уравнения с правой частью

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

Что ломается, когда у задачи появляется правая часть?#

Вспомним конструкцию из урока про метод разделения переменных: решение складывается из мод $T_n(t)X_n(x)$, каждая мода живёт независимо, со своим уравнением и своими данными. Вся машина стояла на одном фундаменте — однородности: убери нолик из уравнения или из граничного условия, и прямое разделение переменных рассыпается. Неоднородности, впрочем, живут не где попало, а в двух конкретных местах:

Начальные условия в таблице не упоминаются не случайно: они всегда неоднородны в каком-то смысле, но методу не мешают — просто перекрашивают коэффициенты. «Больными» считаются уравнение и границы. Пойдём по таблице сверху вниз, с полным числовым примером в каждой строке.

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

Как источник попадает в уравнения на моды?#

Схема первая, рабочая лошадка. Идея не отказываться от мод, а разложить по ним всё: и решение, и правую часть. Собственные функции $X_n$ берутся те же, что у однородной задачи, — источник не меняет спектр, он лишь раскачивает уже существующие моды. Для стержня с закреплёнными концами это синусы, и коэффициенты источника вынимаются ортогональностью — тем же механизмом, который доказывает задача Штурма–Лиувилля:

Подстановка в уравнение разваливает задачу на независимые обыкновенные дифференциальные уравнения — по одному на моду. Для тепла выходит уравнение первого порядка $T_n^{\prime }+a^2{\lambda }_n{T}_n=f_n(t)$, для волны — осциллятор $T_n^{\prime \prime }+{\omega }_n^2{T}_n=f_n(t)$. Техника решения — та же, что в курсе дифференциальных уравнений для линейных неоднородных ДУ: общее решение однородного плюс частное решение неоднородного; есть и второй взгляд — источник как поток мгновенных импульсов, который собирает принцип Дюамеля. Проведём её на конкретной задаче до конца.

1. 1
   Параметры: $L=1$, $a^2=1$. Собственные функции — $\sin (\pi nx)$, собственные значения ${\lambda }_n={\pi }^2{n}^2$.
2. 2
   Правая часть уже готова: $8\sin (\pi x)$ — ровно первая мода. Контроль по формуле: $f_1=2{\int }_0^{1}8{\sin }^2(\pi x)dx=8$, остальные $f_n=0$.
3. 3
   ОДУ первой моды: $T_1^{\prime }+{\pi }^2{T}_1=8$, начальное $T_1(0)=0$. Частное решение ищем константой $T_{\ast }=8/{\pi }^2$ — подстановка даёт ${\pi }^2{T}_{\ast }=8$.
4. 4
   Общее решение $T_1=8/{\pi }^2+C{e}^{-{\pi }^{2}t}$; начальное условие отбирает $C=-8/{\pi }^2$. Итог: $T_1(t)=\frac{8}{{\pi }^2}\left(1-e^{-{\pi }^{2}t}\right)$.
5. 5
   Моды $n\ge 2$ спят: их правые части нулевые, начальные условия тоже. Ответ: $u(x,t)=\frac{8}{{\pi }^2}\left(1-e^{-{\pi }^{2}t}\right)\sin (\pi x)$.
6. 6
   Три проверки. Границы: синус нулевой на концах. Старт: скобка при $t=0$ занулилась. Уравнение: $u_t-u_{xx}=8e^{-{\pi }^{2}t}\sin (\pi x)+8(1-e^{-{\pi }^{2}t})\sin (\pi x)=8\sin (\pi x)$ — сошлось.

Цифры. Константа $8/{\pi }^2\approx 0,8106$ — стационарный потолок моды. Значения по времени: $T_1(0,1)\approx 0,508$, $T_1(0,25)\approx 0,742$, $T_1(0,5)\approx 0,805$ — уже $99\%$ потолка, потому что показатель ${\pi }^2\approx 9,87$ работает быстро. Смысл предельного профиля: $\frac{8}{{\pi }^2}\sin (\pi x)$ — стационарное решение уравнения Пуассона $-u_{xx}=8\sin (\pi x)$ с нулями на концах. Источник работает постоянно, отток тепла через торцы растёт вместе с температурой, и на их равновесии рост останавливается.

Заметьте, что произошло с режимом затухания из урока про теплопроводность: экспонента $e^{-{\pi }^{2}t}$ никуда не делась, она просто переехала в скобку $(1-e^{-{\pi }^{2}t})$. Источник не отменил релаксацию — он добавил свой вклад поверх неё. Если бы начальная температура была не нулевой, к решению прибавились бы свободные экспоненты с коэффициентами из её разложения, и снова без всяких новых идей.

Почему при резонансе амплитуда растёт именно линейно?#

Вернёмся к струне. Мода с частотой ${\omega }_n$ качается, и на неё действует внешняя сила. Пусть сила повторяет форму моды и бьёт точно в её частоту: $T_n^{\prime \prime }+{\omega }_n^2{T}_n=F\cos ({\omega }_{n}t)$. Первый рефлекс — искать частное решение в виде $A\cos ({\omega }_{n}t)$ — проваливается молча: подстановка даёт слева ноль, потому что $\cos ({\omega }_{n}t)$ — решение однородного уравнения, а справа стоит $F$. Уравнение требует другой формы ответа.

Правильная проба — $T=Ct\sin ({\omega }_{n}t)$. Дифференцируем дважды: $T^{\prime }=C\sin ({\omega }_{n}t)+C{\omega }_{n}t\cos ({\omega }_{n}t)$, затем $T^{\prime \prime }=2C{\omega }_n\cos ({\omega }_{n}t)-C{\omega }_n^{2}t\sin ({\omega }_{n}t)$. Подставляем в уравнение: слагаемые с $t\sin$ гасятся против ${\omega }_n^{2}T$, и остаётся $2C{\omega }_n\cos ({\omega }_{n}t)=F\cos ({\omega }_{n}t)$. Отсюда $C=\frac{F}{2{\omega }_n}$ — и частное решение выписано окончательно:

Вот и линейный рост: периодический синус, умноженный на $t$. Физика подкладки прозрачна: каждые полпериода сила работает в такт движению и подкачивает энергию, как качели, толкаемые вовремя. Экспоненты нет и не будет: экспоненциальный рост — признак неустойчивости самой системы, а чистый осциллятор нейтрален. Сила, качающая в такт, даёт неубывающую подкачку — выходит прямая. Сравните с нерезонансной силой: уравнение $T_n^{\prime \prime }+{\omega }_n^2{T}_n=F\cos (\Omega t)$ при $\Omega \ne {\omega }_n$ имеет ограниченное частное решение с амплитудой $F/|{\omega }_n^2-{\Omega }^2|$; для $F=4$, ${\omega }_n=\pi$, $\Omega =2\pi$ это $4/(3{\pi }^2)\approx 0,135$. А резонансная раскачка при том же $F$ за десять секунд проходит $\frac{4}{2\pi }\cdot 10\approx 6,4$ — разница в десятки раз и она растёт со временем без предела.

Ту же прямую видно из энергетики. Мощность силы равна произведению силы на скорость моды; скорость резонирующей моды пропорциональна $t$, значит, и мощность подкачки растёт по времени. Энергия при этом накапливается как $t^2$, а амплитуда — корень квадратный из энергии: снова линейный рост. Никакого экспоненциального разгона в нейтральной системе взять неоткуда — подкачка лишь равномерно, такт в такт, вливает энергию.

Как сдвинуть задачу, чтобы границы стали нулевыми?#

Третья схема закрывает границы. Стержень $0<x<1$: левый конец держат при 100°, правый при 0°, начальная температура нулевая. Метод Фурье требует нулей на концах (это однородные граничные условия) — а их нет, и разложение по синусам не поможет: синусы у торцов нулевые, температуру 100 на конце они не соберут. Приём называется сдвигом: представим решение суммой $u=v+w$, где $v$ — простая функция, забирающая на себя граничные условия, а $w$ — остаток с нулевыми границами.

В роли $v$ напрашивается линейный профиль $v(x)=100(1-x)$ — он же стационарное решение этой задачи: вторая производная нулевая, концы согласованы — выполнено условие совместности данных, и ряд не тянет лишних гармоник. Тогда $w=u-v$ имеет на концах нули (из $100-100$ и $0-0$), по времени ведёт себя как $u$, а по $x$ его вторая производная равна $u_{xx}$, потому что $v_{xx}=0$. Уравнение для $w$ осталось однородным, зато начальное условие перекрасилось: $w(x,0)=0-100(1-x)=100(x-1)$ — отрицательное всюду, стержень стартует ниже своего стационара. Дальше — схема первая, уже знакомая.

Выбор $v$ не уникален — подойдут любые функции с нужными границами. Стационарная предпочтительна по двум причинам: уравнение для $w$ остаётся чистым, без добавок, а предел $t\to \infty$ виден сразу — ряд для $w$ гаснет, и остаётся ровно $v$. Поэтому при постоянных температурах на торцах никто не ищет ничего хитрее линейного профиля.

1. 1
   Сдвиг: $u=v+w$ с $v(x)=100(1-x)$; для $w$ границы нулевые, уравнение $w_t=w_{xx}$ прежнее.
2. 2
   Начальные данные: $w(x,0)=100(x-1)$. Раскладываем по синусам: $b_n=2{\int }_0^{1}100(x-1)\sin (\pi nx)dx$.
3. 3
   Интегрирование по частям даёт $b_n=-\frac{200}{\pi n}$ для всех $n$ — не нулевые ни при чётных, ни при нечётных: функция всюду отрицательна.
4. 4
   Ряд: $w(x,t)=-\frac{200}{\pi }{\sum }_{n\ge 1}\frac{1}{n}{e}^{-{\pi }^2{n}^{2}t}\sin (\pi nx)$.
5. 5
   Ответ: $u=100(1-x)-\frac{200}{\pi }\sum \frac{1}{n}{e}^{-{\pi }^2{n}^{2}t}\sin (\pi nx)$. Проверка при $t=0$: ряд $\sum \frac{\sin (\pi nx)}{n}=\frac{\pi (1-x)}{2}$ на $(0,1)$, и $u(x,0)=0$ — стартовый ноль восстановлен.
6. 6
   Числа: $u(0,5;0,1)\approx 26,3^{\circ }$ — середина ещё холодная; к $t=3$ остаётся чистый стационар $100(1-x)$, в середине $50^{\circ }$.

Логика универсальна: сдвиг забирает неоднородность границ, а что остаётся — раскладывается по модам. Если граничные условия зависят от времени, стационарной $v$ не существует — тогда берут любую удобную функцию с нужными границами (обычно полином по $x$, коэффициенты которого подогнаны под условия), и уравнение для $w$ получает добавку $-v_t$: схема три превращается в схему одну. Цепочка приёмов короткая, и замыкается она в самую первую схему.

Что делать с каждым типом неоднородности?#

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

Общая мораль читается из таблицы: спектр однородной задачи не зависит от источников и печек, поэтому неоднородность не ломает метод, а дополняет его — той же ортогональностью, теми же модами. Дальше по курсу: переменные коэффициенты и обобщённые спектры — в уроке про Штурма–Лиувилля, бесконечные области и сплошной спектр — в преобразовании Фурье для УЧП, а стационарные задачи с ненулевой границей без времени — в задаче Дирихле для круга. Формулы темы лежат в шпаргалке по УЧП.