8.3. Конечно-разностные методы решения уравнений в частных производных

8.3.1. Построение конечно-разностных схем

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

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

;

, .

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

;

, ;

или

;

, .

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

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

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

Рассмотрим наиболее известные варианты конечно-разностных схем, используемые при решении задач математической физики (например, [Андерсон,Дьяченко,Емельянов, ГодуРябенький, Роуч, Флетчер.]). Для упрощения будем рассматривать уравнения в частных производных, записанных в нестационарной одномерной постановке.

Уравнение переноса

(8.21)

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

1. Явные схемы первого порядка точности по времени и по пространственной координате (шаблоны, представленные на рисунках 8.7а, 8.7б):

; (8.22)

. (8.23)

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

2. Схема «тренога» - явная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7в):

; (8.24)

3. Схема Лакса - явная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7г):

; (8.25)

4. Схема «чехарда» - явная схема второго порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7д):

; (8.26)

5. Схема Лакса-Вендрофа - явная схема типа «предиктор-корректор» второго порядка точности по времени и второго порядка по пространственной координате:

первый шаг (предиктор) – осуществляется по формуле Лакса (8.25);

второй шаг (корректор) – выполняется по формуле

; (8.27)

6. Схема Мак-Кормака - явная схема типа «предиктор-корректор» второго порядка точности по времени и второго порядка по пространственной координате:

на первом шаге (предиктор) вычисляется промежуточное значение сеточной функции по формуле вида (8.22)

;

на втором шаге (корректор) вычисляется окончательное значение сеточной функции по формуле вида (8.23)

. (8.28)

7. Неявная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.7е):

(8.29)

Уравнение диффузии (теплопроводности)

. (8.30)

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

1. Явная схема первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.8а):

; (8.31)

2. Явная схема Дюфорта-Франкела первого порядка точности по времени и второго порядка по пространственной координате (шаблон представлен на рисунке 8.8б):

; (8.32)

3. Неявная конечно-разностная схема запишется в виде (шаблон представлен на рисунке 8.8в):

; (8.33)

4. Неявная конечно-разностная схема Кранка-Николсона (шаблон представлен на рисунке 8.8г):

. (8.34)

В формуле (8.34) - . В частности, можно принять .

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

Рассмотрим решение уравнения переноса, выполняемого по явной схеме (8.22). Уравнение приводится к виду

. (8.35)

Определение сеточной функции по уравнению (8.35) устанавливается следующим алгоритмом:

• примем значение n=0 и установим шаг интегрирования , обеспечивающий устойчивость последующих вычислений (порядок выбора шага интегрирования уравнений по времени, обеспечивающих устойчивый расчет конечно-разностных уравнений, рассматривается ниже в разделе 8.3.3);

• в соответствии с начальными условиями определяем значения сеточной функции для всех значений ;

• в соответствии с граничными условиями устанавливаем значение сеточной функции ;

• начиная с до в соответствии с уравнением (8.35) вычисляем значения ;

• увеличиваем значение n на единицу и возвращаемся к началу алгоритма.

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

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

. (8.36)

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

;

;

. (8.37)

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

Следует дополнительно отметить следующее.

1. Выше рассматривались варианты конечно-разностной аппроксимации уравнений, в которых присутствуют частные производные по времени и по пространственным координатам от одной и той же функции. Однако на практике возникает необходимость в решении уравнений вида (8.7), в которых каждая из частных производных, входящая в уравнение, определяется от своей функции. Эти функции, как правило, связаны между собой. Для подобных уравнений применимы все рассмотренные выше приемы конечно-разностной аппроксимации. Например, при решении уравнения

(8.38)

можно применить схему Лакса в следующей конечно-разностной форме

. (8.39)

В формулах (8.38), (8.39) g – может быть функцией переменных t, x ( ) или t, x, f ( ).

2. Аналогичные рассмотренным могут быть построены явные конечно-разностные схемы решения многомерных задач, описываемых уравнениями в частных производных. Сеточную функцию в этом случае обозначают, например, в виде . Здесь i, j, k – целочисленные координаты сеточной функции по трем координатным направлениям.

Рассмотрим трехмерный вариант уравнения переноса

и запишем для его решения конечно-разностную явную схему Лакса

.

Решение последнего уравнения остается столь же простым, что и решение уравнения (8.25).

3. Применение неявных схем для многомерных задач существенно усложняет их решение. Рассмотрим, например, многомерный вариант уравнения диффузии

. (8.40)

Для этого уравнения можно записать следующую неявную конечно-разностную схему

(8.41)

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

.

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

Существуют эффективные вычислительные методы решения многомерных задач, называемые методами пространственного расщепления, методами дробных шагов, суть которых состоит в сведении решения многомерной задачи к совокупности одномерных задач. Для уравнения (8.40) применение этого метода требует решения трех систем уравнений – (8.42), (8.43), (8.44).

, (8.42)

, (8.43)

. (8.44)

Каждая из систем уравнений вместе с начальными и граничными условиями может быть приведена к виду (8.33) с трехдиагональной матрицей А.

4. Исходное дифференциальное уравнение в частных производных можно привести к виду системы обыкновенных дифференциальных уравнений. Например, для случая уравнения диффузии (8.30) при явной аппроксимации второй производной по пространственной координате такая система может быть записана в виде

. (8.45)

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

5. В последние годы особый интерес приобретают конечно-разностные схемы высоких порядков точности по времени и по пространственным координатам. Построение таких схем рассматривается, например, в работах [Марчук, МарчукШайдуров, Липанов].