Неявная схема Эйлера

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

Построение неявной разностной схемы

Чтобы построить неявную разностную схему для уравнения диффузии, используем шаблон, изображенный на рис. 13.11, т. е. для дискретизации пространственной производной будем брать значения сеточной функции с верхнего (неизвестного) слоя по времени. Таким образом, разностное уравнение для (i,k)-ro узла будет отличаться от уравнения для явной схемы (8) только индексами по временной координате в правой части:

Рис. 13.11. Шаблон неявной схемы для уравнения теплопроводности

Если привести подобные слагаемые, то получится система уравнений, связывающая для каждого 1-го узла три неизвестных значения сеточной функции (в самом этом узле и в соседних с ним слева и справа узлах). Множители при неизвестных значениях сеточной функции в узлах шаблона показаны на рис. 13.11 в виде подписей, подобно тому как это было сделано для явной схемы (см. рис. 13.6).

Очень важно, что если само уравнение теплопроводности линейно, то с в левой части разностного уравнения является константой, а ф в его правой части может зависеть только от первой степени и. Поэтому система уравнений (10) для всех пространственных узлов 1=1...м-1 является линейной системой, что существенно упрощает ее решение (поскольку известно, что для линейных систем с ненулевым определителем решение существует и является единственным). Напомним, что для получения замкнутой системы линейных уравнений необходимо дополнить данный набор разностных уравнений граничными условиями, т. е. известными значениями сеточной функции для 1=0 и i=M.

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

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

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

Листинг 13.2. Неявная схема для линейного уравнения теплопроводности

Алгоритм прогонки

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

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

Основным вычислительным ядром программы, реализующей на Mathcad неявную разностную схему, было решение (на каждом временном слое) системы линейных алгебраических уравнений, задаваемых матрицей А. Заметим, что эта матрица, как говорят, имеет диагональное преобладание, а точнее, является трехдиагоналъной (рис. 13.13). Все ее элементы, кроме элементов на главной диагонали и двух соседних диагоналях, равны нулю. С точки зрения оптимизации быстродействия алгоритма, применение встроенной функции isolve является весьма расточительным, поскольку основной объем арифметических операций, выполняемых компьютером (а он составляет, как нетрудно убедиться величину порядка M²), сводится к непроизводительному перемножению нулей.

Рис. 13.12. Решение линейного уравнения теплопроводности при помощи неявной схемы на первом слое по времени (листинг 13.2)

Для отыскания решения линейных систем алгебраических уравнений имеется чрезвычайно эффективный алгоритм, называемый прогонкой, который позволяет снизить число арифметических операций на целый порядок, т. е. до значения порядка м. Это означает, что при использовании пространственных сеток с юоо узлами выигрыш во времени вычислений составит величину порядка 10³! Реализация данного алгоритма приведена в листинге 13.3, который является продолжением листинга 13.2, используя определенные в нем коэффициенты матрицы А, а также начальное условие.

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

Рис. 13.13. Матрица системы линейных разностных уравнений для неявной схемы (листинг 13.2 для М=10)

Листинг 13.3. Алгорктм прогонки (продолжена листинга 13.2)