Владимир Дьяконов - Maple 9.5/10 в математике, физике и образовании

M;

LUDecomposition: "calling external function"

LUDecomposition: "NAG" hw_f07adf

Конечной целью большинства матричных операций является решение систем линейных уравнений. Для этого пакет LinearAlgebra предлагает ряд методов и средств их реализации. Основными методами решения являются следующие:

• обращением матрицы коэффициентов уравнений и решением вида Х=А -1 *В;

• применением метода LU-декомпозиции (method='LU');

• применением метода QR-декомпозиции (method='QР');

• применением метода декомпозиция Холесского (method='Cholesky');

• метод обратной подстановки (method='subs').

Решение с применением обращения матрицы коэффициентов левой части системы уравнений А уже не раз рассматривалось и вполне очевидно. В связи с этим отметим особенности решения систем линейных уравнений другими методами. Любопытно отметить, что указание метода может быть сделано и без его заключения в одинарные кавычки.

Зададим матрицу А левой части системы уравнений и вектор свободных членов В:

> restart; with(LinearAlgebra): UseHardwareFloats := false:

> A:=<<4|.24|-.08>,<.09|3|-.15>,<.041-.08|4>>; B:=<8, 9, 20>;

Прямое решение этим методом выполняется одной из двух команд, отличающихся формой записи:

> х := LinearSolve(А, В, method= 'LU');

х := LinearSolve(<���А|B>, method='LU');

Проверим решение данной системы уравнений:

> А.х-В;

В данном случае решение точно (в пределах точности вычислений по умолчанию).

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

> P,L,U:=LUDecomposition(A);

> V2:=Transpose(Р).В;

> V3:=ForwardSubstitute(L,V2);

> x:=BackwardSubstitute(U,V3);

> A.x-B;

Выполним теперь решение для тех же исходных данных методом QR-декомпозиции, обозначив метод в функции LinearSolve:

> х := LinearSolve(А, В, method='QR');

> A.x-B;

Другой, более явный, но и более громоздкий метод решения представлен ниже:

> Q,R := QRDecomposition(А);

> V2:=Transpose(Q).B;

> x:=BackwardSubstitute(R,V2);

> A.x-B;

Тут, пожалуй, любопытно, что погрешность вычислений оказалась несколько выше, чем при использовании функции LinearSolve. Однако погрешность не выходит за рамки допустимой по умолчанию.