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

Spline Toolbox

Обзор средств MATLAB и ToolBox'ов для приближения данных

Приближение полиномами в смысле наименьших квадратов

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

Задача о приближении данных полиномом в смысле наименьших квадратов

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

требуется найти такой полином степени

коэффициенты которого являются решением следующей задачи минимизации

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

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

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

Поставим задачу приблизить данные, которые заданы массивами x и y:

>> x = [0.1 0.3 0.45 0.5 0.79 1.1 1.89 2.4 2.45];
>> y = [-3 -1 0.9 2.4 2.5 1.9 0.1 -1.3 -2.6];

полиномами первой, третьей и пятой степени. Используем функцию polyfit, указав в ее входных аргументах вектора с данными и степени полиномов. Сами коэффициенты запишем в вектора p1, p3 и p5, соответственно:

>> p1 = polyfit(x, y, 1)
p1 =
   -0.6191    0.6755
>> p3 = polyfit(x, y, 3)
p3 =
    2.2872  -12.1553   17.0969   -4.5273
>> p5 = polyfit(x, y, 5)
p5 =
   -6.0193   33.9475  -62.4220   35.9698    4.7121   -3.8631

Итак, мы получили полиномы:

Для построения графиков этих полиномов следует найти их значения в промежуточных точках, принадлежащих интервалу, на котором заданы данные, т.е. между x(1) и x(end). Сгенерируем 100 точек, равномерно расположенных на области определения данных, при помощи функции linspace

>> xx = linspace(x(1), x(end), 100);

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

>> yy1 = polyval(p1, xx);
>> yy3 = polyval(p3, xx);
>> yy5 = polyval(p5, xx);

Для наглядности построим графики полиномов и разместим заданные массивами x и y точки круглыми маркерами

>> plot(x, y, 'o', xx, yy1, xx, yy3, xx, yy5)

Кроме этого, при помощи функции legend разместим справа от осей легенду, в которой укажем соответствие линий полиномам:

>> legend('DATA', '{\itp}^{(1)}({\itx})', '{\itp}^{(3)}({\itx})', '{\itp}^{(5)}({\itx})',-1)

При создании легенды использовались возможности интерпретатора TeX (см. Вывод математических формул). В результате получается график, приведенный ниже

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

>> [p3, S3] = polyfit(x, y, 3)

приводит к выводу в командное окно того же самого вектора коэффициентов кубического полинома, что получился в предыдущем примере, и структуры S3 с полями R, df и normr:

p3 =
    2.2872  -12.1553   17.0969   -4.5273
S3 = 
        R: [4x4 double]
       df: 5
    normr: 1.7201

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

(если работа со структурами незнакома, то достаточно понять, что для записи содержимого поля структуры в некоторую переменную надо отделять имя структуры от имени поля точкой, т.е. >> r = S3.normr). В поле df содержится разность между числом точек и числом коэффициентов полинома.

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

>> [p8, S8] = polyfit(x, y, 8)
S8 = 
        R: [9x9 double]
       df: 0
    normr: 2.2382e-010

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

>> [p7, S7] = polyfit(x, y, 7)
>> yy7 = polyval(p7, xx);
>> yy8 = polyval(p8, xx);
>> plot(x, y, 'o', xx, yy7, xx, yy8)
>> legend('DATA', '{\itp}^{(7)}({\itx})', '{\itp}^{(8)}({\itx})',-1)

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

Дальнейшее увеличение степени полинома ни к чему хорошему не приведет, например при построении полинома девятой степени

>> [p9, S9] = polyfit(x, y , 9)

выведется предупреждение о том, что такой полином не единственный

Warning: Polynomial is not unique; degree >= number of data points.

Приближение полиномом девятой степени окажется еще хуже, чем интерполяция полиномом восьмой степени:

>> yy9 = polyval(p9, xx);
>> plot(x, y, 'o', xx, yy8, xx, yy9)
>> legend('DATA', '{\itp}^{(8)}({\itx})', '{\itp}^{(9)}({\itx})',-1)

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

Метод наименьших квадратов

Напомним, что задача приближения данных

полиномом некоторой степени

состоит в решении задачи минимизации

Запишем полином в виде

Тогда задача о построении полиномиального приближения сводится к поиску минимума следующей функции

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

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

   

или в матричном виде

где

а элементы матрицы и вектора правой части системы выражаются формулами

      

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

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

Аналогично, вектор правой части системы уравнений так же выражается через вектор

следующим образом:

В функции polyfit при помощи функции qr пакета MATLAB разыскивается QR-разложение матрицы , т.е. находятся такая ортогональная матрица (квадратная матрица называется ортогональной, если произведение является единичной матрицей) и верхняя треугольная матрица , что (более точно, ищется только первые столбцов матрицы и не вся матрица , а только первые ее строк и столбцов). Далее система линейных уравнений

решается так:

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

Ошибка приближения, т.е. значение

для найденного полинома

вычисляется по формуле

где евклидова норма вектора находится при помощи функции norm.

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

  • R - содержит множитель QR-разложения матрицы ;
  • df - содержит значение
  • normr - содержит норму ошибки полиномиального приближения .

Например, для записи евклидовой нормы ошибки приближения в переменную e можно использовать следующие команды:

[p, S] = polyfit(x,y,4);
e = S.normr

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

Когда функция polyfit выдает предупреждения, плохая обусловленность, центрирование и масштабирование данных

Функция polyfit выдает предупреждения в нескольких случаях.

  1. Если степень полинома , приближающего данные в смысле наименьших квадратов больше или равна числу данных, то выводится сообщение о том, что полином не единственный.
  2. Если данные такие, что при нахождении коэффициентов полинома при решении системы линейных алгебраических уравнений возможна большая вычислительная ошибка, то выводится предупреждение. О том, что возможна ошибка, свидетельствует число обусловленности матрицы (см. разд. Метод наименьших квадратов), которое оценивается в функции polyfit при помощи функции condest (она делает оценку числа обусловленности по отношению к первой матричной норме). Большое число обусловленности (в функции polyfit порог равен 1010) может привести к большим ошибкам при решении системы и нахождении коэффициентов полинома.

Большое число обусловленности может быть в нескольких случаях. Во-первых матрица

при , как уже было замечено, является матрицей Вандермнода (с точностью до перестановки столбцов). Если точки , в которых заданы значения данных , равномерно распределены на отрезке [0, 1], то число обусловленности матрицы с ростом N ведет себя как число обусловленности матрицы Гильберта, которая является классическим примером плохо обусловленной матрицы с экспоненциально растущим числом обусловленности. Тогда число обусловленности матрицы , совпадающее по значению с числом обусловленности множителя , входящего в QR-разложение , будет вести себя как квадратный корень из числа обусловленности матрицы . Для N равномерно распределенных на отрезке [0, 1] точек зависимость числа обусловленности матрицы от N (при ) в полулогарифмической шкале представлена на следующем графике

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

	>> N = 13;
	>> x = linspace(0,1,N);
	>> y = sin(x);
	>> p = polyfit(x, y, N-1)
	Warning: Polynomial is badly conditioned. Remove repeated data points
	         or try centering and scaling as described in HELP POLYFIT.

Для получения предыдущего графика число обусловленности прямоугольной матрицы вычислялось как отношение ее наибольшего сингулярного числа к наименьшему, т.е. , где и - максимальные и минимальные числа на диагонали матрицы , входящей в сингулярное разложение матрицы , где и - ортогональные матрицы. Сами сингулярные числа определялись при помощи функции svd. Поскольку в функции polyfit используется сверху оценка для числа обусловленности, возвращаемая функцией condest, то предупреждение о превышении порога 1010 оценкой для числа обусловленности может наступить раньше. Для сравнения числа обусловленности матрицы , полученного при помощи svd, и его оценки сверху, возвращаемой функцией condest, приведен следующий график

Если степень полинома не в точности равна , но близка к числу точек N, в которых заданы данные, то число обусловленности тоже будет велико, как видно из следующего графика

К примеру, приближение данных, заданных в двадцати равноотстоящих точках на отрезке [0, 1] полиномом тринадцатой степени уже приведет к предупреждению о плохой обусловленности:

	>> N = 20;
	>> x = linspace(0, 1, N);
	>> y = sin(x);
	>> p = polyfit(x, y, 13)
	Warning: Polynomial is badly conditioned. Remove repeated data points
	         or try centering and scaling as described in HELP POLYFIT.

Следующий пример демонстрирует, к чему может привести большая обусловленность, если исходные данные заданы с ошибкой 10-3. Данные задаются на отрезке [0,1] в двадцати равномерно расположенных точках , , . Заданные значения в них являются значениями функции . Сначала производится приближение этих данных полиномом восемнадцатой степени по методу наименьших квадратов, затем в данные вносится ошибка порядка 10-3 и снова строится приближение полиномом восемнадцатой степени по методу наименьших квадратов. При этом функция polyfit каждый раз выдает предупреждение о плохой обусловленности.

% определяем данные в 20-и точках
N = 20;
x = linspace(0,1,N);
y = sin(x);
% приближаем данные по методу наименьших квадратов полиномом 18-ой степени
n=18;
p = polyfit(x, y, n)
% строим график данных и полинома
xx = linspace(x(1), x(end), 100);
yy = polyval(p, xx);
plot(x, y, 'o', xx, yy)
	% вносим в данные ошибку порядка 0.001
ye = y + 0.001*rand(size(y));
% приближаем данные по методу наименьших квадратов полиномом 18-ой степени
pe = polyfit(x, ye, n)
% строим график полинома
yy = polyval(pe, xx);
hold on
plot(xx, yy, 'm')
legend('data', '{\itp}^{(18)} for data', '{\itp}^{(18)} for data + 10^{-3}',-1)

Более того, разница между коэффициентами полиномов, применяемых для приближения исходных и возмущенных данных составляет величину порядка 1010,

>> max(abs(p-pe))
ans =
  3.0239e+010

Приближение же полиномом третьей степени даст разницу коэффициентов порядка 10-3, в чем несложно убедиться, присвоив в предыдущем примере переменной n значение 3.

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

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

В качестве примера рассмотрим приближение данных

x = [1 2 2 3 3 3 4];
y = [1 5 4 3 4 3 4];
в смысле наименьших квадратов полиномом четвертой степени
p = polyfit(x, y, 4);
Warning: Polynomial is badly conditioned. Remove repeated data points
         or try centering and scaling as described in HELP POLYFIT.

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

for k=1:length(x) - 1  % проходим по массиву x
    % находим индексы элементов массива x, совпадающих с x(k)
    ind = find(x == x(k));
    if length(ind)>1         % есть элементы, совпадающие с x(k)
        % удаляем их из массива x
        x(ind(2:end))=[];
        % усредняем соответствующие элементы в массиве y
        y(ind(1)) = mean(y(ind));
        % удаляем лишние элементы из массива y
        y(ind(2:end)) = [];
    end
end

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

y = [1   5    5.03   3   2.99   3    4];

заданных в двух наборах точек:

x1 = [1 2  2.01         3    3.01         3.02               4];
x2 = [1 2  2.000001 3    3.000001 3.0000002     4];

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

 

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

x = [51  52  53  54  55  56  57];
y = [1.2 3.4 2.9 4.4 4.5 5.1 4.2];

полиномом четвертой степени получаем предупреждение о плохой обусловленности

>> p = polyfit(x, y, 4)
Warning: Polynomial is badly conditioned. Remove repeated data points
         or try centering and scaling as described in HELP POLYFIT.

хотя степень полинома невелика и в массиве x нет совпадающих или близких точек, в которых определены значения данных, как было описано выше.

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

понимаемое в смысле отношения ее наибольшего сингулярного числа к наименьшему, т.е., где и - максимальные и минимальные числа на диагонали матрицы , образующей сингулярное разложение матрицы с ортогональными матрицами и . Пусть данные заданы в семи равномерно отстоящих друг от друга на отрезке [0, 1] точках , . Исследуем, как зависит число обусловленности матрицы от передвижения отрезка, т.е. как будет вести себя обусловленность метода наименьших квадратов, применяемого для семи равномерно отстоящих друг от друга точках на отрезке с увеличением . Как видно из следующего графика, число обусловленности растет достаточно быстро

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

где и соответственно - среднее значение и среднеквадратичное отклонение, вычисляемые по формулам

 

Это преобразование не нужно делать самому, его выполняет функция polyfit, если обратиться к ней с тремя выходными аргументами. В третьем аргументе вернется вектор из двух компонент, первой из них будет среднее , а второй - среднеквадратичное отклонение . Этот вектор понадобится в дальнейшем при вычислении значений полинома при помощи функции polyval, его следует указывать в качестве четвертого входного аргумента (третий входной аргумент задается, если требуется получить оценку ошибки приближения, что относится к статистическим методам анализа данных, см. пример в справочной системе в разделе MATLAB: Mathematics: Data Analysis and Statistics: Case Study: Curve Fitting). Для того, чтобы пропустить третий входной аргумент функции polyval, будем использовать пустой массив.

В нашем предыдущем примере для приближения данных

x = [51  52  53  54  55  56  57];
y = [1.2 3.4 2.9 4.4 4.5 5.1 4.2];

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

x = [51  52  53  54  55  56  57];
y = [1.2 3.4 2.9 4.4 4.5 5.1 4.2]
[p, S, mu] = polyfit(x, y, 4)
xx = linspace(x(1), x(end), 200);
yy = polyval(p, xx, [], mu);
plot(x, y, 'o', xx, yy)

Применение масштабирования в этом примере существенно улучшило вычислительные свойства метода наименьших квадратов и в ходе работы функции polyfit сообщения о плохой обусловленности выведено не было.

Операции с полиномиальными приближениями

Функция polyfit возвращает вектор коэффициентов полинома, приближающий данные в смысле наименьших квадратов. В этом векторе коэффициенты полинома записываются начиная с коэффициента при старшей степени, т.е. в соответствии с соглашением, принятым в MATLAB. Это позволяет применить ряд функций для дальнейших операций с полученными полиномиальными приближениями. Одна из этих функций polyval, предназначенная для вычисления значений полинома в заданных точках была описана в предыдущих разделах, в том числе и для случая, когда к данным применялось центрирование и масштабирование (см. Задача о приближении данных полиномом в смысле наименьших квадратов и Когда функция polyfit выдает предупреждения, плохая обусловленность, центрирование и масштабирование данных). В MATLAB возможно производить следующие операции над полиномами:

Умножение полиномов и деление с остатком при помощи функций conv и deconv. Пусть, например, требуется найти произведение и отношение следующих полиномов:

и

Задаем их векторами коэффициентов

>> p = [1 0 0 2 -1 5];
>> q = [1 0 -1];

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

>> P = conv(p, q)
P =
     1     0     1    -1     3     1    -5 

т.е.

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

 >> [R, r] = deconv(p, q)
R =
     1     0     3
r =
     0     0     0    -1     8

т.е.

, где , , что несложно проверить перемножением при помощи conv и сложением:

>> conv(R, q) + r
ans =
     1     0     2    -1     5

Получаются в точности коэффициенты исходного полинома. Кстати, функция deconv возвращает коэффициенты остатка от деления в векторе (в нашем примере r), длина которого всегда совпадает с длиной вектора коэффициентов делимого (в нашем примере p), а если степень остатка меньше степени делимого, то вектор с коэффициентами остатка дополняется нулями. Это позволило нам выполнить операцию сложения в выражении conv(R, q) + r без ошибки, поскольку складывались два вектора одинаковой длины.

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

>> d = length(q)-length(p);
>> s = [zeros(1,d) p] + [zeros(1,-d) q]
s =
     1     0     3    -1     4

Аналогично производится и вычитание полиномов. MATLAB является объектно-ориентированной системой, поддерживающей создание классов. Можно было создать собственный класс для полиномов и перегрузить операции сложения и умножения (см. раздел справочной системы MATLAB: Programming: Classes and Objects: Example A Polynomial Class).

Интегрирование и дифференцирование полиномов. Для нахождения коэффициентов полиномов, являющихся производной и неопределенным интегралом от заданного полинома предназначены функции polyder и polyint. Для полинома

заданного вектором коэффициентов

>> p = [1 0 0 2 -1 5];

функция polyder возвращает вектор коэффициентов первой производной:

>> p1 = polyder(p)
p1 =
     4     0     4    -1

т.е.

Кроме этого, функция polyder умеет вычислять производную от произведения и отношения двух полиномов. Если функция polyder вызывается с двумя входными аргументами (векторами с коэффициентами полиномов), то она возвращает производную от произведения полиномов, например:

>> dP = polyder(p, q)
dP =
     6     0     4    -3     6     1

возвращает то же самое, что и

>> polyder(conv(p, q))
ans =
     6     0     4    -3     6     1

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

и

заданных их векторами коэффициентов

>> p = [1 0 0 2 -1 5];
>> q = [1 0 -1];
>> [P,Q] = polyder(p, q)
P =
     2     0    -4     1   -14     1
Q =
     1     0    -2     0     1

возвращается

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

>> I = polyint(p)
I =
    0.2000         0    0.6667   -0.5000    5.0000         0

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

а при обращении

>> I = polyint(p,100)
I =
    0.2000         0    0.6667   -0.5000    5.0000  100.0000

вычисляются коэффициенты полинома

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

>> p = [2 -3 4 4];

возвращается

>> r = roots(p)
r =
   1.0545 + 1.4739i
   1.0545 - 1.4739i
  -0.6090

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

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

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

При помощи функции poly создадим вектор коэффициентов полинома с заданными корнями, например:

>> p = poly([1e-16 1.3 2.3 3.3 4.3 5.3 6.3 7.3 1e16]);

и решим обратную задачу о нахождении его корней.

>> format short e
>> r = roots(p)
r =
  1.0000e+016              
  8.1590e+000 +1.3736e+000i
  8.1590e+000 -1.3736e+000i
  5.2418e+000 +3.1452e+000i
  5.2418e+000 -3.1452e+000i
  1.6492e+000 +2.5042e+000i
  1.6492e+000 -2.5042e+000i
            0              
            0

Результат не соответствует действительности. Однако, если решать эту задачу при помощи функции fzero, входящей в MATLAB, то корни найдутся верно. Такое решение реализовано в функции test с подфункцией polynom, которая вычисляет значение исследуемого полинома.

function test
r = fzero(@polynom, 0)
r1 = fzero(@polynom, 1)
r2 = fzero(@polynom, 2)
r3 = fzero(@polynom, 3)
r4 = fzero(@polynom, 4)
r5 = fzero(@polynom, 5)
r6 = fzero(@polynom, 6)
r7 = fzero(@polynom, 7)
R = fzero(@polynom,[10 1e20])

function y=polynom(x)
p = poly([1e-16 1.3 2.3 3.3 4.3 5.3 6.3 7.3 1e16]);
y = polyval(p,x);

В результате работы функции test все корни [1e-16 1.3 2.3 3.3 4.3 5.3 6.3 7.3 1e16] находятся c хорошей точностью:

>> test
r =
    9.999999548218181e-017
r1 =
    1.299999999999999e+000
r2 =
    2.299999999999998e+000
r3 =
    3.299999999999994e+000
r4 =
    4.300000000000304e+000
r5 =
    5.300000000000961e+000
r6 =
    6.299999999997263e+000
r7 =
    7.300000000000616e+000
R =
    1.000000000000000e+016

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

характеристический полином равен

>> A = [5 2 3; 1 6 1; 2 2 9];
>> c = poly(A)
c =
    1.0000  -20.0000  119.0000 -216.0000

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

>> polyvalm(c, A)
ans =
  1.0e-012 *
    0.0853    0.1137    0.1990
    0.0426    0.0284    0.1137
    0.1421    0.1990    0.3126

Т.е. функция polyvalm находит значение

>> c(1)*A^3+c(2)*A^2+c(3)*A+c(4)*eye(3)
ans =
  1.0e-012 *
    0.0853    0.1705    0.2274
    0.0284   -0.0284    0.0853
    0.1705    0.1705    0.3126

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

Упомянем еще одну функцию residual, которая осуществляет преобразование между представлениями отношения двух полиномов.

Пример. Интегрирование данных и поиск их корней. Рассмотрим простой пример, требуется проинтегрировать табличные данные, найти их корни и отметить их на графике. В качестве данных возьмем таблицу значений функции на отрезке [-1, 8] в равномерно отстоящих точках с шагом 0.5

x = -1:0.5:8;
y = (x-4).^2.*sin(x);

Приблизим наши данные полиномом седьмой степени в смысле наименьших квадратов, используя функцию polyfit

p = polyfit(x, y, 7);

и вычислим значение полинома в 100 промежуточных точках функцией polyval, задав сами промежуточные точки функцией linspace

>> xx = linspace(x(1), x(end), 100)
>> yy = polyval(p, xx);

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

>> plot(x, y, 'o', xx, yy)

Найдем все корни полинома, приближающего наши данные в смысле наименьших квадратов, при помощи функции roots

>> r = roots(p)
r =
   8.6530          
  -2.0428          
   6.3088          
   4.2151 + 0.6798i
   4.2151 - 0.6798i
   3.1963          
  -0.0115

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

R = r((r==conj(r))&(r>=x(1))&(r<=x(end)))
R =
    6.3088
    3.1963
   -0.0115

Отобразим их на графике круглыми маркерами красного цвета

hold on
plot(R, polyval(p, R), 'ro')



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

Система Orphus

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