Обратите внимание, что Maple теперь выдает информационные сообщения о новых условиях реализации операции инвертирования матриц с вещественными элементами и, в частности, об использовании алгоритмов NAG и арифметики, встроенной в сопроцессор.
Следующий пример иллюстрирует создание двух случайных матриц M1 и М2 и затем их умножение:
> M1:=RandomMatrix(2,3); М2:=RandomMatrix(3,3);
Multiply(M1,M2,'inplace'); M1;
Параметр inplace в функции умножения обеспечивает помещение результата умножения матриц на место исходной матрицы М1 — излюбленный прием создателей быстрых матричных алгоритмов NAG. Поскольку матрицы M1 и М2 заданы как случайные, то при повторении этого примера результаты, естественно, будут иными, чем приведенные.
Другой пример иллюстрирует проведение хорошо известной операции LU-разложения над матрицей М, созданной функцией Matrix:
> M:=Matrix([[14,-8,1],[-11,-4,18],[3,12,19]], datatype=float);
LUDecomposition(М,output=['NAG'],inplace);
ipiv:=%[1];
M;
LUDecomposition: 'calling external function'
LUDecomposition: 'NAG' hw_f07adf
6.3.3. Методы решения систем линейных уравнений средствами пакета LinearAlgebra
Конечной целью большинства матричных операций является решение систем линейных уравнений. Для этого пакет LinearAlgebra предлагает ряд методов и средств их реализации. Основными методами решения являются следующие:
• обращением матрицы коэффициентов уравнений и решением вида Х=А-1*В;
• применением метода LU-декомпозиции (method='LU');
• применением метода QR-декомпозиции (method='QР');
• применением метода декомпозиция Холесского (method='Cholesky');
• метод обратной подстановки (method='subs').
Решение с применением обращения матрицы коэффициентов левой части системы уравнений А уже не раз рассматривалось и вполне очевидно. В связи с этим отметим особенности решения систем линейных уравнений другими методами. Любопытно отметить, что указание метода может быть сделано и без его заключения в одинарные кавычки.
6.3.4. Решение системы линейных уравнений методом LU-декомпозиции
Зададим матрицу А левой части системы уравнений и вектор свободных членов В:
> 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;
6.3.5. Решение системы линейных уравнений методом QR-декомпозиции
Выполним теперь решение для тех же исходных данных методом QR-декомпозиции, обозначив метод в функции LinearSolve:
> х := LinearSolve(А, В, method='QR');
> A.x-B;
Другой, более явный, но и более громоздкий метод решения представлен ниже:
> Q,R := QRDecomposition(А);
> V2:=Transpose(Q).B;
> x:=BackwardSubstitute(R,V2);
> A.x-B;
Тут, пожалуй, любопытно, что погрешность вычислений оказалась несколько выше, чем при использовании функции LinearSolve. Однако погрешность не выходит за рамки допустимой по умолчанию.
