MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Вычисления и приближение данных в MATLAB

Решение систем линейных алгебраических уравнений и спектральных задач

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

Обзор возможностей MATLAB для решения систем линейных алгебраических уравнений

Среда MATLAB изначально разрабатывалась для работы с матрицами (MATLAB является сокращением от Matrix Laboratory), поэтому арсенал средств MATLAB для решения систем линейных алгебраических уравнений достаточно богат и включает в себя:

  • решение систем с квадратными и прямоугольными матрицами;
  • решение систем прямыми и итерационными (в том числе с возможностью предобусловливания) методами;
  • матричные разложения;
  • хранение больших разреженных матриц в компактной форме и специальные алгоритмы для решения систем с такими матрицами.

Решение систем при помощи знака обратной косой черты

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

Для этого заполняем матрицу и вектор-столбец правой части (правая часть должна быть именно столбцом, иначе выведется ошибка о несовпадении размерностей)

A = [4 1 2
        3 7 1
        2 2 8];
f = [7; 11; 12];

и используем знак обратной косой черты

x = A\f
x =
     1
     1
     1

Вместо знака обратной косой черты можно было вызвать функцию mldivide

x = mldivide(A, f)

Результат будет тем же самым

Вычисляя невязку, убеждаемся в том, что решение найдено верно (вообще говоря, не всегда малые значения компонент невязки свидетельствуют о верно найденном решении)

f - A*x
ans =
     0
     0
     0

Очень важно не перепутать местами A и f, т.к. при выполнении операции

x = f\A

никакой ошибки не будет, а выведется

x =
    0.2707    0.3439    0.3854

что не имеет ничего общего с решением рассматриваемой системы. Разберем, почему так происходит, записав систему с «матрицей» f и «правой частью» A. Неизвестная при этом будет только одна, т.к. f размера 3 на 1:

Если в правой части стоит матрица, то решатся столько систем, сколько столбцов у матрицы, при этом каждый столбец является вектором правой части соответствующей системы (т.е. оператор обратной косой черты позволяет решить сразу несколько систем). Итак, решаются системы

Строго говоря, решения ни у одной из этих систем нет. Однако, если матрица системы прямоугольная и число ее строк больше числа столбцов (т.е. число уравнений больше числа неизвестных), то такая система называется переопределенной (оverdetermined systems) и ее приближенное решение ищется минимизацией евклидовой нормы невязки. Так для первой системы решением будет x , доставляющий минимум следующему выражению

(4-7x1)2+(3-11x1)2+(2-12x1)2

Несложно убедиться, что минимум доставляет как раз . Воспользуемся, например, средствами Symbolic Math Toolbox

syms x
R = (7*x-4)^2 + (11*x-3)^2 + (12*x-2)^2
dRdx = diff(R)
x = solve(dRdx)
x =
85/314
double(x)
ans =
0.2707

Итак, полученный при записи x = f\A вектор [0.2707 0.3439 0.3854] содержит решение трех переопределенных систем, правая часть каждой из которых является соответствующим столбцом матрицы A.

Разумеется, матрица переопределенной системы не обязательно содержит только один столбец. В качестве примера рассмотрим задачу о подборе параметров a и b линейной модели по заданным табличным данным

k12345
xk1.01.52.02.53.0
yk2.992.812.893.033.21

Эта задача сводится к переопределенной системе относительно a и b, если мы потребуем равенства y(xk)=yk для k=1,2,...,5:

Составляем матрицу системы и вектор правой части

x = (1:0.5:3)'
A = [1./x log(x)];
f = [2.99;  2.81;  2.89;   3.03; 3.21];
и решаем ее 
c = A\f
получаем
c =
    2.9903
    2.0100
т.е.    и   . График данных и приближения подтверждает правильность полученных результатов
figure
fun = @(x)c(1)/x + c(2)*log(x);
fplot(fun, [1 3])
hold on 
plot(x, y, 'or')

Если в системе больше уравнений, чем неизвестных, то такая система называется недоопределенной (underdetermined system) и решается в MATLAB так же при помощи знака обратной косой черты. У нее будет бесконечно много решений и находится решение, содержащее как можно больше нулевых значений. Рассмотрим такой пример

A = [1 2 3 4
        5 5 7 8
        9 8 7 6];
f = [10; 25; 30];
x = A\f
x =
    1.3333
    1.0000
             0
    1.6667

Кроме знака обратной косой черты используется и знак обычной косой черты. Они связаны следующим правилом

B/A эквивалентно (A'\B')'
где апостроф означает транспонирование. Вместо знака обычной косой черты можно использовать также функцию mrdivide.

Пока мы не будем углубляться в алгоритмы решения систем, заложенные в операции обратной косой черты. Этому посвящен раздел Алгоритм решения системы линейных уравнений при помощи знака обратной косой черты. В частности, для решения систем с квадратной матрицей общего вида используется LU-разложение. Кроме того, проверяется обусловленность матрицы. Классическим примером плохообусловленной матрицы является матрица Гильберта, элементы которой определяются по формуле . Для создания матрицы Гильберта заданного размера в MATLAB служит функция hilb. Решим, например, систему 13-го порядка, правая часть которой такова, что ее решением должны быть все единицы (k-ый элемент правой части есть сумма элементов k-ой строки матрицы).
A = hilb(13);
f = sum(A, 2);
x = A\f

В результате получаем сообщение о плохой обусловленности

Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 2.339949e-018.

(RCOND есть оценка для величины, обратной числу обусловленности) и неверное решение

x =
    1.0000
    1.0000
    1.0004
    0.9929
    1.0636
    0.6552
    2.2023
   -1.7860
    5.3352
   -3.4773
    3.9431
   -0.1145
    1.1851

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

Матричные разложения Холецкого, LU и QR

В MATLAB есть функции для следующих разложений:

  • разложения Холецкого — функция chol;
  • LU-разложения — функция lu;
  • QR-разложения — функция qr.

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

Ax=f(k), k=1, 2,...,N.

Предположим, что матрица A представлена в виде произведения двух матриц, причем системы с матрицами A=BC причем системы с матрицами B и C решаются значительно быстрее, чем с матрицей A. Тогда решения исходных систем получаются следующим образом. Поскольку Ax=f(k) и A=BC, то BCx=f(k) и решая последовательно By=f(k) и Cx=y находим решение k-ой системы. Для заполненных матриц разложение выполняется за O(n3) операций, где n — размер матрицы (т.е. за то же время, что и решение системы), а решение системы с каждым множителем разложения за O(n2) операций. Поэтому решение N систем c с предварительным разложением занимает O(n3)+O(n2)N. Решение же каждой системы без разложения матрицы потребовало бы O(n3)N арифметических действий.

Разложение Холецкого для заданной матрицы A заключается в нахождении такой верхней треугольной матрицы R с положительными диагональными элементами, что A=RTR. Известно, что если матрица A симметрична и положительно определена (т.е. для любого вектораx верно: xTAx0, или, что тоже самое, все собственные числа матрицы положительны), то разложение Холецкого существует и при том единственно.

Рассмотрим в качестве примера разложение Холецкого матрицы

A = [4 1 1
        1 4 1
        1 1 4];
R = chol(A)
R =
    2.0000    0.5000    0.5000
             0    1.9365    0.3873
             0            0     1.8974

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

A - R'*R 
ans =
1.0e-015 *
         0         0         0
         0         0         0
         0         0         0.4441

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

A = [1 4 4
        4 1 4
        4 4 1];
R = chol(A)
??? Error using ==> chol
Matrix must be positive definite.

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

A = [4 1 2
        1 4 1
        1 1 4];
R = chol(A)
R =
    2.0000    0.5000    1.0000
         0        1.9365    0.2582
         0                0    1.7127
A - R'*R
ans =
     0     0     0
     0     0     0
    -1     0     0

Для произвольной квадратной матрицы можно выполнить LU-разложение, т.е. найти нижнюю треугольную матрицу L и верхнюю треугольную матрицу U такие, что A=LU. Более строгое утверждение формулируется следующим образом.

Если Ak — главный минор квадратной матрицы A размера n, составленный из первых k строк и столбцов (т.е. А(1:k, 1:k)) и det(Ak)0 для k=1,2,..., n-1, то существует единственная нижняя треугольная матрица L, диагональ которой состоит из единиц, и единственная верхняя треугольная матрица U такие, что A=LU.

При вычислении LULU-разложения может понадобится перестановка строк для обеспечения численной устойчивости процесса разложения, поэтому функция lu может вернуть матрицу L, которая является нижней треугольной с точностью до перестановки строк, например:

A = [0 1 1
        2 3 1
       1 2 5];
[L, U] = lu(A)
L =
             0    1.0000             0
    1.0000             0             0
    0.5000    0.5000    1.0000
U =
     2     3     1
     0     1     1
     0     0     4

Произведение LU где L — треугольная с точностью до перестановки строк, и будет равно матрице A (вообще говоря, с учетом погрешности при операциях с вещественными числами).

A-L*U
ans =
     0     0     0
     0     0     0
     0     0     0

Если для матрицы A сделано разложение Холецкого или LU -разложение, то решение системы с матрицей A, как уже было сказано выше, сводится к решению двух систем с треугольными матрицами. Это решение может быть выполнено при помощи операции обратной косой черты, поскольку заложенный в ней алгоритм определяет треугольные матрицы и применяет эффективный способ решения системы с ними, требующий O(n2) арифметических действий. Рассмотрим пример решения системы

с предварительным LU-разложением матрицы.

A = [4 1 2
        3 7 1
        2 2 8];
f = [7; 11; 12];

Выполняем LU-разложение

[L, U] = lu(A);

и решаем последовательно две системы с треугольными матрицами, сначала с L, затем с U

y = L\f;
x = U\y
x =
     1
     1
     1

Решение двух систем можно записать одним выражением

x = U\(L\f)

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

x = U\L\f

приведет к совсем другому результату, поскольку сначала выполняется U\L, что эквивалентно (как было разобрано выше) решению систем с матрицей U и столбцами L в качестве векторов правой части. В результате получится матрица, каждый столбец которой является решением соответствующей системы, затем с этой матрицей и вектором правой части f решится система. Этот процесс, очевидно, не имеет ничего общего с решением исходной системы.

Рассмотрим последнее матричное разложение — QR-разложение, выполняемое функцией qr. QR-разложение может быть проделано для прямоугольных матриц, именно если A матрица размера m на n, то существует ортогональная матрица Q размера m на m (т.е. такая, что Q-1=QT) и матрица R размера m на n со всеми нулевыми элементами под главной диагональю, что A=QR .

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

[Q, R] = qr(A)
Q =
   -0.7428    0.6000   -0.2971
   -0.5571   -0.8000   -0.2228
   -0.3714    0.0000    0.9285
R =
   -5.3852   -5.3852   -5.0138
             0   -5.0000    0.4000
             0         0        6.6108

После того, как сделано QR-разложение, решение системы Ax=f может быть найдено как решение системы с треугольной матрицей R, т.к. поскольку QR=f, то Rx=QTf:

>> x = R\(Q'*f)
x =
    1.0000
    1.0000
    1.0000

Более подробно про матричные разложения, соответствующие функции MATLAB chol, lu и qr и их аргументы написано в разделе Матричные разложения.

Алгоритм решения системы линейных уравнений при помощи знака обратной косой черты

При решении в MATLAB системы линейных алгебраических уравнений при помощи знака обратной косой черты
x = A\b
работает алгоритм, который в зависимости от типа матрицы решает систему при помощи наиболее подходящего метода, реализованного в одной из процедур пакета LAPACK или UMFPACK. Здесь b может быть и матрицей, число строк которой совпадает с числом строк матрицы A. Тогда возвращается матрица x, каждый i-ый ее столбец содержит решение системы Ax(i) = b(i), i=1,2,...,k, где b(i) — i-ый столбец матрицы b = [ b(1) | b(2)...| b(k) ].

Алгоритм, реализованный в операции \ состоит из следующих шагов:

  • 1. В тривиальном случае, если A разреженная и диагональная (разреженным матрицам в MATLAB посвящен раздел Разреженные матрицы), решение находится по простой формуле xk = bk/akk, где k=1,2,...n.
  • 2. Если A — квадратная, разреженная и ленточная матрица, то используется солвер для ленточных матриц. При этом могут быть два варианта:

    a. Если A — трехдиагональная матрица и b — один столбец из вещественных чисел, тогда система решается гауссовым исключением (его операции совершаются только с элементами на диагоналях). Если в процессе исключения для сохранения устойчивости потребуется сделать перестановки строк, или если матрица A не является трехдиагональной, то работает следующий пункт.
    b. Если A — ленточная с плотностью ненулевых элементов больше параметра bandden, по умолчанию равного 0.5, или не выполняются условия предыдущего пункта, тогда, в зависимости от типа A и b (double или single) вызываются следующие процедуры библиотеки LAPACK:

    A и b вещественные типа double — процедуры DGBTRF, DGBTRS
    A и b комплексные типа double — процедуры ZGBTRF, ZGBTRS
    A или b вещественные типа single — процедуры SGBTRF, SGBTRS
    A или b комплексные типа single — процедуры CGBTRF, CGBTRS

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

    nz = 1 (на диаг.№-2) + 6 (на диаг.№-1) + 7 (на диаг.№0) + 6 (на диаг.№1) + 1 (на диаг.№2) + 1 (на диаг.№3) = 22,

    а число всех элементов в ленте

    band = 5 (на диаг.№-2) + 6 (на диаг.№-1) + 7 (на диаг.№0) + 6 (на диаг.№1) + 5 (на диаг.№2) + 4 (на диаг.№3) = 33

    и плотность ненулевых элементов будет 2/3. Поскольку 2/3 > 0.5, то будет применен солвер для ленточных матриц
    Параметр bandden, как и другие параметры, управляющие алгоритмами для разреженных матриц в MATLAB, устанавливается при помощи функции spparms.

    Смысл данного шага алгоритма в том, что решение системы с ленточной матрицей не требует действий вне ленты. Если лента достаточно сильно заполнена, то действий с нулевыми элементами будет не очень много. Напротив, если лента достаточно сильно разрежена, то больший эффект могут принести подходы, приведенные в дальнейших шагах алгоритма.

  • 3. Если A — верхняя или нижняя треугольная матрица, то применяется метод обратной подстановки, или, соответственно, метод прямой подстановки, в котором из последнего (или первого уравнения) находится компонента решения и далее найденные компоненты подставляются в следующие уравнения для нахождения очередных компонент решения системы уравнений.
  • 4. Если A — может быть приведена перестановками к треугольной форме, то выполняется соответствующее приведение, а далее система уравнений решается как и в п.3.
  • 5. Если A — симметричная матрица (или, в комплексном случае, эрмитова), на ее диагонали находятся только положительные вещественные элементы, то в зависимости от того, разрежена A или нет, выполняется п. а) или п. б).

    a. Если А не является разреженной, тогда делается попытка выполнить разложение Холецкого A = RT R с последующим решением систем с треугольными матрицами RT и R: RTy = b и Rx = y. При этом вызываются следующие процедуры библиотеки LAPACK:

    A и b вещественные и типа double — DLANGE, DPOTRF, DPOTRS, DPOCON
    A и b комплексные и типа double — ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
    A и b вещественные и типа single — SLANGE, SPOTRF, SPOTRS, SDPOCON
    A и b комплексные и типа single — CLANGE, CPOTRF, CPOTRS, CPOCON

    Разложение Холецкого закончится неудачно, если матрица не является положительно определенной. Но предварительной проверки на положительную определенность не производится, а она определяется именно в процессе разложения, поскольку довольно часто достаточно выполнить несколько шагов разложения Холецкого для выяснения того факта, что матрица не положительно определенная (в процессе вычислений возникает необходимость извлечь квадратный корень из отрицательного числа). Если разложение Холецкого не выполнилось, то переходим к следующему пункту.

    b. Если А разрежена, то предварительно выполняются симметричные перестановки строк и столбцов по симметричному алгоритму минимальной степени (функцией symmmd) для уменьшения заполнения множителя разложения Холецкого, т.е. для уменьшения числа новых ненулевых элементов, возникающих в процессе заполнения:

    p = symmmd(A) — в вектор p записаны перестановки

    Далее производится разложение Холецкого матрицы со сделанными перестановками

    R = chol(A(p, p));

    и решаются две системы с множителями Холецкого, первая с транспонированным множителем и соответствующими перестановками в векторе правой части

    y = R' \ b(p);

    и вторая, с множителем разложения и с занесением компонент решения на соответствующие позиции в векторе x
    x(p, :) = R \ y

  • 6. Если A — хессенбергова матрица, то она приводится к верхней треугольной (с соответствующей модификацией правой части) и затем система решается подстановками
  • 7. Если A — квадратная матрица, не удовлетворяющая условиям пунктов 1-6, то в зависимости от того, разрежена она или нет, выполняются следующие действия

    a. Если A не является разреженной матрицей, то при помощи исключения Гаусса с выбором ведущего элемента (для обеспечения устойчивости разложения) выполняется LU-разложение матрицы A = LU, где

    U — верхняя треугольная матрица
    L — матрица, перестановками строк сводящаяся к треугольной

    и решение системы Ax = b находится последовательным решением систем с треугольными матрицами Ly = b, Ux = y.
    Для выполнения LU-разложения вызываются следующие процедуры библиотеки LAPACK:

    A и b вещественные и типа double — DLANGE, DGESV, DGECON
    A и b комплексные и типа double — ZLANGE, ZGESV, ZGECON
    A и b вещественные и типа single — SLANGE, SGESV, SGECON
    A и b комплексные и типа single — CLANGE, CGESV, CGECON

    b. Если A — разреженная матрица, то производится перестановка столбцов с целью уменьшения заполнения множителей L и U в процессе нахождения разложения. Далее перестановкой строк в процессе LU-разложения обеспечивается вычислительная устойчивость, после которого снова решаются системы с треугольными матрицами. Для выполнения LU-разложения разреженной матрицы используются процедуры библиотеки UMFPACK

  • 8. Если предыдущие пункты алгоритма не сработали и, следовательно, A не является квадратной матрицей, то все определяется ее разреженностью и работает один из приведенных ниже пунктов:

    a. Если A не является разреженной матрицей, то выполняется QR-разложение AP = QR, где P — матрица перестановки столбцов, Q — ортогональная матрица (QTQ = I), а R — верхняя треугольная. Если A размера m на n, то Q размера m на m, а R размера m на n. Далее решение системы находится так:

    x = P*(R \ (Q' * b))

    поскольку из Ax = b и AP = QR следует: QRx = bP и Rx = QTbP.

    b. В случае разреженной и прямоугольной матрицы A нет смысла вычислять матрицу Q, поскольку она скорее всего окажется заполненной. Поэтому в ходе алгоритма QR-разложения вычисляется c = QTb (т.е. модифицированная правая часть) и решается система Rx = c с треугольной матрицей для получения решения исходной системы Ax = b.

Кроме всего перечисленного, в приведенном выше алгоритме оценивается обусловленность матрицы и выводится предупреждение о возможной погрешности в решении, если матрица плохообусловлена (см. раздел Влияние обусловленности матрицы на точность решения системы с ней). Приведенный алгоритм содержит всевозможные проверки, которые могут занять ощутимое дополнительное время. В следующем разделе Функция linsolve для решения систем линейных уравнений приводится функция linsolve, которая позволяет задать тип матрицы, таким образом сократив время вычислений.

 

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

Вместо решения системы (или систем с несколькими правыми частями, заданными в матрице) линейных алгебраических уравнений при помощи знака обратной косой черты можно использовать функцию linsolve, которая не делает всех проверок матрицы, заложенных в алгоритме операции \ (см. предыдущий раздел Алгоритм решения системы линейных уравнений при помощи знака обратной косой черты).

Функция linsolve, вызванная в самом простом варианте с двумя входными аргументами и одним выходным аргументом
x = linsolve(A, b)
решает систему Ax = b одним из способов, в зависимости от того, квадратная матрица, или нет.

  • Если A — квадратная матрица, то предварительно вычисляется ее LU-разложение и затем решаются две системы с треугольными матрицами L и U.
  • Если A — прямоугольная матрица, то предварительно вычисляется ее QR-разложение и затем решается система с треугольными матрицами.

(Более подробно про действия с полученными множителями разложения написано в разделе Алгоритм решения системы линейных уравнений при помощи знака обратной косой черты). При этом, если матрица A плохообусловлена или A — матрица неполного ранга, то выводится соответствующее предупреждение, например:

A = [1 2 3
     2 4 6];
b = [6; 12];
x = linsolve(A,b)
Warning: Rank deficient, rank = 1,  tol =   4.4686e-015.
x =
         0
         0
    2.0000

Для подавления вывода данного сообщения в командное окно следует вызывать функцию linsolve с двумя выходными аргументами
[x, r] = linsolve(A, b)
тогда в x запишется полученное решение, а в r — или ранг матрицы (если A является прямоугольной матрицей), или обратная величина к оценке для ее числа обусловленности (см. раздел Влияние обусловленности матрицы на точность решения системы с ней), например

A = [1 2 3
     4 5 6 
     7 8 9+1e-10];
b = [1;1;1] ;
[x, r] = linsolve(A,b)
x =
 -1.0000e+000
  9.9999e-001
  3.3307e-006
r =
  6.9444e-013

Основные преимущества функции linsolve над операцией \ проявляются при указании типа матрицы системы уравнений. Для этого перед вызовом функции linsolve надо заполнить структуру с информацией о матрице со следующими полями

  • SYM — симметричная;
  • LT — нижняя треугольная;
  • UT — верхняя треугольная;
  • UHESS — хессенбергова;
  • POSDEF — симметричная, положительно определенная;
  • RECT — прямоугольная;
  • TRANSA — надо ли решать систему с матрицей, транспонированной к заданной.

Каждое поле может принимать значение либо true, либо false. Разумеется, не все комбинации значений полей допустимы, например матрица не может быть треугольной и в то же самое время симметричной и положительно определенной. Верные комбинации значений полей собраны в следующей таблице

LTUTUHESSSYMPOSDEFRECTTRANSA
truefalsefalsefalsefalsetrue/falsetrue/false
falsetruefalsefalsefalsetrue/falsetrue/false
falsefalsetruefalsefalsefalsetrue/false
falsefalsefalsetruetruefalsetrue/false
falsefalsefalsefalsefalsetrue/falsetrue/false

Если матрица системы положительно определена, то это обязательно учесть при решении, поскольку для положительно определенных матриц решение основано на разложении Холецкого, которое требует меньше операций, чем LU-разложение, используемое при решении систем с квадратными матрицами общего вида. В этом несложно убедиться при помощи следующего примера, в котором создается симметричная положительно определенная матрица (матрица из случайных чисел складывается с транспонированной к ней и на диагональ добавляются большие числа) и система уравнений с этой матрицей сначала решается как система с матрицей общего вида (opts.SYM = false и opts.POSDEF = false), а затем — как с симметричной и положительно определенной матрицей (opts.SYM = true и opts.POSDEF = true).

% задаем все поля структуры ops, кроме SYM и POSDEF
opts.TRANSA = false; opts.UHESS  = false;  opts.RECT  = false;
opts.UT = false; opts.LT = false;
% создаем вектор размеров матриц
N = 2.^(8:12);
% создаем пустые массивы для записи времен решения
TimeGeneral = []; TimePosSym = [];
% в цикле создаем матрицы и сравниваем времена решения
for n = N
    % создаем симметричную и положительно определенную матрицу
    % и вектор правой части
    A = rand(n); A = A + A' + 2*n*eye(n);   b = sum(A, 2);
    % решаем систему как систему с матрицей общего вида
    opts.SYM  = false; opts.POSDEF  = false; 
    Tstart = cputime;
    x = linsolve(A,b, opts); 
    Tend = cputime;  TimeGeneral = [TimeGeneral Tend-Tstart];
    % решаем систему как систему с симметр и полож. опред матрицей   
    opts.SYM  = true; opts.POSDEF  = true; 
    Tstart = cputime;
    x = linsolve(A,b, opts); 
    Tend = cputime;  TimePosSym = [TimePosSym Tend-Tstart];
end
% выводим графики временных затрат
figure
loglog(N, TimeGeneral, N, TimePosSym)
legend('TimeGeneral', 'TimePosSym')

Получающиеся вычислительные затраты этих способов решения системы приведены на графике ниже

Разумеется, если матрица треугольная, то очень важно это указать, поскольку решение системы с треугольной матрицей выполняется за O(n2) операций, а если к системе с треугольной матрицей применить алгоритм решения, предназначенный для матрицы общего вида, то понадобится порядка O(n3) операций.

Функция linsolve не проверяет, удовлетворяет ли матрица указанным свойствам, что может привести к ошибке. Например, если при решении системы со следующей матрицей и правой частью

A = [1 2 3
     2 4 5 
     3 5 6];
b = [6; 9; 6];
задать true для поля UT (а все остальные установить в false)
opts.UT = true; 
opts.TRANSA = false;
opts.LT = false;
opts.UHESS  = false;
opts.SYM  = false;
opts.POSDEF  = false;
opts.RECT  = false;
то при решении системы функция linsolve будет рассматривать матрицу как верхнюю треугольную и выберет нужные ей элементы из верхнего треугольника A
x = linsolve(A,b, opts)
При этом получится решение 
x =
     1
     1
     1
соответствующее матрице 
A = [1 2 3
     0 4 5 
     0 0 6];

 

Влияние обусловленности матрицы на точность решения системы с ней

В этом разделе мы рассмотрим, как ошибки в элементах вектора правой части и в матрице системы линейных уравнений Ax = b могут повлиять на решение этой системы. Обратимся сначала к случаю внесения возмущений в вектор правой части b. Итак, мы решаем две системы

причем предполагается, что каждую из систем мы решаем точно, а основной возникающий вопрос — насколько будет отличаться решение x системы (1) от решения системы (2) с возмущением в правой части δb. Это довольно важный вопрос, поскольку элементы правой части могут быть получены как результат некоторых измерений, соответственно содержащий погрешность δbk в каждой компоненте bk. Хотелось бы, чтобы при измерении одной и той же величины (каждый раз со своими, небольшими погрешностями) соответствующие решения систем с мало отличающимися правыми частями также не очень сильно отличались друг от друга. К сожалению, так бывает не всегда. Причиной этому — характеристика матрицы, называемая ее обусловленностью. Об этом и пойдет дальше речь.

Сначала требуется ввести меру близости векторов и x, т.е. вектора ошибки . Мерой величины вектора является норма (определяться она может по-разному). Возьмем пока в качестве нормы обычную евклидову норму вектора (квадратный корень из суммы квадратов его компонент), т.е.

 
, сответственно

где n — число неизвестных и уравнений в системе. Кроме векторной нормы для оценки величины вектора нам понадобится еще и матричная норма для оценки величины матрицы. Возьмем матричную норму, которая определяется следующим образом (она называется спектральной):

т.е. спектральная матричная норма равна квадратному корню из максимального собственного числа матрицы ATA. В MATLAB имеется функция norm для вычисления норм матриц и векторов, которая, в частности, умеет вычислять приведенные выше нормы. Почему мы выбрали именно такие нормы векторов и матриц, подробно объясняется в разделе Получение неравенства для ошибки, согласованная норма матрицы, где приведено немного выкладок и определений. Это связано с оценкой, которой мы будем пользоваться для ошибки в решении системы (вывод этого неравенства также приведен в упомянутом разделе):

Здесь через обозначено число обусловленности матрицы , которое определяется так:

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

Рассмотрим классический пример плохообусловленной матрицы — матрицу Гильберта — и выясним, как ошибка в правой части системы влияет на ошибку в решении. Матрица Гильберта определяется следующим образом

Для создания матрицы Гильберта в MATLAB предусмотрена функция hilb, входной аргумент которой задает размер матрицы. Возьмем небольшую матрицу 6 на 6:

n = 6;
H = hilb(n)
H =
    1.0000    0.5000    0.3333    0.2500    0.2000    0.1667
    0.5000    0.3333    0.2500    0.2000    0.1667    0.1429
    0.3333    0.2500    0.2000    0.1667    0.1429    0.1250
    0.2500    0.2000    0.1667    0.1429    0.1250    0.1111
    0.2000    0.1667    0.1429    0.1250    0.1111    0.1000
    0.1667    0.1429    0.1250    0.1111    0.1000    0.0909

Далее вычислим такую правую часть системы b, что все компоненты ее решения x равны единице:

x = ones(n, 1);
b = H*x
b =
    2.4500
    1.5929
    1.2179
    0.9956
    0.8456
    0.7365

Видим, что ни в матрице, ни в векторе правой части нет ничего «подозрительного», все числа не сильно отличаются друг от друга. Теперь сформируем возмущенную правую часть b + δb, добавив в вектор b небольшие числа порядка 10-5 и решим систему с возмущенной правой частью для получения вектора .

delta_b = 1e-5*(1:n)';
x_tilda = H\(b + delta_b)
x_tilda =
    0.9978
    1.0735
    0.4288
    2.6632
   -1.0160
    1.8593

Видно, что полученное решение далеко от точного, где должны быть все единицы. Вычислим относительную ошибку в решении и в правой части (функция norm по умолчанию вычисляет евклидову норму вектора):

delta_x = x - x_tilda;
LEFT =  norm(delta_x)/norm(x)
LEFT =
    1.1475
RIGHT = norm(delta_b)/norm(b)
RIGHT =
  2.7231e-005

Итак, ошибка в решении порядка единицы, хотя изменения в правой части были порядка 10-5. Это полностью укладывается в приведенное выше неравенство для ошибки. Действительно, вычислим число обусловленности cond(H) при помощи функции MATLAB, которая называется cond и по умолчанию вычисляет для спектральной нормы матрицы

c = cond(H)
c =
  1.4951e+007

Итак

= LEFT ≈ 1.1475 меньше
= c* RIGHT ? 1.4951e+07 * 2.7231e-05 ≈ 407


и неравенство выполнено (даже с некоторым запасом).

При увеличении размера матрицы Гильберта ошибка в решении будет только расти (это несложно проверить, задавая n = 7, 8, …). Причем, при n = 12 выведется сообщение о том, что матрица плохообусловлена и решение может оказаться неверным

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.409320e-017.

В качестве меры обусловленности здесь выбрана величина RCOND равная единице, деленной на оценку числа обусловленности (число обусловленности оценивается по достаточно быстрому алгоритму, который работает намного быстрее, чем cond, о чем более подробно написано в справке по функции rcond). Если значение RCOND мало, то матрица считается плохообусловленной.

Увеличение ошибки в решении при увеличении размера матрицы Гильберта происходит из-за того, что число обусловленности матрицы Гильберта растет очень быстро с ее размером. В этом несложно убедиться при помощи простого цикла и функции semilogy (масштаб по оси ординат — логарифмический):

N = 1:20; C = zeros(1, 20);
for n = N
    H = hilb(n);
    C(n) = cond(H);
end
semilogy(N, C)
grid on, title('cond(H)'), xlabel('n')

Кстати, поскольку функция cond находит число обусловленности при помощи численного метода (а именно сингулярного разложения для нахождения сингулярных чисел), то число обусловленности после n = 12 вычисляется уже неверно, на самом деле оно должно расти и дальше, в чем можно убедиться при помощи символьных вычислений в MATLAB и операций с заданным числом значащих цифр

N = 1:20; C = zeros(1, N);
digits(60)
for n = N
    H = vpa(sym(hilb(n))); % вычисление матрицы Гильберта с точностью до 60-го знака
    sigma = svd(H);            % нахождение сингулярных чисел матрицы Гильберта
    % вычисление числа обусловленности матрицы Гильберта
    C(n) = max(double(sigma))/min(double(sigma));
end
semilogy(N, C)
grid on, title('cond(H)'), xlabel('n')

Рассмотрим теперь три важных момента, касающихся данного примера.

Первый из них — решение системы Hx = b (с вектором правой части, соответствующим решению из единиц) при помощи оператора обратной косой черты не даст точные единицы, погрешность будет уже в десятом знаке (хотя в MATLAB все вычисления по умолчанию выполняются с двойной точностью)

format long
H\b
ans =
   0.99999999999916
   1.00000000002391
   0.99999999983793
   1.00000000042209
   0.99999999953366
   1.00000000018389

Так происходит из-за того, что при вычислении вектора b при умножении матрицы H на вектор из всех единиц мы уже заложили в него некоторую погрешность. Кроме того, ошибки округления в процессе решения системы также сыграли свою роль и плохая обусловленность матрицы (даже небольшого размера) привела к таким ошибкам в решении. Действительно, для матриц небольшого размера с небольшим числом обусловленности такого эффекта наблюдаться не будет. Возьмем матрицу 6 на 6 из случайных чисел, вычислим ее число обусловленности

A = rand(n);
cond(A)
ans =
  57.35245279907571

Затем сформируем правую часть, соответствующую точному решению из всех единиц

x = ones(n, 1);
b = A*x;  

и решим систему, получив в результате хорошую точность

>> A\b
ans =
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000

Второй важный момент связан с определителем матрицы. Вычислим определитель матрицы Гильберта размера 6 на 6, с которой мы решали систему

>> det(hilb(6))
ans =
  5.3673e-018

Определитель оказывается достаточно мал, откуда можно сделать неверный вывод, что источник проблемы большой ошибки в решении при небольшом возмущении правой части системы — малость определителя. В действительности это не так, если определитель системы мал, то это не означает плохой обусловленности матрицы и, следовательно, большой ошибки в решении. Действительно, возьмем матрицу

A = [1e-40 1e-42 1e-42
     1e-42 1e-40 1e-42
     1e-42 1e-42 1e-40];

Ее определитель очень мал, порядка 10-121

det(A)
ans =
  9.9970e-121

(слово «очень», конечно условно, но он меньше определителя матрицы Гильберта размера 6 на 6, при решении системы с которой у нас возникли проблемы). Однако, малость определителя никак не повлияет на ошибку в решении при возмущении правой части системы, что несложно показать, сформировав систему с известным решением, внеся возмущения в правую часть и решив систему:

x = [1e40; 1e40; 1e40];
b = A*x;
delta_b = 1e-5*[1; 2; 3]; 
x_tilda = A\(b+delta_b);
delta_x = x - x_tilda;
LEFT =  norm(delta_x)/norm(x)
RIGHT = norm(delta_b)/norm(b)
LEFT =
  2.1272e-005
RIGHT =
  2.1179e-005

Итак, относительная ошибка в решении практически равна относительной ошибке в правой части, поскольку число обусловленности приведенной выше матрицы A немного больше 1, а именно:

c = cond(A)
c =
    1.0303

И, наконец, рассмотрим третий вопрос, касающийся достижения знака равенства в неравенстве для ошибки в решении

Т.е., говоря другими словами, выясним: может ли быть такая ситуация, когда мы делаем небольшое относительное возмущение правой части системы, скажем 10-5, число обусловленности матрицы системы равно 1010, а в решении получается относительная ошибка 105. Можно проделать массу примеров, перебирая различные матрицы, не получить достижения равенства и ошибочно заключить, что это всего лишь завышенная оценка сверху для ошибки в решении. Однако, это не так, в чем нас убеждает следующий пример

A = [2.1005632211074903e+009  2.4405728346208658e+009  2.6288320556007657e+009  3.9543845780862803e+009
     5.6384551998346877e+008  6.5511289861702812e+008  7.0564654322518432e+008  1.0614591390422670e+009
     1.9531852931707768e+009  2.2693394431935925e+009  2.4443901805698233e+009  3.6769404175323849e+009
     2.2044117415857959e+009  2.5612308917984776e+009  2.7587973499916158e+009  4.1498830928366508e+009]; 
b = [ -5.7373057243726833e-001
      -1.5400413072907607e-001
      -5.3347697688693385e-001
      -6.0209490373259589e-001];
delta_b = [-0.71685839091451e-5
            0.54786619630709e-5
            0.37746931527138e-5
            0.20850322383081e-5];

в котором относительное возмущение правой части равно 10-5

RIGHT = norm(delta_b)/norm(b)
RIGHT =
  1.0000e-005

Относительная ошибка в решении системы оказывается 105

x = A\b;
x_tilda = A\(b+delta_b);
delta_x = x - x_tilda;
LEFT =  norm(delta_x)/norm(x)
LEFT =
  1.0000e+005

И происходит так потому, что
1) в этом примере число обусловленности матрицы A равно 1010;

c = cond(A)
c =
  1.0000e+010

2) неравенство для ошибки в решении переходит в равенство.

Если установить формат long, то мы увидим, что LEFT, RIGHT и число обусловленности не в точности 105, 10-5 и 1010, соответственно, но это вызвано ошибками округлений. При решении в точной арифметике, равенство получилось бы точным, что можно показать аналитически (см. раздел Получение неравенства для ошибки, согласованная норма матрицы, подробнее про norm и cond).

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

Ax = b

мы решаем систему

В этом случае мы так же можем оценить погрешность в решении, и в оценку снова входит число обусловленности матрицы

В самом общем случае ошибки могут быть как в матрице, так и в векторе правой части, и соответствующая оценка погрешности в решении имеет более сложный вид

 

Получение неравенства для ошибки, согласованная норма матрицы, подробнее про norm и cond

В этом разделе рассказывается:

  1. почему в оценках, приведенных в разделе Влияние обусловленности матрицы на точность решения системы с ней, нельзя использовать любые матричные и векторные нормы (вычисляемые функцией norm) и как они должны быть связаны между собой;
  2. как получается неравенство для ошибки в решении при возмущении правой части;
  3. почему число обусловленности матрицы не может быть меньше единицы;
  4. откуда берутся системы, для которых в оценке выполняется точное равенство (наихудший случай).

Итак, мы рассматриваем две системы: исходную (с точным решением x)

Ax = b

и вторую систему (с точным решением ), в правую часть которой внесены возмущения

Продемонстрируем вывод оценки

Из того, что Ax = b следует, что || Ax ||V = || b ||V (здесь введено обозначение || ||V для векторной нормы, и пока не так важно, какая именно выбрана). Далее требуется совершить такой переход

|| Ax ||V = || b ||V ⇒ || A ||M || x ||V ≥ || b ||V

где || ||M — некоторая матричная норма. Для того, чтобы это можно было сделать, матричная норма || ||M и векторная норма || ||V должны удовлетворять следующему неравенству

|| A ||M || x ||V ≥ || Ax ||V

для любых матриц A и векторов x. Если это неравенство выполняется, то матричная норма || ||M называется согласованной с векторной нормой || ||V. Известно, например, что спектральная норма

(квадратный корень из максимального собственного числа матрицы AT) согласована с евклидовой векторной нормой

Именно поэтому в разделе Влияние обусловленности матрицы на точность решения системы с ней мы проводили все эксперименты с привлечением данных норм.

Получив неравенство || A ||M || x ||V ≥ || b ||V далее заметим, что из Ax = b и следует, что . Поскольку матрица A невырожденная, то отсюда следует, что δx = A-1δb и || δx ||V = || A-1δb ||V. Опять используем согласованность норм и получаем неравенство || δx ||V ≤ || A-1 || M || δb || V. Далее в двух полученных неравенствах

|| A ||M|| x ||V ≥ || b ||V и || δx ||V ≤ || A-1 || M || δb ||V

меньшую величину одного из неравенств делим на большую величину другого неравенства, а большую, соответственно, на меньшую

и простым преобразованием окончательно получаем требуемое неравенство

где cond(A) = || A ||M * || A-1 ||M.

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

согласована с первой векторной нормой

а максимальная строчная матричная норма, вычисляемая по формуле

согласована с бесконечной векторной нормой

Функция norm вычисляет не только евклидову норму вектора и спектральную норму матрицы, но и перечисленные выше векторные и матричные нормы. Для этого ее требуется вызвать с дополнительным вторым аргументом:

  • q = norm(A, 1) — максимальная столбцевая норма матрицы A
  • q = norm(A, inf) — максимальная строчная норма матрицы A
  • q = norm(a, 1) — первая векторная норма a
  • q = norm(a, inf) — бесконечная векторная норма a

Число обусловленности матрицы cond(A) = || A ||M * || A-1 ||M относительно различных матричных норм может быть вычислено при помощи функции cond. Если cond вызывается с одним входным аргументом (матрицей), то вычисляется число обусловленности относительно спектральной матричной нормы. При вызове функции cond с дополнительным аргументом возвращаются числа обусловленности относительно указанной матричной нормы:

  • с = cond(A, 1) — число обусловленности относительно максимальной столбцевой нормы матрицы
  • с = cond(A, inf) — число обусловленности относительно максимальной строчной нормы матрицы

В качестве примера, демонстрирующего важность согласованности норы матрицы с нормой вектора в неравенстве для ошибки, приведем пример с матрицей A, вектором правой части b и вектором ошибок в правой части delta_b.

A = [2.1005632211074903e+009  2.4405728346208658e+009  2.6288320556007657e+009  3.9543845780862803e+009
     5.6384551998346877e+008  6.5511289861702812e+008  7.0564654322518432e+008  1.0614591390422670e+009
     1.9531852931707768e+009  2.2693394431935925e+009  2.4443901805698233e+009  3.6769404175323849e+009
     2.2044117415857959e+009  2.5612308917984776e+009  2.7587973499916158e+009  4.1498830928366508e+009]; 
b = [ -5.7373057243726833e-001
      -1.5400413072907607e-001
      -5.3347697688693385e-001
      -6.0209490373259589e-001];
delta_b = [-0.71685839091451e-5
            0.54786619630709e-5
            0.37746931527138e-5
	        0.20850322383081e-5];

Вычислим правую и левую часть оценки с использованием первой векторной нормы, а число обусловленности матрицы по отношению к спектральной матричной норме

RIGHT = norm(delta_b,1)/norm(b,1);
c = cond(A);
x = A\b;
x_tilda = A\(b+delta_b);
delta_x = x - x_tilda;
LEFT =  norm(delta_x,1)/norm(x,1);
format short e
disp([LEFT c*RIGHT])

Получаем для левой и правой частей неравенства следующие значения
9.9484e+004 9.9323e+004

Левая часть не должна превосходить правую (по доказанному выше), однако получилась больше из-за выбора несогласованных норм.

Разберем теперь, почему число обусловленности матрицы не может быть меньше единицы. Число обусловленности матрицы A определяется как cond(A) = || A ||M * || A-1 ||M, где || A ||M — какая-то матричная норма A. Матричная норма (т.е. правило, по которому каждой матрице сопоставляется число) не может быть произвольной, она должна удовлетворять следующим четырем аксиомам:

A1. || A || ≥ 0 для любой матрицы A и || A || = 0 тогда и только тогда, когда A — нулевая матрица.

А2. || αA || = | α | * || A || для любой матрицы A и числа α.

А3. || A + B || ≤ || A || + || B || для любых матриц A и B

А4. || AB || ≤ || A || * || B || для любых матриц A и B.

Из последней аксиомы видно, что норма определена только для квадратных матриц (хотя, в приведенных выше формулах для вычисления различных норм, в принципе, нет такого ограничения). Кроме того, из последней аксиомы следует, что любая норма единичной матрицы I не меньше единицы, действительно

|| I || = || I*I || ≤ || I ||2 ⇒ || I || ≥ 1.

Тогда, опять с привлечением четвертой аксиомы, получаем, что число обусловленности матрицы всегда больше единицы (верно для числа обусловленности матрицы по отношению к произвольной матричной норме)

1 ≤ || I || = || AA-1 || ≤ || A || * || A-1 || = cond(A).

Завершим данный раздел исследованием причины появления точного равенства в оценке ошибки в решении при возмущении правой части:

и построением соответствующих примеров в MATLAB. (Приводимое ниже рассуждение содержится, например, в книге Дж. Форсайт, К Молер. Численное решение систем линейных алгебраических уравнений. М: «Мир», 1969.)

Существенную роль в рассуждениях играет теорема о сингулярном разложении матрицы, согласно которой для любой вещественной матрицы размера на существуют две ортогональные матрицы U и V размера n на n (UTU=UUT и VTV = VVT ) такие, что произведение D = UTAV является диагональной матрицей, причем можно выбрать U и V так, что

где μ1 ≥ μ2 ≥…≥μr ≥ μr+1 =…=μn=0, ,

а r — ранг матрицы A. Числа μk называют спектральными числами матрицы A. Для невырожденной матрицы A верно:

μ1 ≥ μ2 ≥ … ≥μn ≥ 0.

Следующий важный факт заключается в том, что умножение на ортогональную матрицу не изменяет евклидовой нормы вектора, т.е. для любого вектора x с n элементами и любой ортогональной матрицы U размера n на n верно равенство

|| Ux || = || x ||.

Поскольку умножение на ортогональную матрицу не изменяет спектральной нормы, то

следовательно для числа обусловленности матрицы верно равенство

У нас есть две системы: Ax = b (с точным решением x) и (с точным решением ). Очевидно, что ошибка δx является решением системы, правая часть которой является возмущением δb, т.е. системы Aδx = δb. Пусть D = UTAV есть сингулярное разложение матрицы A, тогда UDVT = A из-за того, что U и V — ортогональные матрицы. Далее

Ax = b ⇒ UDVTx = b.

Введем обозначения

x' = VTx, b' = UTb,

тогда следующие системы эквивалентны

Ax = b ⇔ Dx' = b'

Совершенно аналогично рассмотрим систему относительно ошибки

Aδx = δb ⇒ UDVTδx = δb

Вводим обозначения

δx' = VTδx, δb' = UTδb,

при которых следующие системы оказываются эквивалентны

Aδx = δb ⇔ Dδx' = δb'

Т.е. мы получили простые эквивалентные системы с диагональной матрицей, на диагонали которой стоят спектральные числа матрицы A.

Выберем теперь специальным образом правые части этих систем, именно, пусть

где β > 0, тогда соответствующее ей решение системы Dx' = b' будет

где μ1 — максимальное сингулярное число матрицы A. Схожим образом поступим и с системой Dδx' = δ b' , именно, пусть

где γ > 0, тогда соответствующее ей решение системы Dδx' = δb' будет

Из сохранения нормы вектора при его перемножении с ортогональной матрицей следует, что

β/μ1 = || x' || = || VTx || = || x || и γ/μn = || δx' || = || VTδx || = ||δx ||.

Совершенно аналогично получаем равенства:

β = || b' || = || UTb || = || b || и γ = || δb' || = || UTδb || = || δb ||.

Тогда

и

а поскольку

то окончательно получаем

Итак, для каждой матрицы A можно построить вектор правой части системы и его возмущение такие, что ошибка в решении будет произведением числа обусловленности матрицы на ошибку в правой части. В MATLAB соответствующее построение может быть выполнено и без использования сингулярного разложения (хотя оно может быть найдено при помощи функции svd).

Сначала задаем n и получаем две ортогональные матрицы U и V, выполняя QR-разложение матриц n на n из случайных чисел:

n = 4;
[U, R] = qr(rand(n));
[V, R] = qr(rand(n));

Далее создаем диагональную матрицу D заданными числами на диагонали, возьмем, к примеру, первый элемент 1010, а остальные единицы.

D = diag([1e10 ones(1,n-1)]);

После этого конструируем матрицу A, используя диагональную матрицу D и ортогональные матрицы U и V

A = U*D*V';

Число обусловленности матрицы A будет совпадать с числом обусловленности матрицы D, которое в нашем примере равно 1010

c = cond(A)
c =
  1.0000e+010

Дальше задаем β > 0 и γ > 0, которые фигурировали в приведенных выше рассуждениях:

beta = 1;
gamma = 1e-5;

и строим вектора b' и δb'

b1 = [beta zeros(1,n-1)]';
db1 = [zeros(1,n-1) gamma]';

по которым находим вектора b и δb

b = U*b1;
db = U*db1;

Далее решаем системы Ax = b и

x = A\b;
x_tilda = A\(b+db);
вычисляем левую и правую части неравенства

dx = x - x_tilda;
RIGHT = norm(db)/norm(b);
LEFT =  norm(dx)/norm(x);

и выводим их

format short e
disp([LEFT c*RIGHT])

Получаем равенство

1.0000e+005 = 1.0000e+005

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

Система Orphus

Яндекс.Метрика