В простейшем случае уравнения формируются полиостью. При таком подходе требуется Л ячеек памяти, где Л -количество уравнений, и его можно использовать только для небольших задач и при применении больших ЭВМ. При 100 уравнениях, например, требуется 10 000 ячеек памяти.
Для решения систем линейных уравнений хорошо разработаны два метода: а) прямой метод, позволяющий получить точное (в пределах ошибки округления) решение; б) итерационный метод, использующий сходящийся к точному решению процесс последовательных приближений.
В качестве прямого метода часто используется метод исключения Гаусса, а в качестве итерационного - метод Гаусса - Зей-деля.
20.5.1. Метод исключения Гаусса
Пусть система Л/ уравнений представлена а блочном виде
Ки К\2
- 21 Кп -
-б2.
(20.1)
где АГц -матрица размерности 1 X 1,
K\i - матрица размерности 1 X (Л - 1). К21 - матрица размерности (Л/ - 1) X U К22 - матрица размерности (Л/ - 1) X (ЛГ - 1), 6 -вектор неизвестных величин, F - вектор известных величин. Метод исключения Гаусса позволяет уменьшить размерность матрицы К и получить матричное уравнение размерности - 1 в виде
Гд = Г, (20.2)
К* = - К2\К\1 К\2, F = р2 - KiiKwFi
(20.3) (20.4)
Разбивая точно так же на блоки матрицу К*, этот прием можно првторить. Основной операцией прн этом будет вычисление тройного произведения KiiKTiKn. Так как матрица Км имеет размерность 1 X 1, количество операций пропорционально (Л- 1). Когда размерность матрицы К будет сведена к 1 X 1, последнюю неизвестную 8„ можно будет найти непосредственно.
Используя обратный ход), остальные неизвестные можно
) Процесс нсключення Гаусса обычно называют прямым ходом, а решение системы с треугольной матрицей- обратным ходом. -Ярнл. ред.
определить из уравнения типа
Si = KiiFi - КчКА.
(20.5)
Описанная процедура представляет собой простейший прямой метод, при использовании которого все элементы матрицы хранятся в оперативной памяти и участвуют в вычислениях. Используя свойство симметрии матриц жесткости, можно сократить число операций почти вдвое, производя все действия только с частью матрицы, расположенной выше диагонали. При этом для решения необходимо выполнить около /еЛ операций.
Типичная матрица жесткости ансамбля элементов содержит много нулей, в частности, на некотором расстоянии от диагонали находятся только нулевые члены. Такая матрица называется ленточной, а расстояние от диагонального члена элемента до последнего ненулевого элемента этой же строки называется половиной ширины ленты. Ленточный характер этой матрицы можно продемонстрировать, записывая ее в блочном виде, как показано ниже, где нулевая подматрица заменяет часть подматрицы Кц, содержащую нули:
Кп Ki2 О
Ки О
К22 К23
К23 Кзг
(20.6)
В силу симметрии матрицы будут иметь размерности
К,,(1X1), КАХт), KiiimXm):-
К23 (m X (ЛГ - 1 - /п)), KaMN - 1 - m) X (iV - 1 - т)).
При исключении лодматрицы Ки изменяется только К22, так как нулевая подматрица не изменяет Кгз и Кзз- Количество операций исключения пропорционально т, а общее количество операций - величине
или приблизительно ЧгМт, где /Птах - максимальная величина половины ширины ленты.
На практике половина ширины ленты обычно меньше /ю размерности матрицы, и описанный способ использования свойства ленточности матрицы позволяет уменьшить число арифметических операций почти до 3% от числа выполняемых операций прн использовании метода, не учитывающего ленточного характера матрицы.
Другим преимуществом является экономное использование памяти машины, так как матрица может быть помещена в прямоугольный массив размерности X "max, как показано иа фиг. 20.1. Однако при этом общее число уравнений, которое может быть решено, ограничивается размерами прямоугольного массива.
/С24
55
Kll К12
К22 К23
Кзз Кз4
/Cg5 О
О 2i О
Фиг. 20.1. Хранение матрицы при решении уравнений методом ленточных
матриц.
Ленточная матрица, задаваемая формулой (20.6), особенно удобна при решении больших систем уравнений, так как в процессе исключения операции производятся только над К22 и поэтому только эта матрица нужна в оперативной памяти. Однако истинная эффективность достигается, если составление уравнений для ансамбля и исключение производятся параллельно, т. е. если уравнение исключается сразу после его составления. Это легко осуществить, если в качестве первого узла каждого элемента всегда использовать узел с наименьшим номером, а элементы располагать так, чтобы номера их первых узлов располагались последовательно. При таком способе процесс составления уравнений ансамбля будет автоматически прекращаться, как только номер первого узла станет больше номера рассматриваемого узла. После исключения Кц оставшаяся матрица сдвигается так, что первый элемеит измененной матрицы К22 занимает позицию (1,1), а остальные - соответствующие им места. Исключенные уравнения могут временно пересылаться в промежуточное запоминающее устройство, а затем переписываться в виде блока в накопитель до следующего вызова. Для небольших задач объем промежуточной памяти может оказаться достаточным для того, чтобы вместить в себя все исключенные уравнения, при этом делать пересылку нет никакой необходимости.
Новая матрица может быть вновь выражена через подматрицы Км, Ki2 и К22, после чего процесс исключения повторяется.
На каждом этапе исключения требуется (m+l)X(m + 2) .ячеек памяти или при наличии симметрии - половина этого
числа. С помощью этого метода можно решить практически неограниченное число уравнений при условии конечности максимальной ширины ленты. В большинстве случаев при решении за-дачиЧ1етодом конечных элементов точность ие является проблемой. При необходимости можно решать уравнения с двойной точностью, однако при этом увеличиваются используемый объем памяти и время решения. В некоторых случаях точность можно повысить, если подставить полученные значения в первоначальную систему уравнений и вычислить невязки, а затем, используя эти невязки в качестве новых правых частей, снова решить систему. Сумма этих двух решений даст уточненное решение. Однако чаще эти невязки используются для оценки ошибок округления.
20.5.2. Экономное распределение памяти для решения ленточных систем
Чтобы уменьшить используемый объем памяти вычислительной машины, в программе решения больших систем FESS (в настоящей главе ие приводится) для храиерия в памяти верхней
X-элементы, не хранящиеся вптттц
Фнг. 20.2. Порядок размещения в памяти уравнений в программе FESS.
треугольной части матрицы жесткости, необходимой при исключении самого верхнего уравнения, используется одномерный массив. Пример расположения в памяти элементов матрицы показан иа фиг. 20.2. В этом примере элементы матрицы располагаются по столбцам, так как это удобнее, чем построчное расположение. В программе FESS матрицы, необходимые для обратного хода [уравнение (20.5)], хранятся в верхней части одномерного массива. Если матрица жесткости и эти матрицы перекрывают друг друга, содержимое остатка переписывается иа ленту и процесс повторяется. Такой способ достаточно эффективен, так как позволяет использовать длинные записи иа ленту. В некоторых случаях удается решить уравнения без использования виешией памяти.
20.5.3. Итерационный метод Гаусса - Зейделя
В общем случае п-е уравнение системы N уравнений может быть записано в виде
а-1 N
1=\ i=n+l
Из этого уравнения можно найти
(20.7)
(20.8)
Если процесс итераций таков, что в правой части используются последние приближения б,-, то для те-й итерации имеем
:i==k-AF,-tknibT- т jn,6r-.
I 1 = 1 i=n-l )
(20.9)
В этом заключается итерационный процесс Гаусса - Зейделя.
Часто для уточнения решения используется прием, состоящий в умножении разности между двумя итерациями для б на некоторый коэффициент и представлении уточненной величины б в виде
б7 = бГ+р(бГ-бГ).
(20.10)
где б™* - величина, вычисленная в соответствии с (20.9), а р - коэффициент верхней релаксации, значение которого обычно лежит между 1 и 2. Установлено, что во многих практических случаях самым подходящим является значение, близкое к 1,8.
Итерационный метод Гаусса - Зейделя легко программируется. Матрица жесткости хранится в компактной форме без нулевых членов вместе с матрицей-указателем номеров столбцов, в которых находятся ее элементы. Каждое уравнение итерируется в соответствии с (20.9), и найденное значение уточняется в соответствии с (20.10). Процесс повторяется столько раз, сколько необходимо для получения приемлемого решения, причем сходимость обычно оценивается путем вычисления разности между двумя последовательными приближениями.
Итерационные методы удобны для решения нелинейных задач, так как при их использовании обычно требуется решение ряда сходных задач, поэтому для вектора решения всегда есть удачное начальное приближение. Кроме того, если на промежуточных этапах ограничиться невысокой точностью решения, то число арифметических операций уменьшается.
Недостатком использования итерационных методов является необходимость повторения основного цикла для всех уравнений.
Поэтому при использовании для хранения уравнений внешней памяти процесс решения занимает много времени. Однако компактное расположение матрицы жесткости дает возможность решать большое число уравнений с использованием только оперативной памяти, например 1000 дополнительных уравнений на 32К слов памяти. Главный же недостаток итерационных методов состоит в том, что при решении, как правило, может быть рассмотрен только один вектор нагрузки, так как перемещения занимают почти всю память.
В приложении 20А приведены две подпрограммы: FORMK (формирование матрицы жесткости и соответствующей матрицы-указателя) и SOLVE (решение системы уравнений итерационным методом Гаусса - Зейделя). Это упрощенные подпрограммы, в которых заданы число циклов и пределы сходимости. Кроме того, в них используется взаимно однозначное соответствие между матрицей жесткости и матрицей-указателем. Эти подпрограммы совместимы с остальными подпрограммами, приведенными в этой главе. Заметим, что подпрограмма FORMK вычисляет члены, лежащие как выше, так и ниже диагонали компактной матрицы жесткости, и этим она отличается от другой подпрограммы FORMK (программа 20-5), которая используется при применении прямого метода решения.
20.5.4. Некоторые другие прямые методы
Видоизменением метода исключения, используемого в случае ленточных матриц, является так называемый метод редкозапол-ненных матриц, часто применяемый в экономических расчетах. В соответствии с этим методом подматрица записывается в компактной форме вместе с матрицей-указателем, назначение которой состоит в указании номера строки и столбца полной матрицы. При решении задач методом конечных элементов с одинаковым числом степеней свободы в каждом узле матрица-указатель может быть записана в компактной форме для того, чтобы соответствовать узловым подматрицам.
В методе редкозаполненных матриц исключение производится точно так же,-как описано выше, но при этом в качестве индекса используется матрица-указатель. Преимущество, получаемое за счет избавления от некоторых операций, необходимо сравнить с дополнительными затратами времени на отыскание элемента в полной матрице с помощью матрицы-указателя. Требуемый объем памяти зависит от числа ненулевых подматриц на каждом шаге процесса исключения и не зависит от ширины ленты. Следует отметить, что при исключении на местах, где раньше стояли нули, появляются ненулевые подматрицы. Для максимально эффективного использования памяти и времени
