Повторяя описанную процедуру исключения N - 1 раз, приходят к следующей преобразованной системе уравнений:

(1.15)

коэффициенты которой хранятся в соответствующем массиве машинной памяти.

После завершения преобразований из посаеднего уравнения системы (1.15) находят 11/, затем, используя его, нз предпоследнего уравнения определяют г/л/-1 и т. д., т. е. из к-го уравнения системы (1.15) определяют по формуле

Л"

yk-t

fck

(1.16)

в которой значения и,,-: ыуже найдены ранее, а коэффициенты Qft/" И 6 берутся из массивов, сформированных при проведении исключения.

Внешне процедура выглядиточень просто. Сложности, возникающие при ее реализации, связаны с накоплением погрешностей округления, обусловленных ограниченным числом знаков в ЭВМ. В частности, при вычислении преобразованных коэ([)фициентов по формуле (1.11) при малом значении второе слагаемое может стать столь большим, что при вычитании значения а и Ь вообще «потеряются». Например, если ЭВМ оперирует числами с шестью десятичными знаками и 1, а Qtajla -10 то a\j auj - akidiyldu =

1 ООО 001 будет округлено до Ю**.

Чтобы уменьшить погрешности округления при реализации fe-ro шага исключения, берут соответствую1цее уравнение и неизвестное не в «естественном» порядке, как это бьгло в рассмотренном выше алгоритме, а находят их в результате специального поиска. Цель поиска определить уравнения с максимальным коэффициентом а Такой прием называют выбором ведущего элемента. При этом усложняется алгоритм пересчета коэффициентов уравнений, поскольку приходится как бы переставлять строки и столбцы в матрице линейной системы, чтобы найденный максимальный коэффициент оказался на ее главной диагонали. Эта процедура реализована в стандартных подпрограммах. Поэтому для решения линейной системы по методу Гаусса не следует самому составлять программу, используя простейнше формулы типа (1.11), а целесообразно брать какую-нибудь стандартную программу, в которой разработчики уже предусмотрели меры для уменьшения влияния погрешностей округления.

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

например, симметричные матрицы, у которых ~ flji, и ленточные -матрицы, у которых отличны от нуля только коэффициенты, лежащие внутри «ленты», расположенной вдоль диагонали (рис. 1.2). Для таких случаев можно повысить эффективность метода Гаусса, если -учесть особенности матрицы. Так, например, в ленточных матрицах

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

-учтены. Они будут рассмотрены в § 1.3.

Итдмционные методы. Для изложения итерационных методов

перепишем систему (1.8) в следующем виде, выразив из каждого i-ro

уравнения неизвестное Ы; через оста.пьные

неизвестные:

/ = 1. ITJ

(1.17)

Любой итерационный процесс иачинает--ся с задания по некоторым соображениям начального приближения {un-jnUi- На- Р"*- -

Т!ример, часто полагают w, = 0, п - \, ...

N. В методе простой итерации иа s-m шаге итерационного процесса приближения «** рассчитывают из выражений (1.17), в которых все остальные неизвестные в правой части берут с предыдущей итерации:

Л"

(1.18)

Можно доказать, что достаточным условием сходимости этого метода является условие

Причем хотя бы для одного с-го уравнения должно выполняться строгое неравенство.

Анализ показывает, что итерационный процесс сходится быстрее,

если прн расчете uf в правой части (1.18) для неизвестных с номерами /<: i подставлять не значения wj-, а значения уже найденные на данном шаге. Такой итерационный метод называется методом Гаусса-Зейделя и записывается в виде

2 ctjuy--+du i = \, N. (1.20)

Существует еще одна модификация итерационного метода, позволяющая ускорить сходимость. В методе Гаусса-Зейделя изменение неизвестного при переходе от (s - 1)-й итерации к s-й равно

(1.21)

Скорость изменения значений uf в итерационном процессе можно увеличить или уменьшить, если итерационную схему записать в виде

/1-1

Л"

(1.22)

где множитель a, как показывает дополнительный анализ, следует брать из интервала [О, 2]. При а <; 1 значение иР изменяется медленнее, чем в методе Гаусса-Зейделя, а при а> 1 - быстрее.

На первый взгляд кажется, что для повышения скорости сходимости целесообразно использовать только а> 1. Однако часто пове-

дение и

при увеличении номера итерации s носит колебательный

характер, при котором и"\ переходит от значений, меньших Uj. к значениям, больишм точного решения. При этом чрезмерное увеличение а может лишь увеличить такие колебания. Поэтому в ряде случаев для ускорения сходимости решения имеет смысл «тормозить» изменение значений и*;\ задавая а<; 1. В большинстве реальных задач оптимальное значение множителя а подбирают путем числеи- ных экспериментов. ;

Метод, задаваемыйформулой (1.22), называется методом после- довательной верхней релаксации при а >> 1 или методом последовательной нижней релаксации при а < 1. При а - I получаем как частный случай метод Гаусса-Зейделя.

К достоинствам рассмотренных итерационных методов следует отнести простоту их программной реализации и отсутствие принципиальной необходимости хранения в памяти ЭВМ всех коэффициентов матрицы. Действительно, при вычислении очередного приближения и\ согласно (1.22) нужны только отличные от нуля коэффициенты 1-й строки Qij, bi, которые в принципе могут каждый раз вычисляться заново по исходным данным решаемой задачи. Это обстоятельство обусловливает широкое применение итерационных методов для систем с сильно разреженными матрицами большой размерности, в которых большинство элементов нулевые. Причем это делается как для матриц иеленточной структуры, у которых ненулевые коэффициенты «разбросаны» по всему полю, так и для некоторых ленточных

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

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

Решение систем нелинейных алгебраических уравнений. Ограничимся изложением только двух методов решения, рассматривая их применительно к нелинейным системам частного, но наиболее часто встречающегося в разных теплофизических задачах «квазилинейного» вида. Такие системы записываются аналогично (1,8), но имеют коэффициенты зависящие от искомых величин {и.}-. а/- = aij{Ui, .... Uftf). Они возникают, например, при решении стационарных уравнений теплового баланса (1.2), в которых тепловые проводимости зависят от температур Т, Tj. Для решения этих нелинейных систем обычно применяют итерационные методы, в которых на каждой итерации решается линеаризованная система, т. е. некоторая линейная система, полученная из исходной нелинейной задачи. Наиболее часто применяют два подхода к линеаризации.

Первый из них состоит в том, что на каждом s-m шаге итерационного процесса коэффициенты линеаризованной системы aJ вычисляют по значениям неизвестных, найденным на предыдущей (s - - ])-й итерации: af) = (н<7\ .... "*), а затем путем решения полученной линейной системы с известными a<f} находят новое приближение {и<р}-. Описанный метод называется л/ето(?ож/1/70С/ШЙ итерации для решения квазилинейных систем (не путать с рассмотренным выше методом простой итерации для линейных систем).

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

(a„(Ui. %)u, + o,2(Ui, Щ)ЧгЬх, 1

предположим, что найдены приближенные значения искомых величин u<y> u->na(s-])-й итерации. Построим линейную систему уравнений для приближенных значений ц<*, и*2> следующей й итерации. Для этого представим новые значения в виде = w<«-) 4- Ды(5) u(s) = ы(«-0 + Дц<2), где Ди<»)--уточняю-

щие приращения для приближении неизвестных на s-й итерации. Очевидно, что определение Дц<р эквивалентно определению ир. В методе Ньютона система линейных уравнений обычно записывается относительно приращений Ди*?>. Для ее получения значения коэффициентов а.) представим, используя разложение в ряд Тейлора в точке (ы(у), и ограничиваясь его первым членом, в следую-

щем виде:

(1:24)

Подставив теперь выражения для коэффициентов вида (1.24) и = u(s--i) -j- Дц(р в исходную систему (1.23) и отбросив члены второго порядка малости (произведения Лы<) • Аи)), получим нужную нам систему линейных уравнений относительно приращений

s-1)

(Ь25>

"з 1

"22

й«1

"21 "1 "22 "г

Au(p =

Линейность системы (1.25) относительно приращений Au(f\ Ди(> следует из того, что все ее коэффициенты рассчитаны по значениям искомых величин на предыдущей итерации. На каждом шаге итераций можно использовать для решения (1-25) какой-либо прямой метод. Легко увидеть, что выражения в правых частнх уравнений системы (1.25) представляют собой невязки для исход1ЮЙ системы (1.23) при значениях ы<-) и Если итерационный процесс со-

шелся, то в правых частях получим нули и решение системы (1.25) даст Ды<р Ды<) = 0.

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

1 3. ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ

В математическом обеспечении ЕС ЭВМ имеется пакет прикладных программ, предназначенных для решения систем линейных алгебраических уравнений 115]. Подпрограммы написаны на ФОРТРАНе и могут быть использованы не только на ЕС ЭВМ, но и иа других типах ЭВМ. Эти подпрограммы реализуют прямые методы какдля матриц общего вида, так и для матриц специального вида (симметричных, ленточных). Ниже рассмотрим несколько широко применяемых подпрограмм, которые далее будут использованы при решении задач теплопроводности, лучистого и конвективного теплообмена.

Чтобы решить систему линейных уравнений с помощью какой-либо стандартной подпрограммы, пользователь должен составить головную (вызывающую) программу, в которой элементы матрицы А и стапбца правых частей В линейной системы АХ = В записываются в некоторые массивы, а затем выполняется вызов стандартной подпрограммы. При работе со стандартными подпрограммами из пакета И51 начинающие программисты часто допускают некоторые типичные ошибки, связанные с формированием массивов, в которые .записываются элементы матриц. Например, такие ошибки возникают, когда в вызывающей программе матрица формируется в виде двумерного массива А, предельные размеры которого, установленные в операторе DIMENSION, превышают фактические размеры М X М. Остановимся на данном вопросе подробнее.

Формы представления матриц. Входными параметрами рассматриваемых ниже подпрограмм являются массивы А и В, содержащие элементы матрицы А и столбца В, расположенные в строго опреде-лен}1ой последовательности, число уравнений М, а также некоторые дополнительные параметры для матриц специального вида. Эти стандартные подпрограммы позволяют решать системы с произвольным числом уравнений Л1, поскольку число М и массивы А, В входят в число формальных параметров подпрограммы, а фaктгчecкнe размеры массивов устанавливаются в головной программе. Таким образом стандартная подпрограмма оперирует с матрицей А как с массивом переменной длины /И х М и «не знает» о предельных размерах массивов, опред&тенных в головной программе пользователем в операторе DIMENSION. Прн этом элементы матрицы А должны быт расположены в массиве А подряд в определенной последовательности. Например, для матрицы общего вида в соответствующей области машинной памяти последовательно по столбцам должны быть записаны М элементов: а,;. a2i, ам-i, амм-

В головной программе матрицу часто бывает удобно формировать в виде двумерного массива А, первый индекс которого соответствует номеру строки, а второй - номеру столбца. Рассмотрим форму хранения двумерных массивов в памяти ЭВМ. Пусть, например, для 1*"ения матрицы А предусмотрен двумерный массив, описанный