Архитектура Аудит Военная наука Иностранные языки Медицина Металлургия Метрология Образование Политология Производство Психология Стандартизация Технологии |
Решение систем линейных уравнений итерационными методами.
Помимо прямых методов решения систем линейных алгебраических уравнений MATLAB обладает мощными современными итерационными методами решения больших систем уравнений. К итерационным методам относятся те, решение которых получается в ходе ряда шагов – итераций, постепенно ведущих к получению результата с заданной погрешностью. В качестве примера рассмотрим двунаправленный метод сопряженных градиентов, который реализуется с помощью функции bigcg.
Построим матрицу системы с три диагональной структурой и размером 10 x 10:
n = 10; on = ones(n, 1); A = spdiags([-2*on 4*on -on], -1: 1, n, n);
В качестве элементов вектора правой части возьмем сумму элементов матрицы A в соответствующих строках:
b = sum(A, 2);
Для решения системы Ax=b применим один из одиннадцати существующих форматов вызова функции bigcg: x = bicg(A, b); bicg converged at iteration 10 to a solution with relative residual 8.6e-014
В результате система выдает сообщение о сходимости итерационного процесса и приводит относительную погрешность.
Обратная матрица и определитель. Определитель квадратной матрицы находится с помощью функции det:
Det(A) ans = Обратная матрица находится с помощью функции inv:
Inv(A) ans = 3 -3 1 -3 5 -2 1 -2 1 Факторизация Холецкого. Если матрица системы является симметричной (или эрмитовой) и положительно определенной, то ее можно представить в виде произведения двух треугольных матриц: A = R’R, где R – верхняя треугольная матрица, R’ – транспонированная. Факторизация Холецкого выполняется с помощью функции chol:
R = chol(A) R = 1 1 1 0 1 2 0 0 1 Легко проверить, что произведение R’R дает исходную матрицу Паскаля. Факторизация Холецкого приводит систему Ax=u к виду R’Rx=u. Поскольку оператор ‘ \ ’ распознает системы с треугольными матрицами, приведенную систему можно решить очень быстро:
x = R\(R'\u) x = -12 LU факторизация. LU факторизация, или гауссово исключение, выражает любую квадратную матрицу A как произведение перестановки нижней треугольной матрицы L и верхней треугольной матрицы U, где матрица L имеет единичную главную диагональ. Перестановки важны как с теоретической, так и с практической точки зрения. Так, матрицу невозможно представить в виде произведения двух треугольных матриц без перестановки двух строк. С другой стороны, хотя матрицу можно представить в виде произведения двух треугольных матриц, при малом значении é элементы матриц-сомножителей будут очень большими, что приведет к большой численной погрешности. Для выполнения разложения служит функция lu:
· [ L, U ] = lu(X) возвращает верхнюю треугольную матрицу U и нижнюю треугольную матрицу L (точнее, произведение нижней треугольной матрицы и матрицы перестановок), так что X = L*U; · [ L, U, P] = l u(X) возвращает верхнюю треугольную матрицу U, нижнюю треугольную матрицу L и матрицу перестановок P, так что L*U = P*X; · B = lu(X) возвращает матрицу B такую что, B(i, j) = U(i, j) для всех индексов jSi и B(i, j) = L(i, j) для всех индексов j< i. Поскольку диагональные элементы матрицы L равны единице, их можно не хранить. Рассмотрим факторизацию магической матрицы B:
[L U] = lu(B) L = 1.0000 0 0 0.3750 0.5441 1.0000 0.5000 1.0000 0 U = 8.0000 1.0000 6.0000 0 8.5000 -1.0000 0 0 5.2941 Отметим, что строки матрицы L переставлены местами. Легко проверить, что B = L*U. LU факторизация матрицы B позволяет очень быстро решить систему Bx=u:
x = U\(L\u) x = 0.5528 0.2611 -0.2806
Получим представление матрицы перестановок P:
[L, U, P] = lu(B) L = 1.0000 0 0 0.5000 1.0000 0 0.3750 0.5441 1.0000 U = 8.0000 1.0000 6.0000 0 8.5000 -1.0000 0 0 5.2941 P = 1 0 0 0 0 1 0 1 0
Строки матрицы L стоят на своих местах. Легко проверить, что L*U = P*B. Определитель и обратную матрицу из LU разложения можно найти по формулам: det(B) = det(L)*det(U) и inv(B) = inv(U)*inv(L). QR факторизация. Функция qr выполняет QR разложение матрицы. Эта операция полезна как для квадратных, так и для прямоугольных матриц. Исходная матрица представляется в виде произведения действительной ортонормальной или комплексной унитарной матрицы Q и верхней треугольной матрицы R. Функция используется в следующих формах:
· [Q, R] = qr(A) вычисляет верхнюю треугольную матрицу R того же размера, что и A, и унитарную матрицу Q, так что A = Q*R; · [Q, R, E] = qr(A) вычисляет матрицу перестановок E, верхнюю треугольную матрицу R c убывающими по модулю диагональными элементами и унитарную матрицу Q, так что A*E = Q*R; · [Q, R] = qr(A, 0) и [Q, R, E] = qr(A, 0) вычисляют экономное разложение, в котором E – вектор перестановок такой, что Q*R = A(:, E). Вектор E выбирается так, чтобы диагональные элементы матрицы R убывали по модулю; · X = qr(A) возвращает результат из пакета LAPACK так, что R есть triu(X) – верхняя треугольная часть матрицы X. Матрица R позволяет избежать потери точности при вычислении A'*A с помощью разложения Холецкого: A'*A = R'*R.
В качестве примера рассмотрим QR разложение прямоугольной матрицы C:
C = fix(10*rand(3, 2)) C = 6 8 9 6 8 1 [Q, R] = qr(C) Q = -0.4460 0.7450 0.4961 -0.6690 0.0908 -0.7377 -0.5946 -0.6609 0.4579 R = -13.4536 -8.1762 0 5.8437 0 0 Во многих случаях последние m-n столбцы матрицы Q можно отбросить, так как они умножаются на последние нулевые строки матрицы R. Используя экономный формат функции разложения, получим: [Q, R] = qr(C, 0) Q = -0.4460 0.7450 -0.6690 0.0908 -0.5946 -0.6609 R = -13.4536 -8.1762 0 5.8437 Решение системы Cx = u с факторизованной матрицей находится в два этапа:
y = Q'*u; x = R\y x = 0.3590 -0.0544 Поскольку матрица системы является переопределенной, то найденное решение является наилучшим среднеквадратичным приближением. Отметим, что оператор x = C\u даст такой же результат.
Матричная экспонента. Матричная экспонента expm(X) возводит число e в матричную степень X. Рассмотрим в качестве примера задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами:
d x/ dt= Ax, x(0)=u,
здесь x = x (t) – векторная функция, A – матрица коэффициентов, u – вектор начальных условий. Решение задачи можно выразить через матричную экспоненту x( t )=e t A x(0). Положим:
A=[0 -6 -1; 6 2 -16; -5 20 -10]; u=[1; 1; 1];
Определим решение системы дифференциальных уравнений в узлах равномерной сетки с шагом 0.1 на отрезке [0, 1] и построим трехмерный график фазовой траектории:
X = []; for t = 0:.01: 1 X = [X expm(t*A)*u]; End plot3(X(1,: ), X(2,: ), X(3,: ), '-o') Box on
|
Последнее изменение этой страницы: 2017-03-17; Просмотров: 458; Нарушение авторского права страницы