ЛИНЕЙНЫЕ УРАВНЕНИЯ И МЕТОД ИСКЛЮЧЕНИЯ ГАУССА

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

Пусть имеется система линейных уравнений

Ах = b, (2.4.1)

где А - матрица размера (nxn) с постоянными коэффициентами; b - n-мерный вектор известных констант и х — n-мерный вектор неизвестных:

= (2.4.2)

Формально эту систему уравнений можно решить, обратив матрицу А:

х = A b. (2.4.3)

Очень часто, когда требуется найти только одну "выходную" переменную, используют метод, называемый правилом Крамера. Это правило гласит, что для системы (2.4.1) к- якомпонента х_к вектора х равна отношению определителя матрицы А, в которой k -й столбец заменен вектором Ь, и определителя матрицы А:

det (матрицы А с к-м столбцом, замененным на b)

_X = -------------------------------------------------------. (2.4.4)

det (А)

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

Численное решение систем линейных уравнений часто базируется на методе исключения Гаусса. Это один из лучших среди известных алгоритмов, и он основывается на том факте, что сложение одного уравнения системы с другим, возможно умноженным на константу, не изменяет решения системы.

Рассмотрим матричное уравнение (2.4.2) и перепишем его в координатной форме, обозначив элементы b_i вектора b через a . Это упростит дальнейшие обозначения. Тогда система примет вид:

a x +a x +…+a x =a ,

a x +a x +…+a x =a , (2.4.5)

…………………………………..

a x +a x +…+a x =a ,

Разделим первое уравнение на а и запишем его в виде:

x +a x +a x +…=a ,

где введено обозначение a =a /a и т.д. Умножим это уравнение на a и сложим его со вторым уравнением. Коэффициенты вновь полученного второго уравнения a =a -a a , j = 1, 2,..., п + 1. Такой выбор множителя обеспечивает равенство нулю коэффициента a . Аналогично для других уравнений подстановка

a =a -a a i=2,3,…,n, j=1,2,…,n+1,

обеспечивает равенство нулю всех коэффициентов в первом столбце матрицы А, за исключением а который равен 1. Фактически не нужно вычислять элемент, который становится нулевым. Элементы а теперь не занимают память ЭВМ, и вычисления выполняются, начиная с j = 2. В результате этих операций уравнения примут вид

x +a x +a x +…+a x =a ,

a x +a x +…+a x =a ,

.……………………………………

a x +a x +…+a x =a ,

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

a =a /a ; a =a -a a ; i=3,4,…,n; j=3,4,…,n+1.

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

a =a /a ; a =a - a a ; k=1,2,…,n; i=k+1,…,n; j=k+1,…,n+1. (2.4.6)

В результате система уравнений приводится к виду

x +a x +a x +…+a x =a (2-4.7)

x +a x +…+a x =a

x +…+a x =a

…

x =a .

Для определения вектора неизвестных необходимо произвести обратную подстановку. Последняя составляющая х_п используется в (n - 1)-м уравнении для определения х и т. д. В общем виде обратную подстановку можно записать так:

x =a - , i=n-1, n-2,…,1. (2.4.8)

Здесь верхние индексы опущены.

В вычислительной технике принято оценивать эффективность алгоритмов числом операций, причем каждая операция представляет собой комбинацию умножения и вычитания. Можно показать, что исключение по Гауссу требует выполнения ~n³/3 операций, где п — порядок матрицы, а обратная подстановка может быть выполнена приблизительно за n²/2 операций.

2.5. МЕТОД LU-РАЗЛОЖЕНИЯ

Лучшим методом решения системы линейных алгебраических уравнений является метод разложения на треугольные матрицы, или метод LU-факто-ризации. Алгоритмы этого метода близки к методу исключения Гаусса, хотя вычисления могут проводиться в различной последовательности. Главным преимуществом метода LU-факторизации в сравнении с методом исключения Гаусса является возможность более простого получения решений для различных векторов b в правой части системы (2.4.1), а также для транспонированной системы уравнений, что требуется при расчете чувствительностей.

Допустим, что матрицу системы уравнений (2.4.1) можно разложить на два сомножителя:

A = LU, (2.5.1)

где

L= (2.5.2)

U= (2.5.3)

Здесь матрица L является нижней треугольной, а матрица U - верхней треугольной. Заметим, что на главной диагонали матрицы U стоят единицы. Это означает, что определитель матрицы А равен произведению диагональных элементов l_{{ i} матрицы L.

Рассмотрим алгоритм определения матриц L и U позже, а сейчас допустим, -что такое разложение осуществлено. Перепишем систему уравнений в следующем виде:

LUx=b. (2.5.4)

Определим вспомогательный вектор z как

Ux = z. (2.5.5)

Из этого уравнения вектор z найти невозможно, поскольку неизвестен вектор х. Однако, подставив z в (2.5.4), получим

Lz=b. (2.5.6)

Благодаря специальной форме матрицы L вектор z можно легко определить. Для этого запишем (2.5.6) в виде системы уравнений

l z =b

l z +l z =b

l z +l z +l z =b

……………………………………

l z +l z +…+l z =b

откуда получаем

z =b /l

z = (b -l z )/l

z = (b -l z -l z )/l ,

или в общем виде

z =b /l

z = (b - )/l , i=2,3,…,n.

Этот процесс называют прямым исключением (прямой подстановкой или прямым ходом). Чтобы уравнение (2.5.7) имело смысл, диагональные элементы матрицы L должны быть ненулевыми. Теперь вернемся к (2.5.5) и найдем вектор неизвестных х. Для этого запишем (2.5.5) в координатной форме

x +u x +u x +…+u x =z

x +u x +…+u x =z

………………………….

x +u , x =z

x =z .

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

x =z ,

x =z - , i=n-1,n-2,…,1. (2.5.8)

Этот процесс называют обратной подстановкой или обратным ходом. Число операций, требуемых для выполнения как прямой, так и обратной подстановок, равно примерно n²/2, а в сумме для решения требуется п² oneраций. Изучение уравнений (2.5.7) показывает, что компоненты b используются только для определения величин z и позднее не требуются. Аналогично в (2.5.8) величина z_i не нужна после вычисления переменных х . Следовательно, при такой системе расчетов векторы b, z, и х могут быть размещены в одних и тех же ячейках памяти ЭВМ. Обратите также вниманиена эквивалентность обратной подстановки в уравнениях (2.5.8) и (2.4.8).

Займемся выводом алгоритма LU-разложения. С этой целью рассмотрим

матрицу размера 4X4. Предполагая, что разложение существует, запишем произведение матриц L и U:

Сравним это произведение с матрицей А. Видно, что первый столбец разложения остается неизменным и 1 = a , i = 1, 2, 3,4. Заметим также, что первая строка произведения может быть использована для определения элементов первой строки матрицы U из решения уравнений l u =a ,j =2,3,4. Поскольку во втором столбце элементы и₁₂ и l известны, l =a -l u , i = 2, 3, 4. Так как теперь известны l ₂₁, 1 и и , можно использовать вторую строку матриц для расчета u ,j= 3,4:

u =(a -l u )/l j=3,4.

Далее, чередуя строки и столбцы, можно аналогичным образом найти остальные элементы матриц L и U.

Чтобы получить общие соотношения, запишем произвольный элемент произведения матриц L и U:

a = =

где верхний предел суммы учитывает наличие нулевых элементов в матрицах L и U. Рассмотрим произвольный элемент на или под главной диагональю матрицы А, для которого i > j, и заменим индекс / на индекс к. При этом, положив u = 1, получим

a = =l +

ИЛИ

l = a - i≥k.

(2-5.9)

Аналогичным образом, рассматривая произвольный элемент над главной диагональю, для которого i <j, и используя индекс к вместо i, находим

a = =l u +

После преобразования придем к следующему выражению для элементов матрицы U:

u =(a - )/l (2.5.10)

Уравнения (2.5.9) и (2.5.10) описывают алгоритм разложения на треугольные матрицы, называемый алгоритмом Краута. Его выполнение осуществляется при задании t=l,2,…,n использовании формул (2.5.9) и (2.5.10).

Заметим, что требуемые в этих соотношениях значения элементов матриц L и U рассчитываются на предыдущих этапах процесса. Далее, каждый элемент матрицы А требуется для вычисления только соответствующих элементов матриц L и U. Так как нулевые элементы матриц L и U, а также единичную диагональ матрицы U запоминать не нужно, в процессе вычислений матрицы L и U могут быть записаны на месте матрицы А, причем L расположена в нижнем треугольнике (i>j), a U — соответственно в верхнем треугольнике (i<j) матрицы А. Обобщив все вышесказанное, опишем алгоритм LU-разложения следующим образом:

Шаг 1. Положим к = 1 и перейдем к шагу 3.

Шаг 2. Используя (2.5.9), рассчитываем k- йстолбец матрицы L. Если к =п, закончим процедуру разложения.

Шаг 3. Используя (2.5.10), рассчитываем k- юстроку матрицы U.

Шаг 4. Положим k = k + 1 и перейдем к шагу 2

МАТРИЦА А И ВЕКТОР ПРАВОЙ ЧАСТИ:

4.0000 1.0000 1.0000 1.0000 1.0000

1.0000 4.0000 1.0000 1.0000 2.0000

1.0000 1.0000 4.0000 1.0000 3.0000

1.0000 1.0000 1.0000 4.0000 4.0000

ФАКТ0РИ30ВАННАЯ МАТРИЦА И РЕШЕНИЕ:

4.0000 0.2500 0.2500 0.2500 -0.1429

1.0000 3.7500 0.2000 0.2000 0.1905

1.0000 0.7500 3.6000 0.1667 0.5238

1.0000 0.7500 0.6000 3.5000 0.8571

Рис.2.5.1. Программы разложения матрицы на треугольные сомножители

Можно получить другую форму алгоритма разложения, если последовательно строку за строкой рассматривать произведение матриц L и U. Вначале определяем первую строку матрицы U. Затем находим 1_гг, но расчет остальных элементов второго столбца матрицы L пока откладываем. Потом определяем вторую строку матрицы U, а затем элементы l и l и третью строку матрицы U и т.д. При таком алгоритме можно рассматривать матрицу А построчно. Преимущество такого подхода особенно проявляется при матрицах большого размера. На рис. 2.5.1 приведен текст подпрограммы LUROW, написанной на языке Фортран и реализующей указанный алгоритм разложения матрицы. Аналогичный алгоритм можно построить при последовательном переборе столбцов матрицы.

Еще одну форму алгоритма можно получить, если рассмотреть элементы матрицы А, разложенной на L- и U-сомножители. Представим это разложение в виде одной матрицы, где на месте нижней треугольной матрицы записана матрица L, а на месте верхней треугольной матрицы записаны элементы матрицы U, отличные от нуля и единицы. Эту матрицу можно представить в виде

Первый столбец и строка этой матрицы рассчитываются, как и раньше. Далее заметим, что первое вычитание можно произвести сразу во всей матрице, поскольку все элементы l и u известны. После этого находим все элементы второго столбца, т.е. l_i2, i = 2, 3,..., п. Теперь, чтобы найти оставшиеся элементы второй строки, необходимо произвести деление на рассчитанное значение l . Эта процедура аналогична процедуре, выполненной вначале, за исключением того, что индекс теперь начинается с цифры 2. После этого в оставшихся элементах матрицы можно выполнить второе вычитание и т.д. Алгоритм разложения можно, следовательно, представить в виде:

Шаг 1. Положим k = 1

Шаг 2. l_ik = a_ik; i>k.

Шаг 3.u =a /l ; j>k. (2.5.12)

Шаг4.а =а -1 _ки_к1 i,f>k.,

Шаг 5. Если k=п, закончим процедуру разложения, в противном случае положим k= k + 1 и перейдем к шагу 2.

Описанный алгоритм LU-разложения эквивалентен; приведению матрицы к верхней треугольной форме с помощью процедуры Гаусса. Отличие состоит в том, что элементы под главной диагональю в треугольной матрице замещаются не нулями, а элементами матрицы L. Подпрограмма LUG на языке Фортран, реализующая этот алгоритм, также приведена на рис. 2.5.1.

Для машинных приложений метод Краута или разложение по строкам или столбцам являются более предпочтительными в сравнении с методом исключения Гаусса. Причина этого заключается в структуре внутреннего цикла и его реализации (метки 20, 40 и 80 в подпрограмме CROUT, метка 20 в подпрограмме LUROW и метка 10 в подпрограмме LUG). Подпрограммы CROUT и LUROW реже обращаются к массиву, что сокращает время обращения к памяти. Кроме того, внутренние циклы часто программируют на языке ассемблера, а в подпрограммах CROUT и LUROW возможно представление переменной Г с двойной точностью.

Важные черты метода LU-разложения состоят в следующем:

1.Легко вычисляется определитель матрицы А

detA= (2.5.13)

2.Элементы матриц L и U могут быть записаны на месте элементов матрицы А и занесены в те же ячейки памяти (запоминать единичные элементы на главной диагонали матрицы U нет необходимости).

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

4.Для уравнений с транспонированной матрицей A^fx = c решение находится при том же LU-разложении. Анализ таких систем уравнений необходим при расчете чувствительностей. Этот вопрос подробно обсуждается в гл. 6.

Число операций М, требуемых для проведения LU-разложения, определяется как

M= = -

С точки зрения объема вычислений метод LU-разложения вместе с прямой и обратной подстановками эквивалентен методу исключения Гаусса. Его преимущества заключаются в свойствах, указанных в третьем и четвертом пунктах.

Пример 2.5.1. Вычислим разложение следующей матрицы на треугольные сомножители:

Матрицы, получаемые на каждом шаге разложения при использовании метода Краута, построчного разложения и метода Гаусса, приведены ниже:

L =

U=

Если цепь состоит только из пассивных элементов и независимых источников, то она обычно описывается системой уравнений с симметричной матрицей. В таких случаях матрицу можно разложить на сомножители вида А=U^t DU, где U - верхняя треугольная матрица, a D - диагональная матрица При этом требуемые объем памяти и машинное время можно сократить примерно в 2 раза по сравнению с LU-разложением.

Поиск по сайту: