Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 21
Текст из файла (страница 21)
В нашем примере матрица У— это матрица системы (3.3.6), а в общем случае — матрица системы (3.3.9) . Перемножение матриц 1. и У показьвает, что А =АУ, жения выполняется за 1 мс, а операция умножения за 2 мс (мс = 10 с). -б Тогда время, затрачиваемое на операции сложения и умножения при тре- 1003 уголыюй факторизации, составит примерно — (3 мс) = 10бмс = 1 с. Конечно, понадобится еще время на выполнение других арифметических операций, входящих в полный алгоритм исключения, но оно будет значительно меньше, чем 1 с. Более существенными являются различные "накладные расходы", связанные, например, с пересылкой данных в память и из памяти.
Эти расходы могут легко удвоить или утроить полное время счета, но в любом случае на ЭВМ с таким быстродействием для решения системы размера 100 Х 100 потребуется всего несколько секунд. В приведенных выше рассуждениях предполагалось, что матрица системы является заполненной, т.е.
содержит мапо нулевых элементов. Многие матрицы, возникаюшие в практических задачах н, в частности, при решении дифференциальных уравнений, обладают, однако, тем свойством, что большинство их элементов равно нулю. Возможно„простейший нетривиальный пример матриц такого рода дают трехднагональные матрицы, рассьютренные в разделе 3.2; каждая строка такой матрицы содержит не более трех ненулевых элементов, независимо от величины и. Трехдиагональные матрицы представляют собой частный случай так называемых ленточных матриц, все отличные от нуля элементы которых располагаются на относительно небольшом числе диагоналей, примыкающих к главной диагонали. Структура такой матрицы показана на рис. 3.4.
Рнс. 3.4. Матрнцв е шнрнной ленты р+ о + 1 Матрица А = (а; ) является ленточной с шириной ленты р + д + 1, если = 0 для всех 1 н у таких, что ~' — / > р или >' — 1 > д. Все ненулевые элемсн1ы такой матрицы лежат на главной диагонали, р диагоналях ниже нее и ц янгоналях выше нее. Если р = ц = 1., то имеются только три ненулевых диагонали и матрица является трехднагональной. Трн различные трехдиагональные матрицы, приведенные в разделе 3.2, возникли в результате использования метода конечных разностей для решения линейных двухточечных краевых задач.
Аппроксимация производных разностями более высокого порядка приводит к матрицам с большей шириной ленты. Так, например, определяемая формулой (3.2З6) аппроксимация четвертого порядка ведет к матрице с шириной ленты, равной 5 (р = ц = 2) . Более широкие ленточные матрицы будут рассмотрены в гл. 8 в связи с численным решением уравнений в частных производных. 35 Давайте рассмотрим процесс гауссова исключения для трехдиагональной системы уравнений а! !х! + а! зхз = Ь|, я2 $х! + я22х2 + я2 зхз Ь2 (3.3.18) аззхз +!тззхз + аз4х! = Ьз а„„,х„! +а„„х„= Ь„.
Если, как и раньше, предположить, что а, ! и все последующие делители отличны от нуля, то легко видеть, что для трехдиагональных систем алгоритм гауссова исключения реализуется особенно просто. В первом столбце нужно исключить только один элемент а. !. и при этом в системе изменяются только коэффициенты а! з и Ь,. Полученная в результате первого шага преобразованная система снова является трехдиагональной: (!) ()) !ттрн хз + аз зхз = Ьз аззхз +аззхз +аз зх4 Ьз, 86 а„„!х„! а„„х„= Ь„, и относительно нее, как и относительно систем, получающихся на последую- щих шагах, справедливы те же самые замечания.
Так как каждый шаг преобразования требует только двух сложений, двух умножений и одного деления, а всего необходимо л — 1 таких шагов, то для приведения системы к треугольной форме требуется только 2(п — 1) сложений, 2(п — 1) умно- жений и и — 1 делений. Обратная подстановка реализуется тоже чрезвычайно просто. Каждое уравнение содержит самое большее два неизвестных, причем типичное уравнение имеет вид (! — ! ) — ~.(! — ! ) ал х;+а;;+! х;„, =Ь Таким образом, обратная подстановка требует только л — 1 сложений, и — 1 умножений и л делений, и мы приходим к следующему результату.
Число операций для трехдиагональной системы = 3(п — 1) сложений + + 3(л — 1) умножений+ (2п — 1) делений. (3.3.19) Сравнивая это число с 0(лз) операций, требуемых для заполненной матри- цы, видим, что для трехдиагональных систем алгоритм гауссова исключе- ния оказывается очень эффективным. Если, как и прежде, предположить, что операции умножения и сложения выполняются соответственно за 2 мс и 1 мс, а деление за 5 мс, то полное время выполнения арифметических операций для решения системы размера 100 Х 100 составит только (3) (99) (2 мс) + (3) (99) (мс) + (199) (5 мс) ( 0,002 с. Так как число операций растет линейно по п, то даже при решении систе- мы размера 10000 Х 10000 время выполнения арифметических операций на ЭВМ с такими характеристиками окажется меньше, чем 1 с.
Требования к памяти при решении трехдиагональных систем также ока- зываются минимальными. Естественно, мы не станем использовать для хранения матрицы системы двумерный массив размера л Х и, так как это 1 1г и1 а1г (3.3.20) 1и и„ 1 „ а„ ии (см. упражнение 3.3.б) и алгоритм разложения формулируется следующим образом: и, ч-а11.
Для1 =2,...,и 1~'-а~ ~ 1/и~ (3,3.21) и~'-ап — 1; ' а~ Матрицы вида (3.3.20) иногда называют бидиагональными. Многие соображения относительно трехдиагональных систем применимы и к ленточным системам общего вида. Пусть А — ленточная матрица, структура которой показана на рис. 3.4, причем для простоты будем предполагать, что а = р. В первом столбце подлежат исключению р элементов, и при этом исключении будут изменяться только те элементы матрицы А, которые заключены в строках и столбцах с номерами от 2 до р. Прн этом потребуется р сложений, рг умножений и р делений (не считая операций с правой частью системы) . Полученная в результате первого этапа преобразования матрица снова оказывается ленточной с той же самой шириной ленты, так что для ее преобразования требуется такое же число операций. После п — р — 1 таких шагов приходим к задаче преобразования заполненной матрицы размера (р + 1) Х (р + 1).
Таким образом, число сложений (илн умножений), необходимое для приведения ленточной матрицы к треугольной форме, определяется выражением (и — р — 1) р + (р/б) ( р т +-1) (2р+1) =пр -2р ~3. Хранение ленточной матрицы может быть организовано с помощью р + а + 1 одномерных массивов, алело гично случаю трехдиагональной матрицы. Однако, если р достаточно велико, по-видимому, лучше хранить диагонали матрицы А в виде столбцов двумерного массива размера п Х (р + а + 1), как это проиллюстрировано на рис, 3.5 прил = 5,р = 1 на = 2, Теперь вернемся к матрицам общего вида (не обязательно ленточным) . Во-- первых, отметим, что . определитель о а, ааа а„ Рис. 3.5. Хранение ленточной матрицы 87 бы означало выделение и ячеек памяти для хранения только Зп элементов, Для запоминания ненулевых элементов матрицы и вектора правых частей удобно использовать одномерные массивы.
Тогда, включая массив для решения, нам потребуется всего 5п ячеек памяти. Мы могли бы перефразировать описанный процесс исключения в терминах ( Уразложения. В данном случае матрицьс А, обозначаемый как деСА, легко получается в качестве побочного продукта процесса исключения. Действительно, так как определитель произведения двух матриц является произведением их определителей, а определитель треугольной матрицы равен произведению ее диагональных элементов, в результате Е У-разложения матрицы А имеем бес А = с1ес ЕУ = ссес Е с1ес 0 = и, сиз 2... и„„. (3.3.22) Таким образом, определитель является просто произведением диагональных элементов приведенной к треугольному виду системы и для его вычисления требуется только и — 1 дополнительных умножений. Отметим, что и в том случае, когда нас интересует только определитель матрицы, а не решение линейной системы, приведение к треугольной форме по методу гауссова исключения оказывается наилучшим общим способом его вычисления.
Процесс гауссова исключения является также и наилучшим способом обращения матрицы А. Пусть е; — вектор, с-й элемент которого равен 1, а остальные равны нулю. Тогда вектор ес представляет собой 1-й столбец единичной матрицы У и из определяющего обратную матрицу соотношения АА ' = Е следует, что 1-й столбец матрицы А ' является решением линейной системы уравнений Ах = е;. Таким образом, решая п систем уравнений Ах;=е;, 1=1,...,и, (3.3.23) и используя найденные векторы х„..., х„в качестве столбцов, получаем обратную матрицу А Можно рассмотреть более общую задачу, заключающуюся в решении нескольких систем с одной и той же матрицей коэффициентов Ахг=Ьь г' — 1,...,т.
(3.3.24) Эта задача может быть эффективно решена с помощью следующей модификации процесса (3.3.1б): 1. Факториэовать А = Е У. 2. Решить Еус = Ьь 1 = 1,..., т. (3.3.25) 3. Решить Ухг = у;, 1 = 1,..., пь Обратите внимание на то, что матрица А факторизуется только один раз, независимо от числа правых частей. Таким образом, число операций равно 0(п') + 0(ти'), где второе слагаемое соответствует шагам 2 и 3 в (3.3.25).