Лабораторная работа №2

РЕШЕНИЕ ЗАДАЧ ЛИНЕЙНОЙ АЛГЕБРЫ

1. Подключение пакета решения задач линейной алгебры linalg

> with(linalg);

2. Задание квадратной матрицы 3-го порядка

> A:=matrix([[2,1,3],[5,1,0],[7,8,9]]);

3. Вычисление определителя матрицы

> det(A);

4. Вычисление минора элемента с индексами i и j матрицы A

> minor(A,2,1);

5. Вычисление ранга матрицы

> rank(A);

6. Вычисление следа матрицы

> trace(A);

7. Транспонирование матрицы

> AT:=transpose(A);

8. Вычисление обратной матрицы

> A1:=inverse(A);

9. Умножение матриц (используется специальный символ &*)

> evalm(A&*A1);

То же самое можно было найти по команде

> multiply(A,A1);

10. Вычисление собственных значений матрицы

> eigenvals(A);

Из-за громоздкости результат выполнения этой команды не приводится; для представления в числовой форме рекомендуется использовать функцию evalf, например,

> eigenvals(A): evalf(");

12.23573596, -0.117867978 + 2.422913542 I,-0.117867978 –2.422913542 I

Другой вариант:

> evalf(eigenvals(A));

Собственные числа можно получить напрямую, используя их определение

> Id:=evalm(array(identity,1..3,1..3)); # задание единичной матрицы E

> B:=evalm(A-lambda*Id); # ввод матрицы A –l E

> f:=det(B); # вычисление определителя матрицы B

Это характеристический многочлен матрицы A; его также можно определить с помощью команды

> charpoly(A,lambda);

Корни этого многочлена определяются применением функции solve:

> res:=evalf(solve(f,lambda),4);

res:= 12.24, -0.118 + 2.422 I, -0.118 - 2.422 I

11. Нахождение собственных векторов

> evalf(eigenvects(A),4);

v:= [12.24, 1., {[2.247, 1., 7.340]}],

[-0.118 + 2.422 I, 1., {[-0.2237 + 0.4844 I, 1., -0.5665 – 0.5225 I ]}],

[-0.118 – 2.422 I, 1., {[-0.2237 – 0.4844 I, 1., -0.5665 + 0.5225 I ]}]

Формат вывода ясен из нижеследующего комментария:

12.24 – первое собственное число

1. – кратность этого числа

{[2.247, 1, 7.340]} – собственный вектор, соответствующий первому собственному числу

–0.118+2.422 I – второе собственное число и т. д.

12. Решение систем линейных алгебраических уравнений.

Пусть требуется решить систему уравнений

Ax = b,

где матрица A – введенная в п.2, а вектор правых частей b равен . Сначала осуществим ввод вектора b:

> b:=matrix(3,1,[5,6,1]);

Метод Крамера

> C:=evalm(A); # теперь C то же самое, что и матрица A

> C1:=copyinto(b,C,1,1); # замена в матрице C первого столбца на вектор b

> C:=evalm(A): C2:=copyinto(b,C,1,2);

> C:=evalm(A): C3:=copyinto(b,C,1,3);

Вычисление неизвестных как частных определителей соответствующих матриц:

> x1:=det(C1)/det(A);x2:=det(C2)/det(A);x3:=det(C3)/det(A);

> evalf([x1,x2,x3]);

Метод Гаусса

Гауссово исключение – приведение матрицы к треугольному виду

> H:=gausselim(A);

Если применить гауссово исключение к расширенной матрице системы, получим прямой ход метода Гаусса:

> Ar:=augment(A,b); # получение расширенной матрицы

> H:=gausselim(Ar);

Обратный ход метода Гаусса:

> backsub(H);

> evalf(");

Метод Гаусса–Жордана производит элементарные преобразования над расширенной матрицей так, что ее главный минор становится диагональной матрицей с единицами на главной диагонали, а последний столбец – решением системы.

> gaussjord(Ar);

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

> G:= matrix([[1,2,3],[2,1,3],[1,-1,0]]):

v:= vector([1,2,1]):

H1:= augment(G,v):

F1:= gaussjord(H1);

> backsub(F1);

Обратите внимание, здесь ранг равен 2, и одно неизвестное выбирается произвольно (обозначено ), а остальные выражаются через него.

С методом Гаусса тесно связано факторизация (разложение на сомножители) матрицы:

A = LU,

где L – нижняя треугольная матрица с единицами на главной диагонали, U – верхняя треугольная матрица. Следующая функция осуществляет такую факторизацию

> u:= LUdecomp(A,L='l');

Указанные треугольные матрицы хранятся в переменных u и l; вывести их можно с помощью стандартных команд

> evalm(l);

> evalm(u);

На основе полученного разложения можно получить решение системы Ax = b. Для этого последовательно решаются треугольные системы Lg = b и Ux = g. Для решения этих систем предусмотрены функции forwardsub и backsub:

> g:= forwardsub(l,b); # прямой ход

> x:= backsub(u,g); # обратный ход

> evalm(l&*u&*x - b); # проверка решения: вычисляется LUx – b

Симметричные положительно определенные матрицы факторизуются по методу Холесского, который по сути является одним из вариантов метода Гаусса. При этом разложение имеет вид

W = LL^T

(L – нижняя треугольная матрица). Ниже приведен пример на разложение Холесского матрицы AA^T, которая, как известно, при условии det A ¹0 обладает необходимым свойством.

> W:=multiply(A,AT); # произведение матрицы на свою транспонированную

> definite(W,positive_def); # проверка на положительную определенность

> L:=cholesky(W); # вычисление матрицы Холесского L

> evalm(L&*transpose(L)); # вычисление произведения LL^T для контроля

Решение с помощью встроенной функции linalg

Данная функция позволяет решать системы линейных уравнений вида Ax = b, а также матричных уравнений AX = B. Результат действия функции – вектор x или матрица X.

> linsolve(A,b); # решение линейной системы

Вектор правых частей здесь можно задавать как матрицу размером 3´1 (см. начало п. 12), так и непосредственно как вектор (vector):

> b:=vector([5,6,1]): linsolve(A,b);

13. Скалярное и векторное произведения векторов

> u:=vector([2,3,5]): v:=vector([-2.3,4,10]):

> dotprod(u,v); # скалярное произведение

> crossprod(u,v); # векторное произведение

> psi:=angle(u,v);evalf(psi); # угол между векторами в радианах

14. Норма вектора. Норма матрицы. Число обусловленности матрицы.

Норма – обобщенное понятие абсолютной величины числа. Норма может определяться по-разному (в соответствующей метрике). В качестве нормы вектора чаще понимается его длина. Норма матрицы может определяться по формулам

, , и т.д.

Примеры:

> norm(A,infinity);

> norm(A,1);

> norm(A,frobenius);

> norm(b,frobenius);

> norm(b);

Второй параметр в вызовах функции norm определяет способ вычисления нормы (по умолчанию действует параметр infinity).

Число обусловленности матрицы определяется как .

> cond(A);

> norm(A)*norm(A1);

> cond(A,frobenius);

Поиск по сайту: