§ 7. Решение параболических уравнений с несколькими пространственными переменными

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

Перейдем к построению и исследованию разностных схем. Пусть требуется найти функцию и, являющуюся решением уравнения

в области с начальными и граничными условиями

здесь .

Попытаемся применить к решению задачи (1), (2) методы, разработанные ранее. Введем в рассмотрение квадратную сетку с шагом :

а на отрезке [0, Т] — сетку с шагом . Будем искать приближенное решение задачи (1), (2) в дискретной системе точек (узлов)

Множество будем называть сеточной областью, а множество точек при фиксированном будем называть слоем.

По аналогии с одномерным случаем построим явную и неявную разностные схемы для задачи (1), (2) и попытаемся выяснить, в чем заключается принципиальное отличие от случая одной пространственной переменной. Заменяя в узле разделенной разностью или , а - выражением (см. § 6)

получаем разностные схемы: явную

и неявную

При использовании схемы (3) счет ведется по явным формулам — по известным значениям из (3) находятся значения . Поэтому проблем с реализацией алгоритма на ЭВМ не возникает. Остается лишь исследовать устойчивость этой схемы.

По-другому обстоит дело в случае схемы (4). Относительно , имеем систему линейных алгебраических уравнений, т.е. схема (4) неявна. Структура матрицы уравнений (4) совпадает со структурой матрицы оператора (см. § 6). Поэтому при решении этой системы возникают те же трудности, что и в случае эллиптических уравнений. Напомним, что в одномерном случае проблема численного решения уравнений на верхнем слое не возникала, так как можно было воспользоваться методом прогонки.

Введем понятие экономичной разностной схемы. Разностную схему, аппроксимирующую задачу со временем, называют экономичной, если она безусловно устойчива и при переходе от слоя к слою требуется количество арифметических операций, пропорциональное числу узлов на слое. (Иногда условие безусловной устойчивости в определении экономичной схемы отсутствует.) Из определения следует, что чисто неявная схема для одномерного уравнения теплопроводности является экономичной.

Перед тем как заниматься построением экономичных разностных схем, исследуем устойчивость разностной схемы в общей постановке.

Введем пространство Н функций, определенных на ; - значение функции в узле . Скалярное произведение и норму в Н определим как

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

где В и А — симметричные положительно определенные операторы, отображающие Н в себя. Как и в одномерном случае, иногда будем различать устойчивость по начальным данным и по правой части.

Обозначим . Учитывая равенство , преобразуем (5) к виду

Положим . Так как и оператор симметричен, то

Поэтому, умножая обе части (6) скалярно в на , получим

Напомним, что по определению .

Будем считать, что

Тогда (7) можно переписать в виде

Соотношение (9) является энергетическим тождеством.

Так как пространство конечномерно, то существует постоянная к такая, что для всех или, что то же самое,

Так как , то существует оператор (матрица) симметричный и положительно определенный такой, что . Через обозначим оператор . В рассматриваемом случае

и из (9) следует неравенство

Фиксируя , например, полагая , отсюда получаем соотношение

связывающее нормы функции на соседних слоях. Таким образом,

Имеем последовательность неравенств

Предполагая, что не зависит от шагов сетки , при получаем окончательное соотношение

которое означает устойчивость разностной схемы (5) по начальным данным и по правой части. Таким образом, условие является достаточным для устойчивости разностной схемы по правой части и начальным условиям.

Заметим, что выражение , стоящее в правой части (12), является квадратурной формулой для интеграла , а выражение аналогом нормы .

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

Применим к обеим частям (5) оператор ; получим

Обозначим . Тогда и последнее равенство примет вид

Оператор симметричен, поэтому собственные числа S лежат на отрезке [, ], где

Таким образом, собственные числа S не превосходят , если выполняется соотношение

Преобразуя выражение и обозначая , получаем

Поэтому (13) имеет вид

Правая часть неравенства выполнена всегда, поскольку левая часть эквивалентна выполнению для любого , неравенства

Последнее означает, что

Как видим, условие (14), обеспечивающее устойчивость схемы (5), не совпадает с (8). Однако у схем, обладающих свойством (8), имеется существенное достоинство. В ряде практических задач интервал интегрирования Т уравнения (1) достаточно велик или же требуется вести счет задачи до выхода на стационарный режим. В этих случаях целесообразно использовать разностные схемы, удовлетворяющие более сильной оценке устойчивости (15) при . Обозначим и . Тогда из (11) имеем

Пусть . Тогда из последнего неравенства получаем

Отсюда следует, что при решение задачи (5) будет ограничено на бесконечном промежутке времени при выполнении условия (8). Отметим, что из (14) ограниченность решения на бесконечном промежутке времени при наличии правой части, вообще говоря, не следует.

Разностную схему (5) можно рассматривать как итерационный процесс решения уравнения

где . В этом случае выполнение условия обеспечивает сходимость итерационного процесса. Действительно, записывая уравнение для погрешности , имеем

Из условия следует, что определяет в Н норму, которую будем обозначать . Тогда из (9) получаем

Отсюда

Введем норму ; тогда из (16) следует окончательная оценка

Таким образом, итерационный метод (5) сходится со скоростью геометрической прогрессии. При этом скорость сходимости определяется величинами тики условием . Процесс сходится как в норме, определяемой оператором D, так и в норме, определяемой оператором А.

Выясним, при каких условиях будут устойчивы схемы (3), (4). Схема (3) уже имеет вид (5), при этом . Необходимо, таким образом, проверить выполнение условия . Имеем (см. § 6)

Таким образом, чтобы (14) было справедливо, достаточно выполнения неравенства , т.е. явная схема (3) условно устойчива при .

Представим теперь неявную схему (4) в виде (5). В этом случае

Таким образом, условие (14), равно как и условие (8), выполнено при любых и h, т. е. неявная схема (4) безусловно устойчива.

Если для решения системы уравнений относительно значений решения на верхнем слое применяется так называемый марш-алгоритм в его устойчивой форме, то число арифметических операций при переходе от слоя к слою пропорционально числу неизвестных. Тогда неявная схема является экономичной.

Перейдем к изучению других экономичных разностных схем для уравнения (1). Будем рассматривать задачу (1) с однородными граничными условиями, т. е. при . Пусть и - операторы второй разделенной разности по направлениям соответственно, т. е.

Здесь, как и ранее — функция, совпадающая с на , и равная нулю на .

Задача 1. Проверить, что функция , где —произвольная функция аргумента , является собственной для оператора , а любая функция , где - произвольная функция аргумента собственной для оператора .

Задача 2. Проверить, что функции образуют полную систему собственных функций операторов .

Положим

Оператор В является симметричным и положительно определенным как произведение симметричных положительно определенных и коммутирующих между собой операторов. Операторы такого вида называют . Коммутируемость операторов и можно проверить непосредственно; кроме того, она следует из того факта, что эти операторы имеют общую полную систему собственных функций (см. задачу 2) и, следовательно, записываются в виде

где — диагональные матрицы, матрица С одна и та же. Проверим, при каких , выполняется условие (14). Имеем

поэтому условие (14) приобретает вид

Так как оператор положительно определен, то условие обеспечивает выполнение (14), т.е. при разностная схема (5) безусловно устойчива по начальным данным.

Рассмотрим алгоритм реализации схемы (5) в данном случае. Обозначим и представим (5) в виде

Последнее уравнение разобьем на два:

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

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

Примечание. Если решается неоднородная краевая задача, т. е. то, вообще говоря, не удовлетворяют граничному условию и значения требуется определять специальным образом.

Аналогично второе уравнение (17) расписывается в виде

При фиксированном (19) является системой уравнений с трехдиаго-нальной матрицей относительно неизвестных .

Рис. 10.7.1

Рис. 10.7.2

Таким образом, функция z может быть найдена из (19) за арифметических операций. Группы неизвестных, связанные уравнениями (19), объединены на рис. 10.7.2 символом . Значения находятся по явным формулам:

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

Описанный алгоритм часто называют методом расщепления.

Как мы видели, основная его идея заключается в построении оператора при разностной производной по времени таким образом, чтобы этот оператор являлся произведением операторов , каждый из которых действует только в одном направлении, и полученная схема аппроксимировала исходную задачу. Так, в данном случае .

Задача 3. Доказать, что при достаточно гладком решении для рассматриваемой схемы имеет место аппроксимация порядка . Таким образом, при погрешность аппроксимации имеет порядок и схема абсолютно устойчива.

Близким по своей структуре является метод переменных направлений. Суть его заключается в переходе от к по формулам

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

Построение методов с расщепляющимся оператором на верхнем слое в k-мерном случае можно провести по точно такой же схеме. Пусть , где — одномерный оператор в направлении. Положим . Нужная разностная схема будет иметь вид

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

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

Другая разностная схема может быть получена из следующих соображений. Представим А в виде , где — правая и левая треугольные матрицы. В частности, если , то можно записать следующим образом:

Рассмотрим разностную схему

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

Пусть в качестве примера рассматривается уравнение (1) в прямоугольнике с нулевыми граничными условиями. В этом случае определяются по формулам (23). Тогда . Оператор является положительно определенным, поэтому . Условие в данном случае заведомо будет выполнено при . Таким образом, условие обеспечивает безусловную устойчивость схемы (24).

Найдем порядок аппроксимации разностной схемы (24). Так как

то, используя формулу Тейлора, получим

Используя аналогичную оценку погрешности аппроксимации оператора , приходим к заключению, что выражение — от аппроксимирует единичный оператор Е с порядком . Если считать, что порядка 1, то схема (24) аппроксимирует исходное уравнение (1) с порядком . Таким образом, величина должна быть достаточно мала. Тем не менее схема (24) все же лучше, чем явная. Для устойчивости явной схемы требуется выполнение условия , в то время как при схема (24) будет безусловно устойчивой и шаг по времени может быть выбран существенно большим. Например, можно взять или же . Погрешность аппроксимации по времени в этих случаях будет иметь порядок соответственно. Следует, однако, отметить, что часто решение параболического уравнения по t обладает достаточным количеством производных и эти производные стремятся к нулю при .

Это оправдывает применение схемы (24) при расчете нестационарных задач, поскольку шаг по времени можно брать достаточно большим.

Выпишем алгоритм, соответствующий разностной схеме (24). Обозначим через . Тогда, решая уравнение

можно найти значение по явным формулам. Действительно, в случае, когда - квадрат (что несущественно), уравнение (25) в узле имеет вид

(26)

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

Для вычисления значения требуется число арифметических операций, не зависящее от шагов сетки. Поэтому вычисление значений функции во всех узлах потребует арифметических операций, что по порядку совпадает с количеством узлов на слое. Аналогичным образом решая уравнение

за арифметических операций найдем значение . Решение находится после этого простым пересчетом по формуле

Таким образом, переход от в схеме (24) требует числа арифметических операций, пропорционального количеству узлов сетки на слое, т. е. схема (24) экономична.

Наиболее эффективна схема (24), если ее рассматривать как итерационный метод для решения стационарной задачи . В этом случае требование аппроксимации по t не играет никакой роли и параметр можно выбирать только из соображений наиболее быстрой сходимости итерационного метода. Обычно выбирают , чтобы выполнялось (8). Тогда операторы В и А связаны соотношением

где не зависят от h, и, решая стационарную задачу с помощью итерационного процесса (24), получим решение с точностью за арифметических операций.

Задача 4. Доказать оценку (27) при .

Задача 5. Показать, что при любых отношение ограничено снизу постоянной, отличной от нуля.

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

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

а правую часть — в виде суммы правых частей: . Левая и правая части уравнения (1) равны сумме левых и правых частей уравнения

Опишем переход от к слою. Аппроксимируем первое из уравнений (28) следующим образом:

Второе уравнение (28) заменим соотношением

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

Ясно, что каждое из уравнений (29), (30) не аппроксимирует исходную задачу. Найдем погрешность аппроксимации. Имеем

Аналогично

В общем случае , поэтому и уравнения (29), (30) аппроксимируют (1) с порядком . Однако

В этом случае говорят, что схема (29), (30) аппроксимирует задачу (1) в суммарном (или слабом) смысле, т.е. хотя каждое из уравнений (29), (30) не аппроксимирует исходную задачу, сумма погрешностей аппроксимаций этих уравнений равна .

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

Теорема (без доказательства). Схема (29), (30) устойчива в сеточной норме пространства С и при достаточно гладком решении

где не зависит от h и , а значение решения на слое.

Рассмотренный метод получения разностных схем носит название метода дробных шагов или же метода суммарной аппроксимации. Его можно применять не только в линейных задачах, но и в нелинейных. В общем случае для уравнения

где операторы , вообще говоря, нелинейные и не обязательно одномерные, схема метода дробных шагов заключается в следующем. Решение на шаге уравнения (31) заменяется последовательным решением на шаге уравнений

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

В создание описанных выше экономичных методов решений многомерных нестационарных задач внесли большой вклад Е. Г. Дьяконов, Г. И. Марчук, А. А. Самарский и Н. Н. Яненко.