Решение систем линейных уравнений в MATLAB

Решение систем линейных уравнений (СЛУ) относится к самой массовой области применения матричных методов. В этом разделе вы найдете ответы на вопросы, каким образом применяются указанные методы и какие дополнительные функции имеет система MATLAB для решения систем линейных уравнений.

Как известно, обычная СЛУ имеет вид:

а 11 1 , а 12 ,X 2 ..., а 1n n =b 1

Здесь а 11 , а, 2 ,..., а пп коэффициенты, образующие матрицу А, которые могут иметь действительные или комплексные значения, 1 , х 2 ,..., х п неизвестные, образующие вектор X, и b 1 , b 2 ,..., b п .свободные члены (действительные или комплексные), образующие вектор В. Эта система может быть представлена в матричном виде как АХ=В, где А — матрица коэффициентов уравнений, X — искомый вектор неизвестных и В — вектор свободных членов. В зависимости от вида матрицы А и ее характерных особенностей MATLAB позволяет реализовать различные методы решения.

Для реализации различных алгоритмов решения СЛУ и связанных с ними матричных операций применяются следующие операторы: +,-,*,/, \, *, ' . Как отмечалось ранее, MATLAB имеет два различных типа арифметических операций — поэлементные и для массивов (векторов и матриц) в целом. Матричные арифметические операции определяются правилами линейной алгебры.

Арифметические операции сложения и вычитания над массивами выполняются поэлементно. Знак точки «.» отличает операции над элементами массивов от матричных операций. Однако, поскольку операции сложения и вычитания одинаковы для матрицы и элементов массива, знаки «.+» и «.-» не используются. Рассмотрим другие операторы и выполняемые ими операции.

* — матричное умножение; С = А*В — линейное алгебраическое произведение матриц А и В:

 

 

Для случая нескалярных А и В число столбцов матрицы А должно равняться числу строк матрицы В. Скаляр может умножаться на матрицу любого размера.

/ — правое деление. Выражение Х=В/А дает решение ряда систем линейных уравнений АХ=В, где А — матрица размера тхпи В — матрица размера nxk;

\ — левое деление. Выражение Х=В\А дает решение ряда систем линейных уравнений ХА=В, где А — матрица размера тхп и В — матрица размера nxk. Если А — квадратная матрица, то А\В — примерно то же самое, что и inv(A)*B, в остальных случаях возможны варианты, отмеченные ниже.

Если А — матрица размера пхп, а В — вектор-столбец с п компонентами или матрица с несколькими подобными столбцами, тогда Х=А\В — решение уравнения АХ=В, которое находится хорошо известным методом исключения Гаусса.

Если А — матрица размера тхп и тхп, а В представляет собой вектор-столбец с m компонентами или матрицу с несколькими такими столбцами, тогда система оказывается недоопределенной или переопределенной и решается на основе минимизации второй нормы невязок.

Ранг k матрицы А находится на основе QR-разложения (урок 11) с выбором ведущего элемента. Полученное решение X будет иметь не больше чем k ненулевых компонентов на столбец. Если k<n, то решение, как правило, не будет совпадать с pinv(A)*B, которое имеет наименьшую норму невязок ||Х||.

возведение матрицы в степень. Х А р — это X в степени р, если р — скаляр. Если р — целое число, то степень матрицы вычисляется путем умножения X на себя р раз. Если р — целое отрицательное число, то X сначала инвертируется. Для других значений р вычисляются собственные значения и собственные векторы, так что если [V.D]=eig(X), то X*p=V*D. A p/V. Если X — скаляр и Р — матрица, то Х А Р — это скаляр X, возведенный в матричную степень Р. Если X и Р — матрицы, то Х Л Р становится некорректной операцией и система выдает сообщение об ошибке. Возможный вариант решения матричного уравнения АХ=В с применением оператора ^ можно представить как Х=В*А ^ -1.

' — транспонирование матрицы, то есть замена строк столбцами и наоборот. Например, А' — транспонированная матрица А. Для комплексных матриц транспонирование дополняется комплексным сопряжением. Транспонирование при решении СЛУ полезно, если в матрице А переставлены местами столбцы и строки.

При записи СЛУ в матричной форме необходимо следить за правильностью записи матрицы А и вектора В. Пример (в виде m-файла):






А-[2  1

0

1:


1    -3

2

4;


-5    0

-1

-7:


1    -6

2

6]:


В=[8  9

-5

0]:


Х1=В/А




Х2=В*А ^ -1




X3=B*inv(A)







Эта программа выдает результаты решения тремя способами:

X1 =

3.0000 -4.0000-1.00001.0000

Х2 =

3.0000 -4.0000-1.00001.0000

X3 =

3.0000 -4.0000-1.00001.0000

Как и следовало ожидать, результаты оказываются одинаковыми для всех трех методов. При решении систем линейных уравнений, особенно с разреженной матрицей коэффициентов, полезно применение функций colmmd (colamd), symmmd (symamd).

этого хватит.

 

Функции для решения систем линейных уравнений с ограничениями MATLAB

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

X =lscov(A,B.V) —возвращаетвекторXрешения СЛУ видаА*Х=В+е, гдее — вектор шумов, который имеет ковариационную матрицу V. Реализуется метод наименьших квадратов в присутствии шумов с известной ковариацией. Прямоугольная матрица А должна быть размера тхп, где т>п. При решении минимизируется следующее выражение: (AX-b)'*inv(V)*(AX-b). Решение имеет вид X=inv(A'* inv(V)*A)*A'*inv(V)*B. Алгоритм решения, однако, построен так, что операция инвертирования матрицы V не проводится;

[X.dX] = lscov(A,B,V) — возвращает также стандартную погрешность X, помещая ее в переменную dX. Статистическая формула для стандартной погрешности вычисления коэффициентов имеет следующий вид:

mse = B'*(inv(V)-inv(V)*A*inv(A'*inv(V)*A)*A'*inv(V))*B./(m-n)

dX = sqrt(diag(inv(A'*inv(V)*A)*mse))

X =isqnonneg(A,B) — решение СЛУ АХ=В методом наименьших квадратов с неотрицательными ограничениями. А — действительная прямоугольная матрица, В — действительный вектор. Вектор X содержит неотрицательные элементы X>=Q, где j = 1, 2,... п. Критерий: минимизация второй нормы вектора В=АХ;

X = Isqnonneg(A.B.X0) — решение СЛУ с явно заданным неотрицательным начальным значением X для итераций;

[X,W] = Isqnonneg (...) — вместе с решением возвращает вторую норму вектора остатков в квадрате;

[X.W.W1] = Isqnonneg(..) — вместе с решением возвращает вторую норму вектора остатков в квадрате и вектор остатков W1;

[X.W.Wl.exitflag] = Isqnonneg (...) — вместе с решением возвращает вторую норму вектора остатков в квадрате, вектор остатков W1 и флаг exi tflаg, равный 1, если решение сходится после заданного по умолчанию числа итераций, и 0 — в противном случае;

[X.'W.Wl.exitflag,output] = Isqnonneg(...) — дополнительно возвращает число итераций в output;

[X.W.Wl.exitflag,output,lambda] = lsqnonneg(...) — дополнительно возвращает вектор lambda, минимизирующий произведение lambda W1 на последнем шаге итераций решения, lambda (i) близко к нулю, когда X(i) положительно, lambda (i) отрицательно, когда Х(1) равно 0;

[X.W.Wl.exitflag,output,lambda] = lsqnonneg(A.В,ХО,options) — обычно используется, если при предыдущей попытке решения системы exitflag=l или если необходимо изменить заданный по умолчанию порог отбора по X - tol X. равный 10*max(size(A))*norm(A, l)*eps (произведению первой нормы матрицы, большего из измерений матрицы, машинной точности и 10). Также используется такая форма записи, как X=lsqnonneg(A,B,XO,options). Параметры options должны быть предварительно заданы при помощи функции optimset. Функция Isqnonneg использует только поля 'display' и 'tolX' структуры options. Поэтому наиболее часто используемая вместе с Isqnonneg форма записи функции — options= optimset С tol X'.tol value), где tolvalue — новое значение порога отбора по X.

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

» С=[4 3:12 5:3 12];b=[1.45,41];D=lsqnonneg(C.b')

D =

2.2242

2.6954

 

Точное решение, метод наименьших квадратов и сопряженных градиентов MATLAB

lsqr(A, В)—возвращает точное решение X СЛУ А*Х=В, если матрица последовательная, в противном случае — возвращает решение, полученное итерационным методом наименьших квадратов. Матрица коэффициентов А должна быть прямоугольной размера тхя, а вектор-столбец правых частей уравнений В должен иметь размер т. Условие m>=n может быть и необязательным. Функция 1 sqr начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размером п, состоящий из нулей. Итерации производятся или до сходимости к решению, или до появления ошибки, или до достижения максимального числа итераций (по умолчанию равного min(20, m, n) — либо 20, либо числу уравнений, либо числу неизвестных). Сходимость достигается, когда отношение вторых норм векторов norm(B-Ax)/norm(B) меньше или равно погрешности метода tol (по умолчанию 1е-б);

lsqr(A.B,tol) — возвращает решение с заданной погрешностью (порогом отбора) tol;

lsqr(A,b.tol .maxlt) — возвращает решение при заданном максимальном числе итераций maxit вместо, возможно, чересчур малого числа, заданного по умолчанию;

lsqr(A,b.tol .maxit,M) и lsqr(A,b,tol .maxit.Ml.M2) — при решении используются матрица предусловий М или М=М1*М2, так что производится решение системы inv(M)*A*x=inv(M)*b относительно х. Если Ml или М2 — пустые матрицы, то они рассматривается как единичные матрицы, что эквивалентно отсутствию входных условий вообще;

lsqr(A.B,tol .maxit.Ml.M2.X0) — точно задается начальное приближение Х0. Если Х0 — пустая матрица, то по умолчанию используется вектор, состоящий из нулей;

X = lsqr(A,B.tol .maxit,Ml.M2.X0) — при наличии единственного выходного параметра возвращает решение X. Если метод 1 sqr сходится, выводится соответствующее сообщение. Если метод не сходится после максимального числа итераций или по другой причине, на экран выдается относительный остаток попп(В-А*Х)/ norm(B) и номер итерации, на которой метод остановлен;

[X.flag] = lsqr(A.X.tol.maxit.Ml.M2.X0) — возвращает решение X и флаг flag. описывающий сходимость метода;

[X.flag.relres] = lsqr(A,X,tol.maxit,Ml.M2.X0) — также возвращает относительную вторую норму вектора остатков rel res=norm(B-A*X)/norm(B). Если флаг flag равен 0, то relres<tol;

[X.flag.relres.iter] = bicg(A,B.tol,maxit,Ml,M2.X0) — также возвращает номер итерации, на которой был вычислен X. Значение iter всегда удовлетворяет условию 0<iter<maxit;

[X.flag.relres,iter,resvec]= lsqr(A.B.tol.maxit.Ml.M2.X0) — также возвращает вектор вторых норм остатков resvec для каждой итерации начиная с res-vec(l)=norm(B=A*X0). Если флаг flag равен 0, то resvec имеет длину iter+1 и resvec(end)<tol*norm(B). Возможны значения flag, равные 0, 1, 2, 3 и 4. Значения flag предоставляют следующие данные о сходимости решения:

flag=0 - решение сходится при заданной точности tol  и числе итераций не более заданного maxit;

flag=l - число итераций равно заданному maxit, но сходимость не достигнута;

flag=2 - матрица предусловий М плохо обусловлена;

flag=3 - процедура решения остановлена, поскольку две последовательные оценки решения оказались одинаковыми;

fl ag=4 - одна из величин в процессе решения вышла за пределы допустимых величин чисел (разрядной сетки компьютера).

Если значение flag больше нуля, то возвращается не последнее решение, а то решение, которое имеет минимальное значение отношения вторых норм векторов norm(B-A*x)/norm(B).

Пример:

» А=[0 012; 1300; 0101; 1010];

» В=[11; 7; 6; 4];

Введенные в этом примере матрица А и вектор В будут использованы и в других примерах данного раздела. В примере процесс итераций сходится на пятом шаге с относительным остатком (отношением вторых норм векторов невязки и свободных членов) 1.9 10- 13 .

Пример:

» lsqr(A,B.1e-6.5)

Isqr converged at iteration 5 to a solution

with relative residual

1.9e-013

ans =

1.0000

2.0000

3.0000

4.0000