16.0 3.0 2.0 13.0

5.0 10.0 11.0 8.0

9.0 6.0 7.0 12.0

4.0 15.0 14.0 1.0

Сохраним этот файл под именем magik.dat. Тогда команда load magik.dat прочитает этот файл и создаст переменную magik, содержащую нашу матрицу.

ОБЪЕДИНЕНИЕ

Объединение – это процесс соединения маленьких матриц для создания больших. Фактически, вы создали вашу первую матрицу объединением её отдельных элементов. Пара квадратных скобок – это оператор объединения. Например, начнем с матрицы А (магического квадрата 4×4) и сформируем

В = [А А+32; А+48 А+16]

Результатом будет матрица 8×8, получаемая соединением четырех подматриц

B =

16 3 2 13 48 35 34 45

5 10 11 8 37 42 43 40

9 6 7 12 41 38 39 44

4 15 14 1 36 47 46 33

64 51 50 61 32 19 18 29

53 58 59 56 21 26 27 24

57 54 55 60 25 22 23 28

52 63 62 49 20 31 30 17

Это матрица лишь наполовину является магической. Её элементы представляют собой комбинацию целых чисел от 1 до 64, а суммы в столбцах точно равны значению для магического квадрата 8×8.

sum (В)

ans =

260 260 260 260 260 260 260 260

Однако, суммы в строках этой матрицы (sum(B’)’ ) не все одинаковы. Необходимо провести дополнительные операции, чтобы сделать эту матрицу действительно магическим квадратом 8×8.

УДАЛЕНИЕ СТРОК И СТОЛБЦОВ

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

X = А;

Теперь удалим второй столбец матрицы X .

X(:,2) =

Эта операция изменит X следующим образом

X =

16 2 13

5 11 8

9 7 12

4 14 1

Если вы удаляете один элемент матрицы, то результат уже не будет матрицей. Так выражение

X(1,2) =

результатом вычисления выдаст ошибку. Однако использование одного индекса удаляет отдельный элемент или последовательность элементов и преобразует оставшиеся элементы в вектор-строку. Так

X(2:2:10) =

выдаст результат

X =

16 9 2 7 13 12 1

ПЕРЕМНОЖЕНИЕ МАТРИЦ

При перемножении двух матриц используется оператор ‘*’. Например, если

A =

16 3 2 13

5 10 11 8

9 6 7 12

4 15 14 1

B =

16 4 7 3

5 -7 2 9

0 8 23 65

-7 4 17 9

Тогда С = А*В даст результат

С =

180 111 385 322

74 70 444 892

90 98 440 644

132 27 397 1066

Также в системе MATLAB предусмотрена возможность поэлементного перемножения. Для этой цели используется точка перед знаком умножения. Например:

C = A.*B

в результате

C =

256 12 14 39

25 -70 22 72

0 48 161 780

-28 60 238 9

СОЗДАНИЕ М-ФАЙЛОВ

M-файлы являются обычными текстовыми файлами, которые создаются с помощью текстового редактора. Для операционной среды персонального компьютера система MATLAB поддерживает специальный встроенный редактор/отладчик, хотя можно использовать и любой другой текстовый редактор с ASCII-кодами.

Открыть редактор можно двумя способами:

из меню File выбрать опцию New , а затем M-File .

использовать команду редактирования edit .

М-функции являются M-файлами, которые допускают наличие входных и выходных аргументов. Они работают с переменными в пределах собственной рабочей области, отличной от рабочей области системы MATLAB .

Пример

Функция average – это достаточно простой M-файл, который вычисляет среднее значение элементов вектора:

function y = average (x)

% AVERAGE Среднее значение элементов вектора.

% AVERAGE(X), где X – вектор. Вычисляет среднее значение элементов

% вектора.

% Если входной аргумент не является вектором, генерируется ошибка.

Size(x);

If (~((m == 1) | (n == 1)) | (m == 1 & n == 1))

Error(‘Входной массив должен быть вектором’)

End

Y =sum(x)/length(x); % Собственно вычисление

Попробуйте ввести эти команды в M-файл, именуемый average.m . Функция average допускает единственный входной и единственный выходной аргументы. Для того чтобы вызвать функцию average , надо ввести следующие операторы:

z = 1:99;

average(z)

Получаем результат

ans = 50

СТАТИСТИЧЕСКИЕ ХАРАКТЕРИСТИКИ СИГНАЛОВ

Среднее значение сигнала (его постоянная составляющая) определяется по следующей формуле:

(1.1)

Среднеквадратичное отклонение (СКО, девиация, переменная составляющая) сигнала определяется по следующей формуле:

(1.2)

Значение статистической ошибки принимаемого сигнала определяется по следующей формуле:

(1.3)

Функция нормального распределения описывается следующей формулой:

(1.4)

ЗАДАНИЕ

Проработайте основные команды, изложенные выше, в системе MATLAB.

Создайте М-функцию, которая на входе получает вектор произвольной размерности с данными и возвращает:

1. среднее значение, вычисленное в соответствии с формулой (1.1), а также полученное в результате применения функции mean;

   среднеквадратичное отклонение, вычисленное в соответствии с формулой (1.2), а также полученное в результате применения функции std.

Создайте М-функцию, которая на входе получает вектор произвольной размерности с данными и возвращает значение статистической ошибки TE в соответствии с формулой (1.3).

Самостоятельно исследуйте функцию построения гистограммы hist (вызов справки по данной функции – doc hist).

Постройте график функции нормального распределения в соответствии с формулой (1.4) с помощью функций plot и fplot.

Создайте М-функцию на основе команды randn, которая генерирует случайный шум с нормальным законом распределения с заданным средним значением и среднеквадратичным отклонением.

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

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

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

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

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

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

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

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

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

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

X = A\f x = 1 1 1

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

X = mldivide(A, f)

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

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

F - A*x ans = 0 0 0

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

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

X = 0.2707 0.3439 0.3854

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

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

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

(4-7x 1) 2 +(3-11x 1) 2 +(2-12x 1) 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 вектор содержит решение трех переопределенных систем, правая часть каждой из которых является соответствующим столбцом матрицы A.

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

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

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

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

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

A = ; f = ; 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 (n 3) операций, где n - размер матрицы (т.е. за то же время, что и решение системы), а решение системы с каждым множителем разложения за O (n 2) операций. Поэтому решение N систем c с предварительным разложением занимает O (n 3)+O (n 2)N . Решение же каждой системы без разложения матрицы потребовало бы O (n 3)N арифметических действий.

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

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

A = ; 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 = ; R = chol(A) ??? Error using ==> chol Matrix must be positive definite.

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

A = ; 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(A k )0 для k =1,2,..., n -1, то существует единственная нижняя треугольная матрица L , диагональ которой состоит из единиц, и единственная верхняя треугольная матрица U такие, что A=LU .

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

A = ; = 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 (n 2) арифметических действий. Рассмотрим пример решения системы

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

A = ; f = ;

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

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 =Q T) и матрица R размера m на n со всеми нулевыми элементами под главной диагональю, что A=QR .

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

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=Q T f :

Более подробно про матричные разложения, соответствующие функции 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 посвящен раздел Разреженные матрицы), решение находится по простой формуле x k = b k /a kk , где 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 = R T R с последующим решением систем с треугольными матрицами R T и R: R T y = 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));

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

и вторая, с множителем разложения и с занесением компонент решения на соответствующие позиции в векторе 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 - ортогональная матрица (Q T Q = I), а R - верхняя треугольная. Если A размера m на n, то Q размера m на m, а R размера m на n. Далее решение системы находится так:

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

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

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

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

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

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

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

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

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

A = ; b = ; x = linsolve(A,b) Warning: Rank deficient, rank = 1, tol = 4.4686e-015. x = 0 0 2.0000

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

A = ; b = ; = 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. Разумеется, не все комбинации значений полей допустимы, например матрица не может быть треугольной и в то же самое время симметричной и положительно определенной. Верные комбинации значений полей собраны в следующей таблице

Если матрица системы положительно определена, то это обязательно учесть при решении, поскольку для положительно определенных матриц решение основано на разложении Холецкого, которое требует меньше операций, чем 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 = ; % решаем систему как систему с симметр и полож. опред матрицей opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimePosSym = ; end % выводим графики временных затрат figure loglog(N, TimeGeneral, N, TimePosSym) legend("TimeGeneral", "TimePosSym")

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

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

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

A = ; b = ; задать 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 = ;

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

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

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

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

Сответственно

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

т.е. спектральная матричная норма равна квадратному корню из максимального собственного числа матрицы A T A. В 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

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;

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

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

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

A = ;

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

Det(A) ans = 9.9970e-121

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

X = ; b = A*x; delta_b = 1e-5*; 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 , число обусловленности матрицы системы равно 10 10 , а в решении получается относительная ошибка 10 5 . Можно проделать массу примеров, перебирая различные матрицы, не получить достижения равенства и ошибочно заключить, что это всего лишь завышенная оценка сверху для ошибки в решении. Однако, это не так, в чем нас убеждает следующий пример

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

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

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

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 равно 10 10 ;

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

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

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

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

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

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

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

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

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

Именно поэтому в разделе мы проводили все эксперименты с привлечением данных норм.

Получив неравенство || 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 = ; 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()

Получаем для левой и правой частей неравенства следующие значения
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 (U T U=UU T и V T V = VV T) такие, что произведение D = U T AV является диагональной матрицей, причем можно выбрать 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 = U T AV есть сингулярное разложение матрицы A, тогда UDV T = A из-за того, что U и V - ортогональные матрицы. Далее

Ax = b ⇒ UDV T x = b.

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

x" = V T x, b" = U T b,

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

Ax = b ⇔ Dx" = b"

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

Aδx = δb ⇒ UDV T δx = δb

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

δx" = V T δx, δb" = U T δb,

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

Aδx = δb ⇔ Dδx" = δb"

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

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

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

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

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

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

β/μ 1 = || x" || = || V T x || = || x || и γ/μ n = || δx" || = || V T δx || = ||δx ||.

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

β = || b" || = || U T b || = || b || и γ = || δb" || = || U T δb || = || δb ||.

а поскольку

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

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

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

N = 4; = qr(rand(n)); = qr(rand(n));

D = diag();

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

A = U*D*V";

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

Beta = 1; gamma = 1e-5;

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

B1 = "; db1 = ";

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

X = A\b; x_tilda = A\(b+db);

вычисляем левую и правую части неравенства

и выводим их

Format short e disp()

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

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

Работа с системой носит диалоговый характер и происходит по правилу «задал вопрос - получил ответ». Пользователь набирает на клавиатуре вычисляемое выражение, редактирует его (если нужно) в командной строке и завершает ввод нажатием клавиши ENTER. В качестве примера на рисунке показаны простейшие и вполне очевидные вычисления.

Даже из таких простых примеров можно сделать некоторые поучительные выводы:

* для указания ввода исходных данных используется символ >>;

* данные вводятся с помощью простейшего строчного редактора;

* для блокировки вывода результата вычислений некоторого выражения после него надо установить знак; (точка с запятой);

* если не указана переменная для значения результата вычислений, то MATLAB назначает такую переменную с именем ans;

* знаком присваивания является привычный математикам знак равенства =, а не комбинированный знак:=, как во многих других языках программирования и математических системах;

* встроенные функции (например, sin) записываются строчными буквами, и их аргументы указываются в круглых скобках;

* результат вычислений выводится в строках вывода (без знака >>);

* диалог происходит в стиле «задал вопрос - получил ответ».

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

В большинстве математических систем вычисление sin(V) или exp(V), где V - вектор, сопровождалось бы выдачей ошибки, поскольку функции sin и exp должны иметь аргумент в виде скалярной величины. Однако MATLAB - матричная система, а вектор является разновидностью матрицы с размером 1Чn или nЧ1. Поэтому в нашем случае результат вычислений будет вектором того же размера, что и аргумент V, но элементы возвращаемого вектора будут синусами или экспонентами от элементов вектора V.

Матрица задается в виде ряда векторов, представляющих ее строки и заключенных в квадратные скобки. Для разделения элементов векторов используется пробел или запятая, а для отделения одного вектора от другого - точка с запятой. Для выделения отдельного элемента матрицы M используется выражение вида M(j,i), где M - имя матрицы, j - номер строки и i - номер столбца.

Для просмотра содержимого массивов удобно использовать браузер рабочего пространства Workspace. Каждый вектор и матрица в нем представляются в виде квадратика с ячейками, справа от которого указывается размер массива. Двойной щелчок по квадратику мышью ведет к появлению окна редактора массивов Array Editor. Работа с редактором массивов вполне очевидна - возможен не только просмотр элементов массивов, но и их редактирование и замена.

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

Работа с редактором массивов

To get started, select "MATLAB Help" from the Help menu.

>> type sin

sin is a built-in function.

>> help sin

SIN(X) is the sine of the elements of X.

Overloaded methods

>> V=

0.8415 0.9093 0.1411 -0.7568

Error using ==> ^

Matrix must be square.

Можно обратить внимание на форму ответов при выполнении простых операций без указания переменной, которой присваивается результат. В таких случаях MATLAB сам назначает переменную ans, которой присваивается результат и значение которой затем выводится на экран.

Форма вывода и перенос строки в сессии

Следует отметить особенности вывода в системе MATLAB. Вывод начинается с новой строки, причем числовые данные выводятся с отступом, а текстовые - без него. Для экономии места в данной книге в дальнейшем вывод будет даваться без перевода на новую строку. Например, вывод вектора-строки

будет дан в виде:

Исключением является вывод векторов столбцов и матриц - тут будет сохранена более наглядная и присущая MATLAB по умолчанию форма вывода.

В некоторых случаях вводимое математическое выражение может оказаться настолько длинным, что для него не хватит одной строки. Тогда часть выражения можно перенести на новую строку с помощью знака многоточия «...» (3 или более точек), например:

s = 1 - 1/2 + 1/3 - 1/4 + 1/5 - 1/6 + 1/7 ...

1/8 + 1/9 - 1/10 + 1/11 - 1/12;

Максимальное число символов в одной строке командного режима - 4096, а в m-файле - не ограничено, но со столь длинными строками работать неудобно. В ранних версиях в одной строке было не более 256 символов.

Запуск примеров применения MATLAB из командной строки

MATLAB имеет множество примеров применения, часть из которых можно запускать прямо из командной строки. Например, команда

запускает m-файл bench.m демонстрационного примера тестирования системы.

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

Получение краткой справки (help ) в командном окне MATLAB производится с помощью команд:

>> help - вывод сведений о разделах (topics) MATLAB с возможностью

гипертекстового перехода к выводу списков функций каждого раздела и справочных сведений по необходимой функции.

>> help - вывод сведений об именах (названиях) функций, входящих в раздел.

>> help <имя функции> - вывод справочных сведений по функции.

>> helpwin - вывод окна справки, в котором двойным щелчком мыши можно открыть справочные сведения по нужному разделу или функции.

>> lookfor - вывод справочных сведений по ключевому слову.

>> help demos - вывод списка демонстрационных примеров.

>> hthelp- открывает интерактивное окноMATLAB help .

>> help symbolic– выводит сведения об инструментарии символьной математики (symbolicmathtoolbox)MATLAB.

>> help signal processing toolbox – выводит сведения о функциях пакета Signal Processing Toolbox.Чаще других для справок используется команда help <имя функции>.

Пример .

>> help abs

ABS Absolute value.

ABS(X) is the absolute value of the elements of X. When

X is complex, ABS(X) is the complex modulus (magnitude) of

the elements of X.

See also SIGN, ANGLE, UNWRAP.

Overloaded methods

help iddata/abs.m

Большинство функций имеют несколько вариантов синтаксиса. Имена функций в сообщениях справочной системы выводятся прописными символами, но при вводе имен функций должны использоваться только строчные символы . Вывод справки по необходимой функции сопровождается списком связанных (сопутствующих) функций. Для получения более подробных сведений по нужной функции с примерами вычислений служит командаdoc <имя функции>.

Главным средством для получения развернутой справочной информации служит браузер справочной системы Help browser , который содержит документацию по всем инсталлированным продуктамMATLAB. Доступ к документации осуществляется через меню HELP. На начальном этапе работы для знакомства с пакетом особенно полезен и необходим раздел MATLAB справочной системы.

Из меню Help c помощью команды Demos можно получить доступ к демонстрационным примерам MATLAB. Эти примеры очень разнообразны и полезны для целей обучения и создания приложений в MATLAB. Доступ к демонстрационным примерам можно получить также с помощью команды demo из командной строки.

Доступ к справочным сведениям в Internet: >>webhttp://www . mathworks . com - загружаетWEB- сайт фирмыMathWorksInc. - производителяMATLAB.

1. Простые вычисления

MATLAB имеет следующие базовые арифметические операции:

Сложение (a+b , 15+23),

Вычитание (a-b , 17-3),

Умножение (a*b , 0.18*6.12),

Деление (a/b , 92.4/15),

Возведение в степень (a ^ b , 7.4^4).

Примеры

Name Size Bytes Class Attributes

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

Для очистки рабочего пространства Workspace, т.е. удаления из него переменных, можно применить команду clear. Команда clc используется для очистки рабочего окна Command Window без очистки рабочей области.

Пакет MATLAB поддерживает также математические функции общего назначения, такие как извлечение квадратного корня sqrt (x ), вычисление прямых и обратных тригонометрических функций, экспоненциальных функций и др. Перечень всех этих функций с возможностью перехода к любой из них можно получить, введя в командную строку help elfun . Все элементарные функцииMATLABявляются функциями, аргументами которых могут быть массивы, т.е. в пакете реализованавекторизация вычислений.

Пример

>> v1 = [ 2 4 sqrt(10)]

2.0000 4.0000 3.1623

0.4161 -0.6536 -0.9998

MATLAB вычисляет выражения слева направо в обычном порядке приоритета операции возведения в степень над операциями умножения и деления и последних – над операциями сложения и вычитания. Для изменения порядка вычислений используются скобки.

Пример

>> 7*3+5-12/4

>> 7*(3+5-12/4)

Относительная точность выполнения арифметических операций MATLABсоставляет около 16 десятичных цифр в диапазоне чисел от 10 -308 до 10 308 . По умолчанию вMATLABустанавливается формат выводаshort , позволяющий вывести не более 5 значащих цифр числа. Такой формат вывода не всегда достаточен.

Команды установки форматов вывода

>> format short- короткое представление в фиксированном формате (5 знаков),

>> format short e– устанавливает формат научной (экспоненциальной) нотации с 5 десятичными разрядами,

>> format long– формат длинного представления с фиксированной точкой с 15 десятичными разрядами,

>> format long e– формат научной нотации с 15 десятичными разрядами,

>> format bank- денежный формат вывода с двумя десятичными разрядами справа от десятичной точки,

>> format rat- формат вывода в виде рациональной дроби.

Формат вывода может быть также установлен командой меню Preferences .

Следует иметь ввиду, что при вводе чисел в экспоненциальной форме, например, 15.8e-5промежуточные пробелы не допускаются.

Имена переменных MATLABдолжны начинаться с буквы, максимальная длина имени - 31 символ. Имена не должны совпадать с именами функций и процедур и системных переменных. Имена чувствительны к регистру, например,var отличается отVar .

Для создания переменной используется операция присваивания

<имя переменной> = <выражение>;

При этом оператор “;” подавляет эхо-вывод результатов вычисления (присваивания) на экран.

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

Примеры

Символьная переменная >> string="hello"

Действительные скаляры (числа)

>> y=5.2*x+15

Для сохранения переменных в файле текущего каталога (по умолчанию папка work) можно использовать команду save

>> save myfile x y

Без указания имен переменных команда save сохраняет все переменные рабочего пространства.

Переменные можно удалять из рабочего пространства (Workspace )MATLABкомандойclear

>> clear x y

Undefined function or variable "x".

Undefined function or variable "y".

При необходимости можно загрузить переменные из файла в рабочее пространство командой load

>> load myfile

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

>> sqrt(-3)

Для обозначения мнимой единицы в MATLABзарезервированы переменныеi иj :

3.0000 + 4.0000i

>>y= 2*(1+4*j)

2.0000 + 8.0000i

Специальные функции вычислений с комплексным аргументом:

>> abs(x)% получение модуля числа

>> angle(x)% аргумент (фаза) числа в радианах

>> conj(x)% комплексно сопряженное число

>> imag(x)% мнимая часть числа

>> real(x)% действительная часть числа

Имена предопределенных системных переменных MATLAB нельзя использовать в качестве имен пользовательских переменных. Основные из этих имен:

>> ans имя переменной по умолчанию для результата вычислений.

>> eps переменная машинной точности вычислений, имеет порядок 10 -16 .

>> exit завершение (окончание) работы MATLAB.

>> i или j мнимая единица, т.е. .

>> pi число π.

>> Inf обозначение бесконечности.

>> NaN не числовой результат.

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

>> clear x,y команда удаляет переменные x и y.

>> what вывод списка файлов с расширениями ‘.m’, ‘.mat’, ‘.mex’ из текущего каталога.

>> who отображает переменные текущего рабочего пространства.

>> whos выводит информацию о текущих переменных.

>> dir выводит список файлов текущего каталога.

>> save сохраняет все текущие переменные в файле MATLAB.mat в текущем каталоге.

>> load загружает переменные из MATLAB.mat в текущий сеанс работы.

>> diary сохраняет текст (команды) и результаты вычислений текущего сеанса работы (дневник сессии) в файле с именем diary.

>> diary filename сохраняет текущий сеанс в файле с именем filename.

>> diary off приостанавливает запись в файл.

>> diary on включает запись сессии в файл.