На рис. 6.1 показана графическая визуализация матрицы, полученной как разность матриц A и В. Для усиления эффекта восприятия применяется функциональная закраска разными цветами. Для задания цвета введена процедура F.
6.3. Работа с пакетом LinearAlgebra и алгоритмами NAG
6.3.1. Назначение и загрузка пакета LinearAlgebra
В новых реализациях систем Maple была сделана ставка на использование давно апробированных быстрых алгоритмов линейной алгебры, предложенных создателями Number Algorithm Group (NAG). Эти алгоритмы издавна применяются на больших ЭВМ и суперкомпьютерах, обеспечивая ускорение численных матричных операций от нескольких раз до нескольких десятков раз. Их применение обеспечивает эффективное использование систем символьной математики в решении задач, сводящихся к задачам линейной алгебры. В числе таких задач многочисленные задачи теоретической электротехники, механики многих объектов, моделирования электронных устройств и т.д.
В Maple 9.5/10 использование алгоритмов NAG реализуется пакетом LinearAlgebra. Для его загрузки используются следующие команды (файл NAG):
> restart; with(LinearAlgebra):
> infolevel[LinearAlgebra]:=1;
infolevelLinearAlgebra:= 1
Многие функции этого пакета (их большой список опущен) повторяет по назначению функции более старого пакета linalg, описанного выше. Поэтому мы не будем останавливаться на их повторном описании. Главное то, что эти функции задействуют возможности быстрых алгоритмов NAG и, в отличие от функций пакета linalg, ориентированы на численные расчеты в том формате обработки вещественных чисел, который характерен для применяемой компьютерной платформы. Знающий матричные методы читатель легко поймет назначение функций пакета LinearAlgebra по их составным названиям. Например, DeleteColumn означает удаление столбца матрицы, ToeplitzMatrix означает создание матрицы Теплица, ZeroMatrix — создание матрицы с нулевыми элементами и т.д. Все имена функций этого пакета начинаются с заглавной буквы.
6.3.2. Примеры матричных операций с применением пакета LinearAlgebra
Применение алгоритмов NAG особенно эффективно в том случае, когда используется встроенная в современные микропроцессоры арифметика чисел с плавающей запятой. С помощью специального флага такую арифметику можно отключать или включать:
> UseHardwareFloats := false; # use software floats
Use Hardware Floats := false
> UseHardwareFloats := true; # default behaviour
UseHardwareFloats := true
Матрицы в новом пакете линейной алгебры могут задаваться в угловых скобках, как показано ниже:
> М1:=<<1|2>,<4|5>>; М2:=<<1|2.>, <4|5>>;
После этого можно выполнять с ними типовые матричные операции. Например, можно инвертировать (обращать) матрицы:
> М1^(-1); М2^(-1);
MatrixInverse: "calling external function"
MatrixInverse: "NAG" hw_f07adf
MatrixInverse: "NAG" hw_f07ajf
Обратите внимание, что 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. Однако погрешность не выходит за рамки допустимой по умолчанию.
6.3.6. Решение системы линейных уравнений методом декомпозиции Холесски
Выполним решение еще и методом декомпозиции Холесски:
> x:=LinearSolve(А, В, method='Cholesky');
Приведем еще один пример решения системы из четырех линейных уравнений с применением метода декомпозиции Холесски:
> M_temp := Matrix(4, (i,j)->i+i*j-7, shape=triangular[lower]);
M :=M_temp.Transpose(M_temp);
IsMatrixShape(M, symmetric); IsDefinite(M);
> V := <6,1,3,-2>;
> x:=LinearSolve(M, V, method='Cholesky');
> M.x-V;
> M:=Matrix(3, (i,j)->i+2*j-8, shape=triangular[lower]); V:=<7,8,1>;
> x := ForwardSubstitute(M, V);
x := LinearSolve(M, V);
6.3.7. Одновременное решение нескольких систем уравнений
Мы ограничимся простым примером одновременного решения сразу трех систем уравнений. Дабы не загромождать книгу массивными выражениями, ограничимся решением систем из двух линейных уравнений, матрица коэффициентов у которых одна, а векторы свободных членов разные. Ниже показан пример решения такой системы:
> М:=Matrix([[1.,3],[4,5]],datatype=float);
V1:=<1.,2>;
V2:=<7,-11>;
V3:=<-34,-67>;
> LinearSolve(М,<V1|V2|V3>);
> М: =Matrix([[1.,3],[4,5]],datatype=float);
ipiv, M := LUDecomposition(M,output=['NAG'], inplace);
LinearSolve([ipiv, M], <V1|V2|V3>);
Ha этом мы завершаем обзор пакета LinearAlgebra. Читатель, познающий или знающий методы линейной алгебры, может опробовать в работе любые функции этого пакета самостоятельно или познакомиться с множеством примеров, размещенных в справочной системе Maple и в файле демонстрационных примеров LE_Linear_Solve.mws. Возможности пакетов linalg и LinearAlgebra удовлетворят самых требовательных специалистов в этой области математики.
