Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 12
Текст из файла (страница 12)
3.9. Алгоритмы скалярных произведений (слева) и столбцовый для обратнойподстановки (справа)значение xn , вычисляются и вычитаются из соответствующих элементов ciвеличины произведений xn uin (i=1,. . . ,n − 1); таким образом, вклад, вносимый xn в прочие компоненты решения, реализуется до перехода к следующему шагу. Шаг с номером i состоит из скалярного деления, сопровождаемого триадой длины i−1 (подразумевается,что шаги нумеруются в обратномпорядке: n, n − 1, . .
. , 2, 1). Какой из двух алгоритмов выбрать, диктуетсяспособом хранения матрицы U , если он был определен LU -разложением.Как алгоритм скалярных произведений, так и столбцовый алгоритм легкопереформулировать на случай нижнетреугольных систем, процесс решениякоторых называется алгоритмом прямой подстановки.3.7Задание на лабораторный проект № 2Написать и отладить программу, реализующую ваш вариант методаисключения с выбором главного элемента, для численного решения системлинейных алгебраических уравнений Ax = f , вычисления det A и A−1.Предусмотреть сообщения, предупреждающие о невозможности решенияуказанных задач с заданной матрицей A.Отделить следующие основные части программы:а) подпрограмму факторизации матрицы A, отвечающую вашему варианту метода исключения;б) подпрограмму решения систем линейных алгебраических уравнений;в) подпрограмму вычисления определителя матриц;г) подпрограммы обращения матриц;д) сервисные подпрограммы.693 Векторно-ориентированные алгоритмы LU-разложенияУделить особое внимание эффективности программы (в смысле экономии оперативной памяти).
Предусмотреть пошаговое выполнение алгоритмаисключения с выводом результата на экран.Выполнить следующие пункты задания.1. Провести подсчет фактического числа выполняемых операций умножения и деления при решении системы линейных алгебраических уравнений,сравнить его с оценочным числом (n3/3).2. Определить скорость решения задач (решение систем линейных алгебраических уравнений, обращение матриц) с учетом времени, затрачиваемогона разложение матрицы.
Для этого спроектировать и провести эксперимент,который охватывает матрицы порядка от 5 до 100 (через 5 порядков). Представить результаты в виде таблицы и графика зависимости времени выполнения (в минутах и секундах) от порядка матриц. Таблицу и график вывестина экран.3. Оценить точность решения систем линейных алгебраических уравнений, имеющих тот же самый порядок, что и задачи из пункта 2. Для этогосгенерировать случайные матрицы A, задать точное решение x∗ и образоватьправые части f = Ax∗.
Провести анализ точности вычисленного решения xот порядка матрицы. Результаты представить в виде таблицы и графика.Для заполнения матрицы A использовать случайные числа из диапазонаот −100 до 100. В качестве точного решения взять вектор x∗ = (1, 2, . . . , n)T ,где n — порядок матрицы. Для оценки точности использовать норму вектораkxk∞ = max(|xi|).i(3.12)4.
Повторить пункт 3 задания для плохо обусловленных матриц (см. подразд. 2.6), имеющих порядок от 4 до 40 с шагом 4.5. Вычислить матрицу A−1 следующими двумя способами.Способ 1 — через решение системы AX = I, где I — единичная матрица.Способ 2 — через разложение матрицы A в произведение элементарныхматриц, обращение которых осуществляется аналитически, а их произведение дает матрицу A−1.Сравнить затраты машинного времени и точность обращения матриц прииспользовании указанных способов 1 и 2. Эксперименты провести для случайных матриц порядков от 5 до 100 через 5. Для оценки точности в обоихспособах использовать оценочную формулу−1−1−1kA−1т − Aпр k ≤ kI − AAпр k · kAk .70(3.13)3.7 Задание на лабораторный проект № 2В этом выражении норму матрицы вычислять в соответствии со следующим определением:!nXkAk∞ = max|aij | ,(3.14)ij=1−1где A−1т — точное значение обратной матрицы, а Aпр — приближенное значение, полученное в результате обращения каждым из способов 1 и 2.6.
Провести подсчет фактического числа выполняемых операций умножения и деления при обращении матриц первым и вторым способами, сравнитьего с оценочным числом (n3).Замечание 3.3. По ходу проведения численных экспериментов наэкран дисплея должны выводиться следующие таблицы.Решение систем линейных алгебраических уравнений:ТеоретическоеРеальноеПорядокВремяТочностьчисло операций число операцийАналогичная таблица должна быть построена для плохо обусловленныхматриц.ПорядокОбращение матриц:ВремяТочностьЧисло операцийспос. 1 спос. 2 спос. 1 спос. 2 спос. 1 спос. 2 теорет.Замечание 3.4. Результаты экспериментов необходимо вывести наэкран в форме следующих графиков.Графики решения систем линейных алгебраических уравнений:••зависимость реального и расчетного числа операций от порядка матрицы (для разных графиков использовать разные цвета);зависимости времени и точности решения от порядка матриц.Графики для обращения матриц:••зависимость реального и оценочного числа операций от порядка матрицы (для разных графиков использовать разные цвета);зависимости времени и точности обращения первым и вторым способомот порядка матриц.713 Векторно-ориентированные алгоритмы LU-разложенияТаблица 3.5.
Варианты задания на лабораторный проект № 2ВидразложенияkjijkijikikjijkabcabcaabbA = L̄U12345678910A = LŪ11 12 1314 15 1617181920A = Ū L21 22 2324 25 2627282930A = U L̄31 32 3334 35 3637383940abc3.8kij— выбор ГЭ по столбцу активной подматрицы— выбор ГЭ по строке активной подматрицы— выбор ГЭ по активной подматрицеВарианты задания на лабораторный проект № 2В табл. 3.5 приведены 40 номеров вариантов задания на лабораторнуюработу (проект) № 2 с тремя стратегиями выбора главного элемента.Если нет других указаний преподавателя, выбирайте ваш вариант по вашему номеру в журнале студенческой группы.724Алгоритмы окаймленияв LU -разложении4.1Метод окаймленияХотя ijk-формы дают шесть различных способов организации LU -разложения, имеются и другие способы, потенциально полезные для векторных компьютеров.
Даже тогда, когда та или иная ijk-форма теоретическипригодна для конкретной векторной машины, при ее реализации могут возникнуть проблемы, особенно если применяется язык высокого уровня. Разбираемые ниже способы организации вычислений основаны на операциях сподматрицами; потенциально они проще реализуются и облегчают написание переносимых программ.В основе этих способов организации лежит идея окаймления [14].
Математически ее можно представить следующим образом. Разобьем матрицу Aна блоки и, соответственно, разобьем на блоки сомножители L̄ и U искомогоразложения L̄U = A: L̄11 0 0A11 a1j A13U11 u1j U13 lTj1 1 0 0ujj uTj3 = aTj1 ajj aTj3 .(4.1)A31 a3j A3300U33L31 l3j L̄33Здесь lTj1, uTj3 , aTj1 и aTj3 — векторы-строки, а l3j , u1j , a1j и a3j — векторыстолбцы; каждый из этих элементов находится на j-й позиции.4.2Окаймление известной части разложенияПусть известно разложение L̄11U11 = A11, которое можно рассматриватькак равенство (4.1) для блока A11.
Запишем аналогичные равенства для тех4 Алгоритмы окаймления в LU-разложениитрех блоков, которые окаймляют эту известную часть разложения, следуяправилу перемножения блок-матриц, а именно, для aTj1, a1j и ajj :lTj1 U11 = aTj1,L̄11u1j = a1j ,lTj1 u1j + ujj = ajj .(4.2)TИз первого уравнения (4.2), переписанного в виде U11lj1 = aj1 , находимlj1, из второго находим u1j и затем из третьего находим ujj = ajj − lTj1 u1j .При этом первое и второе уравнения описываются следующими нижнетреугольными системамиTU11lj1 = aj1 ,L̄11u1j = a1j .(4.3)Существуют два естественных варианта реализации окаймления известной части LU -разложения.В первом варианте треугольные системы (4.3) решаются с помощьюстолбцового алгоритма, во втором — с помощью алгоритма скалярных произведений.
Псевдокоды этих двух вариантов приведены на рис. 4.1.Для j = 2 до nДля k = 1 до j − 2Для i = k +1 до j −1aij = aij − lik akjДля k = 1 до j − 1ljk = ajk /akkДля i = k + 1 до jaji = aji − ljk akiДля j = 2 до nДля i = 2 до j − 1Для k = 1 до i − 1aij = aij − lik akjДля i = 2 до jlj,i−1 = aj,i−1/ai−1,i−1Для k = 1 до i − 1aji = aji − ljk akiРис. 4.1. Алгоритмы окаймления известной части L̄U-разложения: столбцовый (слева) и алгоритм скалярных произведений (справа)В первом цикле по i на рис. 4.1 (слева) выполняется модификация j-гостолбца матрицы A и тем самым вычисляется j-й столбец матрицы U .
Вовтором цикле по i модифицируется j-я строка матрицы A и вычисляется j-ястрока матрицы L̄. Заметим, что при i = j во втором цикле по i пересчитывается элемент (j, j) матрицы A; в результате, согласно (4.2), получаетсяэлемент ujj .Во второй форме алгоритма окаймления на рис. 4.1 (справа) первый циклпо i,k вычисляет j-й столбец матрицы U , для чего из элементов aij вычитаются скалярные произведения строк с 2-й по (j − 1)-ю матрицы L̄ с j-мстолбцом U .
Это эквивалентно решению системы L̄11u1j = a1j . Во втором744.2 Окаймление известной части разложенияцикле по i,k модифицируется j-я строка A путем делений (нормировок) элементов этой строки, сопровождаемых опять-таки вычитаниями скалярныхпроизведений j-й строки L̄ и столбцов U . Это эквивалентно решению треTугольной системы U11lj1 = aj1 относительно j-й строки матрицы L̄. Отметим, что здесь при j = i модифицируется элемент (j, j) матрицы A, и этоотносится уже к вычислению j-го столбца матрицы U ; в результате получается элемент ujj . Вычисления по этим формам показаны в табл. 4.1.Таблица 4.1. Вычисления по алгоритмам на рис. 4.1 для примера 3.3.
Позицииэлемента–делителя столбца L̄ показаны выделенным шрифтом21324485A−42106115j=22 4 −41/2223 812 50j=324 −41/2243/22/232506115611521/23/22/2j=44 −4624 −22/23 −61/22/34В обеих формах окаймления обращения к данным производятся одинаково, что показано на рис. 4.2.j⇓UL̄j⇒AРис. 4.2. Доступ к данным в алгоритмах окаймления известной части разложения.L̄, U — вычисление закончено, но обращения происходят. A — обращений не было.Вычисляются: j-й столбец матрицы U и j-я строка матрицы L̄Обратим внимание, что в обоих случаях требуется доступ и к строкам,и к столбцам матрицы A. Поэтому алгоритмы будут неэффективны длявекторных компьютеров, требующих, чтобы элементы вектора находилисьв смежных позициях памяти.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
