Метод Галёркина: проекционные методы для краевых задач

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

Идея проекционных методов: невязка, которая обязана быть ортогональной#

Пусть краевая задача записана операторно: $Lu=f$ с граничными условиями. Точное решение мы не знаем, но умеем вычислять невязку: подставим в левую часть приближение $u_n$ и посмотрим, сколько не хватает до правой. Формально $R(x)=L{u}_n-f$ — функция, которая равна нулю только для точного решения. У приближения с конечным числом параметров невязка в ноль не обратится никогда, и это честно: вопрос лишь в том, какую ошибку считать допустимой и куда её «спрятать».

Проекционная схема отвечает так: выберем конечный набор базисных функций ${\phi }_1,\dots ,{\phi }_n$, каждая из которых удовлетворяет граничным условиям, и будем искать приближение линейной комбинацией $u_n=a_1{\phi }_1+\cdots +a_n{\phi }_n$. Коэффициенты подберём из требования: невязка ортогональна каждой базисной функции, то есть $⟨R,{\phi }_j⟩=0$ при всех $j$. Здесь $⟨g,h⟩$ — скалярное произведение функций, чаще всего интеграл $\int ghdx$ по области; его свойства мы разбирали в уроке про скалярное произведение. Смысл условия прозрачен: у метода нет информации о невязке вне натянутого базисом пространства, поэтому внутри базиса она обнуляется целиком, а наружу остаётся лишь непроектируемый остаток.

Получается $n$ уравнений для $n$ неизвестных коэффициентов — после подстановки комбинации в $⟨R,{\phi }_j⟩=0$ возникает система линейных алгебраических уравнений. Как её решать на компьютере и почему для больших $n$ важна структура матрицы — отдельная история численных методов для СЛАУ; здесь же мы всё проделаем руками на задаче, где каждый интеграл берётся.

Слабая форма: зачем понижать порядок производной#

Для задач математической физики ортогональность формулируют не к сильной форме уравнения, а к слабой. Возьмём модельный стержень: $-u^{\prime \prime }=f$ на отрезке $(0,1)$ с заделкой обоих концов $u(0)=u(1)=0$. Умножим уравнение на пробную функцию $v$ с нулями на концах и проинтегрируем по отрезку, а слева интегрируем по частям:

Что изменилось? Во-первых, старшая производная понизилась с $u^{\prime \prime }$ до $u^{\prime }$: теперь достаточно, чтобы решение имело квадратично интегрируемую первую производную, а не вторую. Во-вторых, граничные условия разделились на роль пробных функций — они вшиты в выбор $v$ и больше не отдельный пункт. Функция, для которой подобная форма имеет смысл, называется пробной (trial), а $v$ — весовой (test); в схеме Бубнова–Галёркина обе роли играют одни и те же базисные функции, и эта симметрия отвечает за главные достоинства метода.

Подставим в слабую форму комбинацию $u_n=\sum a_i{\phi }_i$ и последовательно возьмём $v={\phi }_1,\dots ,{\phi }_n$. Интегралы с производными дают матрицу жёсткости $K_{ij}=\int {\phi }_i^{\prime }{\phi }_j^{\prime }dx$, интегралы с правой частью — вектор нагрузок $F_j=\int f{\phi }_{j}dx$. Метод свёлся к системе:

Матрица $K$ симметрична и положительно определена — прямое следствие слабой формы. Именно эту матрицу в методе конечных элементов собирают по кусочкам-элементам, и именно её решают в каждом расчёте конструкции. Для нашей модельной задачи она считается до конца.

Ручной пример: стержень под равномерной нагрузкой#

Задача: $-u^{\prime \prime }=1$ на $(0,1)$, $u(0)=u(1)=0$. Физически это упругая нить или стержень под собственным весом: равномерная нагрузка, оба конца закреплены. Точное решение легко угадать дважды проинтегрировав: $u=x(1-x)/2$ — парабола с максимумом $1/8=0,125$ в середине. Оно нужно нам как эталон для проверки; сам метод его не использует.

Базис возьмём синусный: ${\phi }_1=\sin (\pi x)$, ${\phi }_2=\sin (2\pi x)$ — те самые собственные функции оператора $-d^2/d{x}^2$ из урока про метод разделения переменных и задачу Штурма–Лиувилля. Каждый синус нулевой на концах, значит, граничные условия выполнены при любых коэффициентах. Ищем $u_2=a_1{\phi }_1+a_2{\phi }_2$.

Считаем матрицу жёсткости. Производные: ${\phi }_1^{\prime }=\pi \cos (\pi x)$, ${\phi }_2^{\prime }=2\pi \cos (2\pi x)$. Тогда $K_{11}={\int }_0^1{\pi }^2{\cos }^2(\pi x)dx={\pi }^2/2$, потому что среднее значение квадрата косинуса по целому числу полупериодов равно $1/2$. Аналогично $K_{22}={\int }_0^{1}4{\pi }^2{\cos }^2(2\pi x)dx=2{\pi }^2$. Смешанные интегралы гасятся ортогональностью: $K_{12}=2{\pi }^2{\int }_0^1\cos (\pi x)\cos (2\pi x)dx=0$ — косинусы с разными частотами на $[0,1]$ ортогональны, поэтому матрица диагональна. Правая часть: $F_1={\int }_0^1\sin (\pi x)dx=[-\cos (\pi x)/\pi {]}_0^1=2/\pi$, а $F_2={\int }_0^1\sin (2\pi x)dx=0$.

Ответ одного приближения: $u_1(x)=\frac{4}{{\pi }^3}\sin (\pi x)\approx 0,1290\sin (\pi x)$. В середине стержня $u_1(0,5)=0,1290$ против точных $0,125$: завышение на $0,0040$, или $3,2\%$. Уже одной базисной функцией метод попал в ответ с точностью до трёх процентов — и это не везение, а свойство задачи, к которому мы вернёмся при разговоре о скорости сходимости.

Заметим тонкость, на которой спотыкаются на экзаменах. Невязка сильной формы здесь равна $R(x)=-u_1^{\prime \prime }-1=(4/\pi )\sin (\pi x)-1$: она нигде не ноль, в середине отрезка она равна $(4/\pi )-1\approx 0,273$. Требование «невязка мала поточечно» метод Галёркина не ставит — он требует ортогональности, и проверка $\int R{\phi }_{1}dx=0$ это подтверждает. Для самосопряжённого оператора с базисом из однородных граничных условий обе формулировки — ортогональность невязки сильной формы и слабая форма — эквивалентны: связывает их именно интегрирование по частям, то же, что в краевых задачах для ОДУ.

Теперь нарушим симметрию: возьмём ту же задачу, но с линейно растущей нагрузкой $f=x$. Вторая мода должна ожить. Считаем нагрузки: $F_1={\int }_0^{1}x\sin (\pi x)dx=1/\pi$ (интегрированием по частям: $-x\cos (\pi x)/\pi +\sin (\pi x)/{\pi }^2$, на концах остаётся $1/\pi$). Для второй моды $F_2={\int }_0^{1}x\sin (2\pi x)dx=-1/(2\pi )$ — знак минус: у правого, тяжёлого конца вторая мода имеет противоположный знак. Коэффициенты: $a_1=F_1/K_{11}=(1/\pi )\cdot (2/{\pi }^2)=2/{\pi }^3\approx 0,0645$, $a_2=F_2/K_{22}=(-1/(2\pi ))/(2{\pi }^2)=-1/(4{\pi }^3)\approx -0,00806$.

Точное решение при $f=x$: $u=(x-x^3)/6$ с максимумом $1/(9\sqrt{3})\approx 0,0642$ в точке $x=1/\sqrt{3}\approx 0,577$. В середине отрезка точное значение $0,0625$, двухчленное приближение даёт $u_2(0,5)=0,0645$ — погрешность снова около $3,2\%$. И ещё одна проверка, самая показательная: коэффициенты $2/{\pi }^3$ и $-1/(4{\pi }^3)$ в точности совпадают с первыми коэффициентами ряда Фурье точного решения по синусам. Для самосопряжённого оператора с собственным базисом Галёркин воспроизводит точное разложение, обрезая его на выбранном номере, — проекционный метод с таким базисом превращается в спектральный.

С какой скоростью сходится проекционное приближение?#

Продолжим таблицу коэффициентов для равномерной нагрузки. Для нечётных $n$ нагрузка даёт $F_n=2/(n\pi )$, а диагональный элемент $K_{nn}=(n\pi {)}^2/2$, откуда $a_n=4/(n^3{\pi }^3)$: коэффициенты падают как куб номера. Чётные моды молчат — нагрузка симметрична. Значит, приближение с $N$ синусами — это частичная сумма знакочередующегося ряда в середине отрезка, и ошибка там ограничена первым отброшенным коэффициентом.

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

Ритц, конечные элементы и гидродинамика: где это работает#

У схемы Бубнова–Галёркина есть двойник. Для самосопряжённых положительных операторов слабая форма совпадает с условием минимума энергетического функционала $J[v]=\frac{1}{2}\int (v^{\prime }{)}^{2}dx-\int fvdx$: продифференцируйте $J$ по коэффициентам комбинации — получите ровно систему $Ka=F$. Это метод Ритца: минимизация энергии на натянутом базисом подпространстве. Для симметричных задач Ритц и Галёркин неразличимы; для несамосопряжённых (конвекция, следящие нагрузки) энергетический функционал не существует, а Галёркин продолжает работать — в этом его преимущество. Если хочется увидеть вариационный взгляд на стационарные задачи со стороны классической постановки, он подробно разобран в уроке про уравнение Лапласа и задачу Дирихле для круга.

Второй потомок — метод конечных элементов. МКЭ — это метод Галёркина, в котором базис собран из кусочно-полиномиальных «шляпок»: каждая равна единице в своём узле и нуля во всех остальных. Матрица жёсткости при этом становится разреженной (каждая шляпка пересекается лишь с соседними), а точность регулируется размером сетки. Вся строительная механика — балки, плиты, оболочки, грунтовые массивы — решается сборкой и факторизацией таких матриц; инженеру важно понимать, что за интерфейсом «сетка и закрепления» стоит та самая система $Ka=F$, которую мы только что получили руками. В вычислительной гидродинамике методами взвешенных невязок дискретизируют уравнения Навье–Стокса, а для конвективно-доминируемых течений весовые функции специально сдвигают вверх по потоку (схема SUPG) — иначе осциллируют фронты.

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

Частые ошибки#

• Базис без граничных условий. Если ни одна базисная функция не удовлетворяет однородным краевым условиям, слабая форма теряет смысл: краевой член при интегрировании по частям не выпадает, и система описывает другую задачу. Неоднородные условия сначала вычищают сдвигом — так же, как в неоднородных задачах метода Фурье.
• Сильная форма вместо слабой. Требовать $⟨R,{\phi }_j⟩=0$ с вторыми производными от приближения — значит завысить гладкость и, что хуже, получить несимметричную систему. Интегрирование по частям — не украшение, а несущая конструкция метода.
• Линейная независимость. Повторяющиеся или почти совпадающие базисные функции делают матрицу $K$ плохо обусловленной: числа растут, ответ плывёт. Базис обязан быть линейно независимым — это базовое требование к любому базису.
• Симметрия как догма. Нулевой коэффициент при симметричной нагрузке — следствие чётности, а не общее правило. Стоит нагрузке стать несимметричной (наш пример $f=x$), чётные моды оживают, и точка максимума смещается с середины отрезка.
• Потерянный масштаб невязки. Ортогональность гарантирует наилучшее приближение в энергетической норме, а не поточечную точность. Контроль качества — сходимость по $n$ и сравнение с известными решениями, а не «маленькая невязка в одной точке».

Сводка алгоритма. Выбрать базис с нужными краевыми условиями и хорошей обусловленностью. Выписать слабую форму и собрать $K_{ij}$ и $F_j$ — аналитически или квадратурой. Решить систему, для больших $n$ — методом Гаусса с учётом разреженности. Проверить ортогональность невязки и сходимость при увеличении базиса. Суперпозиция мод, на которой всё держится, — тот же принцип суперпозиции линейных уравнений, что и в методе Фурье; различие лишь в том, что проекционные методы не требуют, чтобы моды были известны точно.

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