Система математических расчетов MATLAB

plot3(X(1,:), X(2,:), X(3,:), '-o')

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

Собственные значения и собственные векторы

Собственным значением и собственным вектором квадратной матрицы А называются ска-ляр л и вектор v, удовлетворяющие условию

Av = лv

Диагональная декомпозиция

Имея диагональную матрицу Л, составленную из собственных значений л матрицы А и мат-рицу V , составленную из соответствующих собственных векторов v, можно записать

AV = VЛ

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

А = VЛV-1

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

lambda = eig(A)

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

lambda =

-3.0710

-2.4645 + 17.6008i

-2.4645 - 17.6008i

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

При двух выходных аргументах, функция eig вычисляет также собственные векторы и выда-ет собственные значения в виде диагональной матрицы

.

[V,D] = eig(A)

V =

-0.8326 0.2003 - 0.1394i 0.2003 + 0.1394i

-0.3553 -0.2110 - 0.6447i -0.2110 + 0.6447i

-0.4248 -0.6930 -0.6930

D =

-3.0710 0 0

0 -2.4645+17.6008i 0

0 0 -2.4645-17.6008i

Первый собственный вектор (первый столбец матрицы V) является действительным, а два других являются комплексно-сопряженными. Все три вектора являются нормализованными по длине, т.е. их Евклидова норма norm(v,2), равна единице.

Матрица V*D*inv(V), которая в более сжатой форме может быть записана как V*D/V, равна, в пределах погрешностей округления, матрице А. Аналогично, inv(V)*A*V, или V\A*V, рав-на, в пределах погрешностей округления, матрице D.

Дефектные матрицы

Некоторые матрицы не имеют спектрального разложения. Такие матрицы называются дефек-тными или не диагонализируемыми. Например, пусть матрица А имеет вид

A =

6 12 19

-9 -20 -33

4 9 15

Для этой матрицы ввод [V, D] = eig(A) дает

V =

-0.4741 -0.4082 -0.4082

0.8127 0.8165 0.8165

-0.3386 -0.4082 -0.4082

D =

-1.0000 0 0

0 1.0000 0

0 0 1.0000

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

Сингулярное разложение матриц

Сингулярным значением и соответствующими сингулярными векторами прямоугольной ма-трицы A называются скаляр у и пара векторов u и v такие, что удовлетворяются соотноше-ния

Av = уu

ATu = уv

Имея диагональную матрицу сингулярных чисел У и две ортогональные матрицы U и V, сформированные из соответствующих собственных векторов, можно записать

AV = U У

ATU = V У

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

A = U УVT

Полное сингулярное разложение матрицы А размера mхn включает mхm матрицу U, mхn матрицу У, и nхn матрицу V. Другими словами, обе матрицы U и V являются квадратными , а матрица У имеет тот же размер, что и A. Если A имеет намного больше строк чем столб-цов, результирующая матрица U может быть достаточно большой, но большинство ее столб-цов умножаются на нули в У . В таких ситуациях может быть использована так называемая экономичная декомпозиция, которая сберегает как время так и память, за счет вывода матри-цы U размера mхn, матрицы У размера nхn и той же матрицы V.

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

Для матрицы

A =

9 4

6 8

2 7

полное сингулярное разложение задается в форме

[U,S,V] = svd(A)

и приводит к следующим результатам

U =

-0.6105 0.7174 0.3355

-0.6646 -0.2336 -0.7098

-0.4308 -0.6563 0.6194

S =

14.9359 0

0 5.1883

0 0

V =

-0.6925 0.7214

-0.7214 -0.6925

Вы можете убедиться, что матрица U*S*V' равна А с точностью до ошибок округления. Для этого примера экономичная декомпозиция дает незначительный эффект.

[U,S,V] = svd(A,0)

U =

-0.6105 0.7174

-0.6646 -0.2336

-0.4308 -0.6563

S =

14.9359 0

0 5.1883

V =

-0.6925 0.7214

-0.7214 -0.6925

Как и в первом случае, матрица U*S*V' равна A с точностью до ошибок округления.

Полиномы и интерполяция

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

Полиномы и действия над ними

Представление полиномов

MATLAB представляет полиномы как векторы-строки, содержащие коэффициенты полино-мов по убывающим степеням. Например, рассмотрим следующее уравнение

p(x) = x3 - 2x - 5

Это известный пример Валлиса (Wallis), использованный при первом представлении метода Ньютона во Французкой Академии. Мы будем использовать его в дальнейшем при рассмот-рении примеров использования различных функций. Для ввода данного полинома в MATLAB, следует записать

p = [1 0 -2 -5].

Корни полинома

Корни полинома вычисляются при помощи функци roots :

r = roots(p)

r =

2.0946

-1.0473 + 1.1359i

-1.0473 - 1.1359i

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

p2 = poly(r)

p2 =

1 8.8818e-16 -2 -5

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

Характеристические полиномы

Функция poly вычисляет также коэффициенты характеристического полинома матрицы:

A = [1.2 3 -0.9; 5 1.75 6; 9 0 1];

poly(A)

ans =

1.0000 -3.9500 -1.8500 -163.2750

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

Вычисление значений полинома

Функция polyval вычисляет значение полинома в заданных точках. Для вычисления p в точ-ке s = 5, следует записать

polyval(p,5)

ans =

110

Можно также вычислить значение матричного полинома. Так, вместо полинома Валлиса мо-жно записать:

p(X) = X3 - 2X - 5I

где X является квадратной матрицей, а I - единичной матрицей. Например, сформируем сле-дующую квадратную матрицу X

X = [2 4 5; -1 0 3; 7 1 5];

и вычислим значение заданного выше полинома p(X) на данной матрице.

Y = polyvalm(p, X)

Y =

377 179 439

111 81 136

490 253 639

Умножение и деление полиномов

Для умножения и деления полиномов предназначены соответственно функции conv и deconv. Рассмотрим полиномы a(s) = s2 + 2s + 3 и b(s) = 4s2 + 5s + 6. Для вычисления их произведения следует ввести

a = [1 2 3]; b = [4 5 6];

c = conv(a,b)

MATLAB возвращает

c =

4 13 28 27 18

Для получения из с полинома b воспользуемся функцией deconv:

[q, r] = deconv(c, a)

q =

4 5 6

r =

0 0 0 0 0

где r - остаток после деления (в данном случае нулевой вектор). В общем случае для поли-номов q, r , c, a в функции deconv справедливо соотношение

c = conv(q, a) + r

Вычисление производных от полиномов

Функция polyder вычисляет производную любого полинома. Для получения производной от нашего полинома p = [1 0 -2 -5], введем

q = polyder(p)

q =

3 0 - 2

Функция polyder вычисляет также производные от произведения или частного двух полино-мов. Например, создадим два полинома a и b:

a = [1 3 5]; b = [2 4 6];

Вычислим производную произведения a*b вводом функции polyder с двумя входными аргу-ментами a и b и одним выходным:

c = polyder(a, b)

c =

8 30 56 38

Вычислим производную от частного a/b путем ввода функции polyder с двумя выходными аргументами:

[q, d] = polyder(a, b)

q =

-2 -8 -2

d =

4 16 40 48 36

где отношение двух полиномов q/d является результатом операции дифференцирования.

Аппроксимация кривых полиномами

Функция polyfit находит коэффициенты полинома заданной степени n , который аппрокси-мирует данные (или функцию y(x)) в смысле метода наименьших квадратов:

p = polyfit(x, y, n)

где x и y есть векторы, содержащие данные x и y, которые нужно аппроксимировать полино-мом. Например, рассмотрим совокупность данных x-y, полученную экспериментальным пу-тем

x = [1 2 3 4 5]; y = [5.5 43.1 128 290.7 498.4].

Аппроксимация функциональной зависимости y(x) в виде полинома третьего порядка

p = polyfit(x,y,3)

дает коэффициенты полинома

p =

-0.1917 31.5821 -60.3262 35.3400

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

x2 = 1 : 0.1 : 5;

y2 = polyval(p, x2);

plot(x, y, 'o', x2, y2); grid on

где функция grid on служит для нанесения координатной сетки, а экспериментальные дан-ные на графике отмечены маркерами о.

.

Как видно из рисунка, полином третьего порядка достаточно хорошо аппроксимирует наши данные.

Разложение на простые дроби

Функция residue вычисляет вычеты, полюса и многочлен целой части отношения двух поли-номов. Это особенно полезно при представлении систем управления в виде передаточных функций. Для полиномов a(s) и b(s), при отсутствии кратных корней имеем

где r есть вектор-столбец вычетов, p есть вектор-столбец полюсов, а k есть вектор-строка це-лой части дробно-рациональной функции. Рассмотрим передаточную функцию

Для полиномов числителя и знаменателя этой функции имеем:

b = [-4 8]; a = [1 6 8].

Введя

[r, p, k] = residue(b, a)

получим

r =

-12

8

p =

-4

-2

k =

[ ]

Функция residue с тремя входными (r, p, и k) и двумя выходными (b2, a2) аргументами вы-полняет обратную функцию свертки имеющегося разложения на простые дроби, в дробно-рациональную функцию отношения двух полиномов.

[b2, a2] = residue(r, p, k)

b2 =

-4 8

a2 =

1 6 8

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

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

Интерполяция

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

Одномерная интерполяция

Двумя основными типами одномерной интерполяции в MATLAB-е являются полиномиаль-ная интерполяция и интерполяция на основе быстрого преобразования Фурье.

1. Полиномиальная интерполяция

Функция interp1 осуществляет одномерную интерполяцию - важную операцию в области анализа данных и аппроксимации кривых. Эта функция использует полиномиальные методы, аппроксимируя имеющийся массив данных полиномиальными функциями и вычисляя соот-ветствующие функции на заданных (желаемых) точках. В наиболее общей форме эта функ-ция имеет вид

yi = interp1(x, y, xi, method)

где y есть вектор, содержащий значения функции; x - вектор такой же длины, содержащий те точки (значения аргумента), в которых заданы значения y; вектор xi содержит те точки, в ко-торых мы хотим найти значения вектора y путем интерполяции; method - дополнительная строка, задающая метод интерполяции. Имеются следующие возможности для выбора мето-да:

* Ступенчатая интерполяция (method = 'nearest'). Этот метод приравнивает значение функ-ции в интерполируемой точке к ее значению в ближайшей существующей точке имеющихся данных.

* Линейная интерполяция (method = 'linear'). Этот метод аппроксимирует функцию между любыми двумя существующими соседними значениями как линейную функцию, и возвр-ащает соответствующее значение для точки в xi (метод используется по умолчанию).

* Интерполяция кубическими сплайнами (method = 'spline'). Этот метод аппроксимирует ин-терполируемую функцию между любыми двумя соседними значениями при помощи куби-ческих функций, и использует сплайны для осуществления интерполяции.

* Кубическая интерполяция (method = 'pchip' или 'cubic'). Эти методы идентичны. Они ис-пользуют кусочную кубическую Эрмитову аппроксимацию и сохраняют монотонность и форму данных.

Если какой-либо из элементов вектора xi находится вне интервала, заданного вектором x, то выбранный метод интерполяции используется также и для экстраполяции. Как альтернатива,

функция yi = interp1(x, y, xi, method, extrapval) заменяет экстраполированные значения теми, которые заданы вектором extrapval. Для последнего часто используется нечисловое значение NaN.

Все методы работают на неравномерной сетке значений вектора x .

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

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

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

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

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

Относительные качественные характеристики всех перечисленных методов сохраняются и в случае двух- или многомерной интерполяции.

2. Интерполяция на основе быстрого преобразования Фурье _

Функция interpft осуществляет одномерную интерполяцию с использованием быстрого пре-образование Фурье (FFT). Этот метод вычисляет преобразование Фурье от вектора, который содержит значения периодической функции. Затем вычисляется обратное преобразование Фурье с использованием большего числа точек. Функция записывается в форме

y = interpft(x, n)

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

Двумерная интерполяция

Функция interp2 осуществляет двумерную интерполяцию - важную операцию при обработке изображений и графического представления данных. В наиболее общей форме эта команда имеет вид

ZI = interp2(X, Y, Z, XI, YI, method)

где Z есть прямоугольный массив, содержащий значения двумерной функции; X и Y являют-ся массивами одинаковых размеров, содержащие точки в которых заданы значения двумер-ной функции; XI и YI есть матрицы, содержащие точки интерполяции (то есть промежуточ-ные точки, в которых нужно вычислить значения функции); method - строка, определяющая метод интерполяции. В случае двумерной интерполяции возможны три различных метода:

* Ступенчатая интерполяция (method = 'nearest'). Этот метод дает кусочно-постоянную поверхность на области значений. Значение функции в интерполируемой точке равно значе-нию функции в ближайшей заданной точке.

* Билинейная интерполяция (method = 'linear'). Метод обеспечивает аппроксимацию данных при помощи билинейной поверхности (плоскости) на множестве заданных значений двумер-ной функции. Значение в точке интерполяции является комбинацией значений четырех бли-жайших точек. Данный метод можно считать «кусочно-билинейным»; он быстрее и требует меньше памяти, чем бикубическая интерполяция.

* Бикубическая интерполяция (method = 'cubic'). Данный метод аппроксимирует поверх-ность при помощи бикубических поверхностей. Значение в точке интерполяции является комбинацией значений в шестнадцати ближайших точках. Метод обеспечивает значительно более гладкую поверхность по сравнению с билинейной интерполяцией. Это может быть ключевым преимуществом в приложениях типа обработки изображений. Особенно эффек-тивным данный метод является в ситуациях, когда требуется непрерывность как интерполи-руемых данных, так и их производных.

Все эти методы требуют, чтобы X и Y были монотонными, то есть или всегда возрастающи-ми или всегда убывающими от точки к точке. Эти матрицы следует сформировать с исполь-зованием функции meshgrid, или же, в противном случае, нужно убедиться, что «схема» то-чек имитирует сетку, полученную функцией meshgrid. Перед интерполяцией, каждый из указанных методов автоматически отображает входные данные в равномерно распреде-ленную сетку. Если X и Y уже распределены равномерно, вы можете ускорить вычисления добавляя звездочку к строке метода, например, '*cubic'.

Сравнение методов интерполяции

Приведенный ниже пример сравнивает методы двумерной интерполяции в случае матрицы данных размера 7х7.

1. Сформируем функцию peaks на «грубой» сетке (с единичным шагом).

[x, y] = meshgrid(-3 : 1 : 3);

z = peaks(x,y);

surf(x,y,z)

где функция meshgrid(-3:1:3) задает сетку на плоскости x и y в виде двумерных массивов размера 7х7; функция peaks(x,y) является двумерной функцией, используемой в MATLAB-е в качестве стандартных примеров, а surf(x,y,z) строит окрашенную параметрическую повер-хность. Соответствующий график показан ниже.

2. Создадим теперь более мелкую сетку для интерполяции (с шагом 0.25).

[xi,yi] = meshgrid(-3:0.25:3);

3. Осуществим интерполяция перечисленными выше методами.

zi1 = interp2(x,y,z,xi,yi,'nearest');

zi2 = interp2(x,y,z,xi,yi,'bilinear');

zi3 = interp2(x,y,z,xi,yi,'bicubic');

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

surf(xi,yi,zi1) surf(xi,yi,zi2) surf(xi,yi,zi3)

Метод `nearest' Метод `bilinear' Метод `bicubic'

Интересно также сравнить линии уровней данных поверхностей, построенных при помощи специальной функции contour.

contour(xi,yi,zi1) contour(xi,yi,zi2) contour(xi,yi,zi3)

Метод `nearest' Метод `bilinear' Метод `bicubic'

Отметим, что бикубический метод производит обычно более гладкие контуры. Это, однако, не всегда является основной заботой. Для некоторых приложений, таких, например, как об-работка изображений в медицине, метод типа ступенчатой интерполяции может быть более предпочтительным, так как он не «производит» никаких «новых» результатов наблюдений.

Анализ данных и статистика

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

· Optimization - Нелинейные методы обработки данных и оптимизация.

· Signal Processing - Обработка сигналов, фильтрация и частотный анализ.

· Spline - Аппроксимация сплайнами.

· Statistics - Углубленный статистический анализ, нелинейная аппроксимация и

регрессия.

· Wavelet - Импульсная декомпозиция сигналов и изображений.

Внимание ! MATLAB выполняет обработку данных, записанных в виде двумерных массивов по столбцам ! Одномерные статистические данные обычно хранятся в отдельных векорах, причем n-мерные векторы могут иметь размерность 1х n или nх1. Для многомерных данных матрица является естественным представлением, но здесь имеются две возможности для ориентации данных. По принятому в системе MATLAB соглашению, различные пере-менные должны образовывать столбцы, а соответствующие наблюдения - строки. Поэтому, например, набор данных, состоящий из 24 выборок 3 переменных записывается в виде мат-рицы размера 24х3.

Основные функции обработки данных

Перечень функций обработки данных, расположенных в директории MATLAB-а datafun приведен в Приложении 8 .

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

count =

11 11 9

7 13 11

14 17 20

11 13 9

43 51 69

38 46 76

61 132 186

75 135 180

38 88 115

28 36 55

12 12 14

18 27 30

18 19 29

17 15 18

19 36 48

32 47 10

42 65 92

57 66 151

44 55 90

114 145 257

35 58 68

11 12 15

13 9 15

10 9 7

Таким образом, мы имеем 24 наблюдения трех переменных. Создадим вектор времени, t, со-стоящий из целых чисел от 1 до 24: t = 1 : 24. Построим теперь зависимости столбцов матри-цы counts от времени и надпишем график:

plot(t, count)

legend('Location 1','Location 2','Location 3',0)

xlabel('Time')

ylabel('Vehicle Count')

grid on

где функция plot(t, count) строит зависимости трех векторов-столбцов от времени; функция

legend('Location 1','Location 2','Location 3',0) показывает тип кривых; функции xlabel и ylabel надписывают координатные оси, а grid on выводит координатную сетку. Соответству-ющий график показан ниже.

Применим к матрице count функции max (максимальное значение), mean (среднее значение) и std (стандартное, или среднеквадратическое отклонение).

mx = max(count)

mu = mean(count)

sigma = std(count)

В результате получим

mx =

114 145 257

mu =

32.00 46.5417 65.5833

sigma =

Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9