I. Функционал метода Бубнова–Галёркина

Запишем ещё раз функционал для уравнения Эйлера (5.4$\displaystyle I(\upsilon)$ $\displaystyle = (L[\upsilon], \upsilon)$ $\displaystyle - 2 \cdot (f, \upsilon).$)

Линейный оператор $L[\upsilon]$ для уравнения параболического типа, интересующего нас, имеет вид

Скалярное произведение $(L[\upsilon], \upsilon)$ имеет вид

Скалярное произведение $(f, \upsilon)$ имеет вид

Известна формула Остроградского для оператора Лапласа

где $S$ — граница подмногообразия $M$, а $dS$ — её векторный элемент, направленный по внешней нормали. Проинтегрируем по частям второй интеграл

Таким образом, исходный функционал можно записать в виде

Дадим приращение $\widehat{\upsilon}$ $+ \epsilon \cdot \upsilon$, где $\widehat{\upsilon}$ — точное решение, а $\epsilon$ — некоторое малое число. Полное приращение функционала можно записать в виде

Пренебрегая членами второго порядка малости, пропорциональными $\epsilon^2$, получим

Если вынести $\epsilon$ за скобки и принять во внимание, что $\epsilon$ может быть отрицательным, а вот приращение функционала всегда положительно, так как $\widehat{\upsilon}$ — точное решение, то получим в качестве обязательного условия

Учитывая, что $\widehat{\upsilon}$ $\approx \upsilon$, и сокращая общий множитель, мы возвращаемся к исходному уравнению. Это означает, что численное решение краевой задачи сводится к минимизации функционала

Множитель $2$ при линейном члене $\int_M f \cdot \upsilon dM$ сохраняем сознательно: именно он обеспечивает, что условие минимума $\delta I$ $= 0$ возвращает исходное уравнение, а не уравнение с лишним множителем $1/2$ (при дифференцировании квадратичные члены дают множитель $2$, и линейный член должен иметь такой же).

Перейдём к разделению на стационарную и нестационарную задачи. Принципиально важно проводить это разделение и последующую дискретизацию по времени на уровне уравнения, а не функционала: производная по времени не является самосопряжённым оператором, поэтому прямая подстановка $\upsilon$ $= \psi(M) \cdot \phi(t)$ в функционал привела бы к неверным коэффициентам. Запишем слабую форму уравнения $L[u]$ $= f$, полученную выше интегрированием по частям

которая должна выполняться для произвольной пробной функции $\upsilon$.

Стационарное уравнение с граничными условиями Дирихле и/или Неймана

Если $\frac{\displaystyle \partial u}{\displaystyle \partial t}$ $= 0$, временной член исчезает, и в силу симметрии пространственного оператора слабая форма (I.4$\displaystyle \int_M \frac{\displaystyle \partial u}{\displaystyle \partial t} \cdot \upsilon \,dM$ $\displaystyle + a^2 \cdot \int_M \nabla u \cdot \nabla \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla u \,dS$ $\displaystyle = \int_M f \cdot \upsilon \,dM,$) эквивалентна минимизации функционала

Нестационарное уравнение с граничными условиями Дирихле и/или Неймана

Производную по времени дискретизируем неявной схемой Эйлера с шагом $\Delta t$, применяя её к слабой форме уравнения: $\frac{\displaystyle \partial u}{\displaystyle \partial t}$ $\approx \frac{\displaystyle u_n - u_{n-1}}{\displaystyle \Delta t}$, где $u_n$ — решение на текущем временном слое, а $u_{n-1}$ — на предыдущем. Подставив это в (I.4$\displaystyle \int_M \frac{\displaystyle \partial u}{\displaystyle \partial t} \cdot \upsilon \,dM$ $\displaystyle + a^2 \cdot \int_M \nabla u \cdot \nabla \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla u \,dS$ $\displaystyle = \int_M f \cdot \upsilon \,dM,$) и домножив на $\Delta t$, получаем уравнение для шага по времени

Это уравнение, в свою очередь, является условием минимума функционала

Таким образом, нестационарное уравнение сведено к последовательному решению стационарных задач: на каждом шаге по известному $\upsilon_{n-1}$ находится $\upsilon_n$.

Матричная форма

Дискретизируем пробную функцию по узлам сетки: $\upsilon$ $= \sum_i q_i \cdot \phi_i(M)$, где $\phi_i(M)$ — базисные функции, а $\overrightarrow{q}$ $= (q_0, \ldots, q_n)^T$ — вектор узловых значений. Подстановка превращает каждый интеграл функционала в квадратичную или линейную форму по $\overrightarrow{q}$, а их коэффициенты собираются в матрицы:

где $K$ — матрица жёсткости (интеграл от произведения градиентов), $D$ — матрица демпфирования (интеграл от произведения пробных функций), $\overrightarrow{F}$ — вектор нагрузки (интеграл от источника); матрицы $K$ и $D$ симметричны. Граничный интеграл $-a^2 \int_S \upsilon \cdot \nabla \upsilon \,dS$ в матрицу не сворачивается: это поверхностный член, который в дискретной форме входит только в уравнения граничных узлов (у внутреннего узла пробная функция на границе равна нулю) и определяется самим граничным условием.

Дискретную систему получаем по методу Бубнова–Галёркина: подставляем $\upsilon$ $= \sum_i q_i \phi_i$ в слабую форму (I.4$\displaystyle \int_M \frac{\displaystyle \partial u}{\displaystyle \partial t} \cdot \upsilon \,dM$ $\displaystyle + a^2 \cdot \int_M \nabla u \cdot \nabla \upsilon \,dM$ $\displaystyle - a^2 \cdot \int_S \upsilon \cdot \nabla u \,dS$ $\displaystyle = \int_M f \cdot \upsilon \,dM,$) и поочерёдно берём пробные функции $\phi_j$. Для стационарного уравнения это даёт систему

где $\overrightarrow{s}$ — вектор узловых граничных потоков. Сам поток здесь ничем не задаётся — что с ним делать, определяет тип граничного условия. Учёт граничных условий — задача непростая, и подробно она разобрана в отдельных приложениях: «Учёт граничных условий Дирихле» и «Учёт граничных условий Неймана».

Для нестационарного уравнения (неявная схема Эйлера, см. (I.6$\displaystyle \int_M \upsilon_n^2 \,dM$ $\displaystyle + \Delta t \cdot a^2 \cdot \int_M \nabla \upsilon_n \cdot \nabla \upsilon_n \,dM$ $\displaystyle - 2 \cdot \int_M \left[ \Delta t \cdot f + \upsilon_{n-1} \right] \cdot \upsilon_n \,dM$ $\displaystyle - \Delta t \cdot a^2 \cdot \int_S \upsilon_n \cdot \nabla \upsilon_n \,dS$ $\displaystyle \rightarrow \min.$)) система на шаге по времени принимает вид