Главная
Новости
Строительство
Ремонт
Дизайн и интерьер




10.05.2022


09.05.2022


06.05.2022


04.05.2022


03.05.2022





Яндекс.Метрика





LUP-разложение

31.12.2021

LUP-разложение (LUP-декомпозиция) — представление данной матрицы A {displaystyle A} в виде произведения P A = L U {displaystyle PA=LU} где матрица L {displaystyle L} является нижнетреугольной с единицами на главной диагонали, U {displaystyle U} — верхнетреугольная общего вида, а P {displaystyle P} — т. н. матрица перестановок, получаемая из единичной матрицы путём перестановки строк или столбцов. Такое разложение можно осуществить для любой невырожденной матрицы. LUP-разложение используется для вычисления обратной матрицы по компактной схеме, вычисления решения системы линейных уравнений. По сравнению с алгоритмом LU-разложения алгоритм LUP-разложения может обрабатывать любые невырожденные матрицы и при этом обладает более высокой вычислительной устойчивостью.

Алгоритм LUP-разложения

Пусть A = ( a i j ) {displaystyle A=(a_{ij})} , L = ( l i j ) {displaystyle L=(l_{ij})} , U = ( u i j ) {displaystyle U=(u_{ij})} . На практике как правило вместо матрицы перестановок P используют вектор перестановок получаемый из вектора p i = ( 1 , 2 , . . . , n ) {displaystyle p_{i}=(1,2,...,n)} путём перестановки элементов соответствующих номерам строк переставляемых в матрице P. Например, если

P = ( 0 1 0 1 0 0 0 0 1 ) {displaystyle P={egin{pmatrix}0&1&01&0&0&0&1end{pmatrix}}}

то p = ( 2 , 1 , 3 ) {displaystyle p=(2,1,3)} так как матрица P получена путём перестановки первой и второй строки. Вычисление LUP-разложения ведётся в несколько шагов. Пусть матрица C = A. На каждом i-м шаге сначала производится поиск опорного (ведущего) элемента — максимального по модулю элемента среди элементов i-го столбца, находящихся не выше i-й строки, после чего строка с опорным элементом меняется местами с i-й строкой. Одновременно производится такой же обмен в матрице P. При этом, если матрица невырождена, то опорный элемент гарантированно будет отличен от нуля. После этого все элементы текущего i-го столбца, находящиеся ниже i-й строки, делятся на опорный. Далее из всех элементов c j k {displaystyle c_{jk}} находящихся ниже i-й строки и i-го столбца (то есть таких что j>i и k>i) вычитается произведение c j i c i k {displaystyle c_{ji}c_{ik}} . После этого счётчик i увеличивается на единицу и процесс повторяется пока i<n где n — размерность исходной матрицы. После того как все шаги будут выполнены матрица C будет представлять собой следующую сумму:

C = L + U − E {displaystyle C=L+U-E}

где E — единичная матрица.

В алгоритме используется три вложенных линейных цикла так что общую сложность алгоритма можно оценить как O(n³).

Реализация алгоритма на языке С++

Ниже представлен программный код приведённого выше алгоритма на языке С++. Здесь Matrix — некоторый контейнер, поддерживающий операцию индексирования. Обратите внимание, что отсчёт ведётся с нуля, а не с единицы.

void LUP(const Matrix &A, Matrix &C, Matrix &P) { //n - размерность исходной матрицы const int n = A.Rows(); C = A; //загружаем в матрицу P единичную матрицу P = IdentityMatrix(); for( int i = 0; i < n; i++ ) { //поиск опорного элемента double pivotValue = 0; int pivot = -1; for( int row = i; row < n; row++ ) { if( fabs(C[ row ][ i ]) > pivotValue ) { pivotValue = fabs(C[ row ][ i ]); pivot = row; } } if( pivotValue != 0 ) { //меняем местами i-ю строку и строку с опорным элементом P.SwapRows(pivot, i); C.SwapRows(pivot, i); for( int j = i+1; j < n; j++ ) { C[ j ][ i ] /= C[ i ][ i ]; for( int k = i+1; k < n; k++ ) C[ j ][ k ] -= C[ j ][ i ] * C[ i ][ k ]; } } } } //теперь матрица C = L + U - E