Решение системы уравнений с трехдиагональной матрицей

Решение систем с трехдиагональной матрицей

Трехдиагональная линейная система уравнений может быть записана в следующем виде:

Решение ищется в виде:

u_i-1=bet_i*u_i+gam_i В этом случае для bet и gam получаются рекуррентные формулы: bet_i+1=c_i*fi, gam_i+1=(f_i+a_i*gam_i)*fi, fi=1/(b_i-a_i*bet_i); bet₀=gam₀=0. По этим формулам вспомогательные массивы вычисляются для i от 1 до N, затем находится u_N-1=gam_N, а затем при известных вспомогательных массивах вычисляются все значения u.

Алгоритм принципиально не может иметь ведущего элемента, и есть вероятность того, что он не сойдется даже для несингулярной матрицы. Для этого в программе, реализующей прогонку, необходимо контролировать значения fi. Еще два замечания: длины рекуррентных цепочек порядка N, поэтому при экспоненциальном накоплении погрешностей (что происходит при преобладании значений |bet_i|>1) прогонка практически обязательно разойдется. Второе: алгоритм сохраняет исходные значения коэффициентов. Доказано, что для устойчивой работы алгоритма Томаса достаточно диагонального преобладания, однако, во многих случаях он сходится и при отсутствии такового.

Общее описание алгоритма таково.
Если провести исключение из имеющейся системы уравнений всех уравнений, находящихся на нечетных позициях, то от системы N исходных уравнений мы перейдем к системе из (N-1)/2+1 уравнений только для четных значений u, нечетные же значения будут вычисляться через четные. Повторяя эту процедуру для укороченной системы, мы в конце концов придем к единственному уравнению для узла, лежащего в позиции (N-1)/2 (только если N-1 — степень двойки). Затем, проводя обратную подстановку, определяются крайние, а затем и все прочие узлы.

Формулы для пересчета коэффициентов на каждом этапе исключения таковы:

Если число уравнений не является степенью двойки, алгоритм реализуется методом фиктивного восполнения области определения до степени 2. При этом никаких реальных действий, кроме нескольких целочисленных проверок, не производится, а только коэффициенты, индекс которых выпадает из области определения [0,N-1] включительно, считаются равными нулю и в преобразованиях не участвуют (подробности в программе).

Особая устойчивость этого алгоритма по сравнению с алгоритмом Томаса связана с тем, что длина рекуррентных цепочек здесь не превышает log2(N)+1, и в связи с этим даже при экспоненциальном росте погрешности исходный результат будет ограниченным. Кроме того, в алгоритме Томаса в рекуррентных цепочках производятся операции между коэффициентами соседних узлов, при редукции на этапах больших начального — между далеко отстоящими, что уменьшает вероятность образования погрешности при вычитании двух близких чисел. Проверено, что алгоритм редукции успешно работает во многих случаях, когда алгоритм Томаса (прогонка) расходится.

Алгоритм не требует дополнительной памяти, но сохраняет только половину из имеющихся в начале его работы коэффициентов: те, что находятся на нечетных позициях (независимо от четности N). Скорость работы почти не зависит от того, является ли N-1 целочисленной степенью 2.

Ниже приводятся программы для прогонки и редукции.

Линейные уравнения. Решение систем линейных уравнений. Метод прогонки.

Метод прогонки (алгоритм Томаса) используют для решения СЛУ типа Ax=F, где A — трёхдиагональная матрица. Это вариант метода последовательного исключения неизвестных.

Система уравнений Ax=F равноценна соотношению:

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

Выразим x_i-1 и x_i через x_i+1, подставим в уравнение, используя это соотношение, (1):

где F_i — правая часть i-го уравнения.

Это соотношение выполняется не завися от решения, если потребовать:

Получаем из 1-го уравнения:

После того, как нашли прогоночные коэффициенты α и β, используем уравнение (2) и получим решение системы. Причем,

Еще одним вариантом объяснения смысла метода прогонки является такой вариант: преобразуем уравнение (1) к равнозначному ему уравнению:

c надиагональной (наддиагональной) матрицей:

Рассчеты проводим в 2 этапа. На 1-ом этапе вычисляем компоненты матрицы C′_i и вектора F′, начиная с i=2 до i=n:

На 2-ом этапе, для i=n,n−1,…,1 вычисляем решение:

Для применимости формул метода прогонки достаточно свойства строгого диагонального преобладания у матрицы A.

О построении решений трехдиагональных систем линейных алгебраических уравнений специального вида

О ПОСТРОЕНИИ РЕШЕНИЙ ТРЕХДИАГОНАЛЬНЫХ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ СПЕЦИАЛЬНОГО ВИДА

Национальный исследовательский Иркутский государственный технический университет,

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

Ил. 3. Табл. 1. Библиогр. 6 назв.

Ключевые слова: трехдиагональная матрица; метод прогонки; приложения; неявная разностная схема; вычислительный эксперимент.

ABOUT CREATION OF DECISIONS OF SPECIAL THREE-DIAGONAL SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS

T. Batagaeva, Ta Chung Thanh

National Research Irkutsk State Technical University,

83, Lermontov St., Irkutsk, 664074

The article inquires into explicit formulas for decision of special systems of linear algebraic equations applied in engineering and physical tasks. The authors produce results of computing experiment in which the traditional method of such systems decision is compared with the method based on received formulas. The new method is proved to be more effective.

3 figures. 6 sources.

Key words: three-diagonal matrix; Thomas algorithm; applications; implicit incremental diagram; computing experiment

Теория разностных методов является обширной областью вычислительной математики, а системы разностных уравнений представляют собой прекрасный пример приложения трехдиагональных матриц [5]. Задачи, включающие в себя трехдиагональные системы [4] линейных алгебраических уравнений (СЛАУ) в качестве основных или вспомогательных компонент, можно разделить на несколько типов [5]:

· одномерные краевые задачи для дифференциальных уравнений второго порядка (задача Дирихле для уравнения Пуассона, задача Неймана, смешанная краевая задача);

· параболические задачи с одной пространственной переменной (одномерная задача нестационарной теплопроводности или диффузии);

· решение многомерных эллиптических и параболических уравнений (задача Дирихле для уравнения Пуассона в двумерной граничной области);

· решение локально-нелинейных одномерных задач (двухфазная задача Стефана);

Еще одним примером краевой задачи с применением трехдиагональной матрицы может служить задача Штурма-Лиувилля [3], заключающаяся в отыскании нетривиальных решений однородного уравнения на некотором отрезке. В свою очередь, эта задача служит постоянным источником новых идей для спектральной теории операторов и смежных вопросов анализа. Применение трехдиагональных систем не ограничивается задачами математической физики, они встречаются в алгоритмах сплайновых аппроксимаций [4], а так же в задачах линейной алгебры, для решения проблемы собственных значений матриц общего вида.

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

В данной работе предложен новый подход к решению СЛАУ специального вида, матрицы которых имеют одинаковые элементы на главной и побочных диагоналях [6]. Такие матрицы встречаются при решении краевых задач теплопроводности, задачи Дирихле и в некоторых других.

Метод прогонки. Для решения трехдиагональных СЛАУ часто применяют метод прогонки. Основным его преимуществом является экономичность, т. е. то, что этот метод максимально использует структуру исходной системы. Недостатком метода является то, что с каждой итерацией накапливается ошибка округления [1].

Запишем трехдиагональную систему в матричной форме: , где – матрица, у которой все элементы, кроме элементов главной и примыкающих диагоналей, равны нулю, а n – это размерность матрицы, – вектор неизвестных, – вектор свободных членов. Матрица представлена формулой

. (1)

Запишем -ое уравнение системы , где , .

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

, , где и показаны формулами

, , , . (2)

Заметим, что , так как .

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

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

, причем , (3)

При этом хотя бы для одного уравнения должно выполняться строгое неравенство [1].

Доказательство достаточности условия диагонального преобладания. Предположим, что условия (3) выполнены. Тогда

; .

Это показывает устойчивость метода прогонки, поскольку исключена возможность лавинообразного накопления ошибок округления на обратном ходе решения [1]:

Для того, чтобы в самом начале обратного хода знаменатель дроби, определяющей , не обратился в ноль, необходимо, чтобы хотя бы одно неравенство (3) было строгим. Если , т. е. либо , либо , то этого не произойдет.

Достаточность условия диагонального преобладания для устойчивости доказана.

Метод явных формул [6]. Рассмотрим трехдиагональную СЛАУ следующего вида:

(4)

Пусть – определитель матрицы СЛАУ (4). Условие , где вычисляются по формулам

, , , (5)

является, очевидно, необходимым и достаточным для однозначной разрешимости системы (4) [5]. Решение (4) определяется как

(6)

где n – порядок СЛАУ (1); .

Докажем соотношения (6):

СЛАУ (4) можно записать, в векторно-матричном виде , где

Получение явного решения основано на методе Крамера [4]. Для этого найдем определители матриц, в которых один из столбцов заменен на вектор свободных членов . Определитель основной матрицы находится из (5).

Пусть – определитель матрицы , в которой i-ый столбец заменен вектором . Его вычисление проведем разложением по i-му столбцу. Получим, что

Наконец, вычисляется по формуле

Согласно правилу Крамера, решение всякой СЛАУ с ненулевым определителем единственно и определяется как Отсюда (с учетом выражений для полученных выше) следует справедливость (6).

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

Первый или последний элемент системы лучше вычислить по формуле (6), а для вычисления остальных провести ход, аналогичный обратному ходу метода прогонки.

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

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

Уравнение теплопроводности является классическим примером уравнения параболического типа. В одномерном по пространству случае, однородное (без источников энергии) уравнение теплопроводности имеет вид:

, , . (7)

Если на границах и заданы значения искомой функции в виде:

, ; (8)

, , (9)

то (7)–(9) называют первой начально-краевой задачей для уравнения теплопроводности.

В терминах теории теплообмена – распределение температуры в пространственно-временной области , – коэффициент температуропроводности, а (8), (9) с помощью функций , задают температуру на границах и [2]. В качестве примера мы рассмотрим линейные краевые условия:

О вычислительном эксперименте. На примере решения задачи (7)–(9) был проведен вычислительный эксперимент, в котором сравнивалась скорость решения трехдиагональной СЛАУ специального вида по методу прогонки и с помощью полученных явных формул.

Так как скорость вычислений зависит не только от алгоритма, но и от характеристик аппаратного и программного обеспечения машины, на котором производится эксперимент, приведем некоторые данные:

· процессор: Intel Core 2 Duo 2.4GHz;

· ОЗУ: 4.00 GB, DDR3, 1067MHz;

· операционная система: Mac OS x Lion 10.7.3;

Программа написана на языке Java с использованием интегрированной среды разработки NetBeans IDE 7.1.1.

Алгоритм программы. Сначала на вход программы подаются следующие данные:

· коэффициенты левой и правой граничных функций, которые представляют собой полиномы третьей степени:

;

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

· метод решения трехдиагональной СЛАУ – метод прогонки или решение с помощью явных формул;

· границы временного периода;

· количество точек сетки.

Не теряя общности рассмотрения, можно считать, что коэффициент температуропроводности равен единице.

Затем начинается построение трехдиагональной СЛАУ на основе неявных разностных схем (явные схемы не применяются, так как они устойчивы условно). Суть метода заключается в покрытии расчетной области сеткой из точек (рис. 1), что определит шаги по времени и пространству и соответственно. Тем самым определяются узлы, в которых будет осуществляться поиск решения. Затем надо заменить уравнение теплопроводности аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого -го узла сетки, где – порядковый номер точки деления по ; – порядковый номер точки деления по , а – значение функции , соответствующее точке .

Нумерация точек разностной сетки организована следующим образом: по оси ; по оси .

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

Рис. 1. Неявная разностная схема для симметричной задачи

Рис. 2. Шаблон неявной

Неявной разностной схемой называется разностная схема, которая использует четырехточечный шаблон, показанный на рис. 2.

Введем два временных слоя: нижний , на котором распределение искомой функции u(xi, tk) известно (при k = 0 распределение определяется начальным условием (8)–(9), и верхний временной слой , на котором распределение искомой функции подлежит определению.

Затем надо заменить уравнения теплопроводности аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого (i, k)-го узла сетки.

Заменим аппроксимирующие уравнения теплопроводности уравнениями в конечных разностях

, (10)

, (11)

если подставить полученные дифференциальные операторы (10), (11) в исходное уравнение, получим соотношение, аппроксимирующее это дифференциальное уравнение в точке на разностной сетке:

Перенеся все неизвестные в левую часть, получим уравнение неявной разностной схемы:

, где . (12)

Применим схему (12) к симметричной задаче (7)–(9).

Сначала рассмотрим нулевой временной слой, где . Тогда уравнение неявной разностной схемы (12) примет вид:

где известны. Найдем

Перейдем к следующему временному слою, . Получится система уравнения, состоящая из трех уравнений:

Матрица этой системы имеет трехдиагональный вид:

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

Матрица для -ого временного слоя выглядит следующим образом:

, (13)

а вектор свободных членов соответствующей системы линейных алгебраических уравнений примет вид:

. (14)

На основе матрицы (13) и вектора свободных членов (14) представим полученную трехдиагональную СЛАУ:

. (15)

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

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

Анализ результатов. В программе измеряется только время составления матрицы и решения системы. Время чтения входных параметров и вывода результатов не учитывается.

Значения входных параметров:

· коэффициенты левой граничной функции: , , ;