LU — разложение

Теоретическая часть

Пусть дано уравнение Ax=b, тогда LU-разложением называется разложение вида: LUx=b, где Ux=y, а Ly=b. При этом: LU=A, этим равенством можно пользоваться для проверки на правильность разложения.

Этот метод лежит в основе метода Гаусса.

Матрица U верхнетреугольная, по диагонали единицы, поэтому ее определитель равен единице. Матрица L — нижнетреугольная, ее определитель равен произведению элементов, расположенных по диагонали.

Прямой ход LU-разложения представляет собой прямой ход метода Гаусса, но в отличие от него, столбец b не участвует в треугольном разложении матрицы A, так как от него не зависят L и U.

Пусть A0=Матрица А0, тогда

каждый элемент матриц L  и U можно выразить следующими формулами:

n — размерность матрицы A (nxn)

Элементы по диагонали:

u{k,k}=1, (1)

l{k,k}=a{k-1}{k,k},где (k-1) — номер итерации, k=1…n, (2)

Остальные элементы:

l{i,j}=a{k-1}{i,j},(3)

u{j,i}=a{k-1}{j,i}/a{k-1}/{j,j}, где k=1…n, i, j=k+1…n. (4)

Обратный ход метода предполагает вычисление двух уравнений:

1.     Ly=b

2.     Ux=y

Из последнего уравнения получаем ответ.

Разберем первое уравнение:

y{1}=b{1}/l{1,1},y{2}=(b{2}-y{1}*l{2,1})/l{2,1}, значит, общая формула для поиска элементов столбца y выглядит следующим образом:

Формула поиска эле-тов столбца y

эта формула нам потом понадобится для вычисления y в программе, поэтому я убрал из общей формулы нахождение первого элемента y1, так как будут производиться лишние вычисления.

Разберем второе уравнение:

Здесь xn=yn, поэтому нет смысла включать его в общую формулу, по которой будут вычисляться все x, кроме последнего:

Общая формула поиска x{n-1}

Алгоритм:

1.     Вводим str/stlb — количество строк/столбцов, A — элементы расширенной матрицы, копируем A в A0;

2.     n:=1

3.     Если n<=str, тогда п.4, иначе п.11

4.     Формируем матрицы L[n,n]:=A0[n,n], U[n,n]:=1, сначала диагональные элементы, i:=n+1 (по формулам (1), (2))

5.     Если i<=str тогда п.6, иначе п.10

6.     Формируем остальные элементы: L[i,n]:=A0[i,n] и U[n,i]:=A0[n,i]/A0[n,n] (по формулам (3),(4)), присваиваем t:=A0[i,n]/A0[n,n], j:=n+1

7.     Если j<=str тогда п.8, иначе п.9

8.     A0[i,j]:=A0[i,j]-A0[n,j]*t, увеличиваем j на единицу, переходим к п.7

9.     i:=i+1 переходим к п.5

10.     n:=n+1 переходим к п.3

11.     Выводим матрицы L, U

12.     Считаем y1:=A[1,stlb]/L[1,1]; далее y по формуле (5)

13.     Выводим y;

14.     Считаем последний x[str]:=y[str]; далее по формуле (6)

15.     Выводим x;

16.     Проверки на невязку, заканчиваем алгоритм.

Описание входной информации: Str (Stlb) — количество строк (столбцов) в расширенной матрице, A [i, j] — матрица A (i — строки, j — столбцы).

Описание выходной информации: x — матрица-ответ.

Ссылки для скачивания:LU-разложение

Вернуться назад

7 комментариев к “LU — разложение”

Комментарии

  1. AlexR

    Денис, рад, что Вам помогла эта запись.

  2. Денис

    Спасибо. Помогли Ваши материалы разобраться с этим методом. Единственное, что смущает, так это формула для вычисления x. Там не нужно делить на u[i,j]! Проверял этот момент в литературе. В остальном все отлично.

  3. Георгий

    А кто сможет сделать блок схему по этому алгоритму, очень нужно? Я отблагодарю!

  4. Givi

    Помогите решить систему линейных уравнений методом LU-разложений в Scilab. Очень срочно

  5. AlexR

    Только что проверил на нескольких примерах. Мой метод работает безошибочно. Правда, я программу немного подправил для работы в Pascal ABC. NET.
    А что касается l[i,i], то здесь дело в том, что нет разницы (практически) ставить либо i, либо j-1…программа все-равно будет все считать верно.

  6. AlexR

    В ближайшее время проверю. Может, Вы пример какой-нибудь приведете, где алгоритм дает слабину?

  7. Максим

    не хочу портить малину, читал алгоритм, даже пробывал реализовать есть ошибки допустим в 5 формуле деление на l[i,i], а у вас в LU-разложении в коде программы на l[i,j-1]. Далее х не правильно находит, эт уже факт, L,U матрицы строит, в принципе че то похожее на правду, но по ним не правильно находится x.в 6 алгоритме деление u[i,j] не уместно, т.к судя по индексу j он будет j=i, u[i,i]=1. У меня есть своя задмука реализации, х считает правильно, но бывают погрешности большие, было бы не плохо списать и подправить их.

Комментарии