Разностные методы для нестационарного уравнения теплопроводности

Рассмотрим задачу для нестационарного уравнения теплопроводности

(8)

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

. (9)

Здесь - некоторые заданные функции. С физической точки зрения данная задача состоит в поиске зависящего от времени распределения температуры в тонком изолированном стержне с начальной температурой, описываемой функцией . Динамика температуры на концах стержня описывается функциями . Внутри стержня имеются источники (стоки) тепла, плотность которых определяется функцией . Коэффициенты характеризуют физические свойства материала стержня: теплоемкость, плотность и теплопроводность соответственно. Если коэффициент теплопроводности не зависит от координат (стержень выполнен из однородного материала), то с помощью соответствующей нормировки можно преобразовать уравнение (8) к виду

(10)

Формально уравнение (8) преобразуется в уравнение (10) путем замены переменных: , .

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

Для приближенного решения поставленной задачи в области построим равномерную сетку : , . Множество узлов сетки представляет собой множество точек на плоскости с координатами ). (См. Рис. 2).

Рис. 2. Равномерная пространственно-временная двумерная прямоугольная сетка и шаблон явной схемы для уравнения теплопроводности.

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

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

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

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

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

.

Аналогично, для первой производной по времени существует возможность использования двух вариантов разностной аппроксимации, а именно использовать разностные производные " вперед" или " назад":

I. ;

II.

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

(11)

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

.

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

(12)

Здесь - значения приближенного решения задачи в узле сетки . Шаблон данной схеме представлен на рис. 2.

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

(13)

Таким образом, нахождение искомого приближенного решения по схеме (12) сводится к непосредственным вычислениям значений во всех внутренних узлах сетки по явным формулам (13). Решение в граничных точках вычисляются непосредственно из краевых условий.

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

, , (14)

.

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

.

Замечание. Системы уравнение (14) является одним из возможных способов аппроксимации дифференциальной задачи (10), (9), при котором дискретизация производится только по пространственной переменной. На таком способе дискретизации основан метод прямых, в котором задача для уравнений с частными производными сводится к решению задачи Коши для системы обыкновенных дифференциальных уравнений с последующим использованием методов интегрирования ОДУ. Название метода прямых происходит из геометрического смысла компонент решения системы ОДУ вида (14), которым соответствуют приближенное решение исходной задачи на прямых в плоскости (x, t).

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

, . (15)

В данном случае погрешность аппроксимации разностной схемы (15) на решении дифференциальной задачи (10), (9) также имеет порядок

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

, , (16)

, .

Решение системы уравнений (16) находится методом прогонки (условие устойчивости метода прогонки легко проверяются непосредственными вычислениями). Разностные схемы, для которых решения не может быть вычислено по явным формулам, а сводится к некоторой системе алгебраических уравнений (выражается системой ЛАУ) принято называть неявными.

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

,

Сравнение полученной разностной схемы с каноническим выражением двухслойной разностной схемы

(17)

позволяет найти вид оператора . Оператор - разностный оператор второй производной. Отсюда легко получить вид оператора для чисто неявной или явной схем. Для явной схемы , для чисто неявной схемы .

Согласно общей теории устойчивости, разностная схема (17) с самосопряженным оператором A=A*> 0 устойчива, когда , откуда следует абсолютная (при любых положительных шагах сетки) устойчивость схемы с весами при (Здесь использован факт, что , E – единичная матрица.

Схема с весами имеет следующий порядок аппроскимации в зависимости от параметра :

1) при , или , если ;

2) при любом , , например, или , если ;

3) при и , если , ;

[1] Программная реализация перестановок строк матрицы с помощью современных возможностей объектно- ориентированного программирования может быть реализована более эффективно, нежели непосредственное умножение с матрицей перестановок.