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

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

p = polyfit (x, y, n)

[p, s] = polyfit (x, y, n)

[p, s, mu] = polyfit (x, y, n)

Возвращает коэффициенты полинома p(x) степени n, который минимизирует ошибку наименьших квадратов подгонки к точкам [x, y].

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

Коэффициенты полинома возвращаются в строковом векторе.

Необязательный выход s — это структура, содержащая следующие поля:

‘R’

Треугольный множитель R из QR-разложения.

‘X’

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

‘C’

Немасштабированная матрица ковариации, формально равная обратной к x’*x, но вычисленная таким образом, чтобы минимизировать распространение ошибок округления.

‘df’

Число степеней свободы.

‘normr’

Норма остатков.

‘yf’

Значения полинома для каждого значения x.

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

sqrt (diag (s.C)/s.df)*s.normr.

Когда присутствует третий выходной параметр mu, коэффициенты p связаны с полиномом

xhat = (x - mu(1)) / mu(2)
где mu(1) = среднее значение (x), а mu(2) = стандартное отклонение (x).

Это линейное преобразование x улучшает числовую устойчивость подгонки.

См. также: polyval, polyaffine, roots, vander, zscore.

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

pp = splinefit (x, y, breaks)

pp = splinefit (x, y, p)

pp = splinefit (…, "periodic", periodic)

pp = splinefit (…, "robust", robust)

pp = splinefit (…, "beta", beta)

pp = splinefit (…, "order", order)

pp = splinefit (…, "constraints", constraints)

Подгоняет кусочно-кубический сплайн с разрывами (узлами) breaks к шумным данным x и y.

x — вектор, а y — вектор или N-мерный массив. Если y является N-мерным массивом, то x(j) соответствует y(:,…,:,j).

p — положительное целое число, определяющее количество интервалов вдоль x, а p+1 — количество разрывов. Количество точек в каждом интервале отличается не более чем на 1.

Необязательное свойство periodic — логическое значение, которое указывает, применяется ли периодическое граничное условие к сплайну. Длина периода равна max (breaks) - min (breaks). Значение по умолчанию — false.

Необязательное свойство robust — логическое значение, которое указывает, применяется ли робастная подгонка для уменьшения влияния выбросов. Выполняется три итерации взвешенных наименьших квадратов. Веса вычисляются из предыдущих остатков. Чувствительность определения выбросов регулируется свойством beta. Значение beta ограничено диапазоном 0 < beta < 1. Значение по умолчанию — beta = 1/2. Значения, близкие к 0, придают всем данным одинаковый вес. Увеличение значения beta уменьшает влияние выбросов. Значения, близкие к единице, могут вызвать неустойчивость или недостаточность ранга.

Подгоняемый сплайн возвращается как кусочно-полиномиальная функция pp и может быть вычислен с помощью ppval.

Сплайны строятся из полиномов степени order. Значение по умолчанию — кубический, order=3. Сплайн с P частями имеет P+order степеней свободы. При периодических граничных условиях степени свободы уменьшаются до P.

Необязательное свойство constraints — это структура, определяющая линейные ограничения на подгонку. Структура имеет три поля, "xc", "yc", и "cc".

"xc"

Вектор расположений ограничений по оси x.

"yc"

Значения ограничений в позициях xc. По умолчанию массив нулей.

"cc"

Коэффициенты (матрица). По умолчанию массив единиц. Количество строк ограничено степенью кусочно-полиномиальных функций, order.

Ограничения являются линейными комбинациями производных порядка от 0 до order-1 по формуле

cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).

См. также: interp1, unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.

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

x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## Uniform breaks
breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
pp1 = splinefit (x, y, breaks);
## Breaks interpolated from data
pp2 = splinefit (x, y, 10);  % 11 breaks, 10 pieces
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})

Результат которого можно увидеть на рисунке 28.1.

Рисунок 28.1: Сравнение подгонки кусочно-полиномиальной функции с 41 разрывом с подгонкой с 11 разрывами. Подгонка с большим количеством разрывов демонстрирует быстрое колебание, отсутствующее в исходной функции.

Подгонка кусочно-полиномиальной функции, предоставляемая splinefit, имеет непрерывные производные до order-1. Например, кубическая подгонка имеет непрерывные первую и вторую производные. Это продемонстрировано в коде

## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## Piecewise constant
pp1 = splinefit (x, y, 8, "order", 0);
## Piecewise linear
pp2 = splinefit (x, y, 8, "order", 1);
## Piecewise quadratic
pp3 = splinefit (x, y, 8, "order", 2);
## Piecewise cubic
pp4 = splinefit (x, y, 8, "order", 3);
## Piecewise quartic
pp5 = splinefit (x, y, 8, "order", 4);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
y4 = ppval (pp4, xx);
y5 = ppval (pp5, xx);
plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})

Результат которого можно увидеть на рисунке 28.2.

Рисунок 28.2: Сравнение кусочно-постоянных, линейных, квадратичных, кубических и четырёхчленных полиномов с 8 разрывами с шумными данными. Решения более высокой степени точнее представляют исходную функцию, но требуют большей вычислительной сложности.

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

## Data (100 points)
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## No constraints
pp1 = splinefit (x, y, 10, "order", 5);
## Periodic boundaries
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})

Результат которого можно увидеть на рисунке 28.3.

Рисунок 28.3: Сравнение подгонки кусочно-полиномиальных функций к периодической шумной функции с и без периодических граничных условий.

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

## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## Breaks
breaks = linspace (0, 2 * pi, 10);
## Clamped endpoints, y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## Hinged periodic endpoints, y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})

Результат которого можно увидеть на рисунке 28.4.

Рисунок 28.4: Сравнение двух периодических кусочно-кубических подгонок к периодическому шумному сигналу. Одна подгонка имеет закреплённые концы, а вторая — шарнирно закреплённые.

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

## Data
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## Add outliers
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## Fit splines with hinged conditions
con = struct ("xc", [0, 2*pi]);
## Robust fitting, beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## Robust fitting, beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## No robust fitting
pp3 = splinefit (x, y, 8, "constraints", con);
## Plot
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
         "robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto

Результат которого можно увидеть на рисунке 28.5.

Рисунок 28.5: Сравнение двух уровней робастной подгонки (beta = 0,25 и 0,75) к шумным данным с выбросами. Также включена стандартная подгонка без робастной подгонки (beta = 0).

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

[num, den] = padecoef (T)

[num, den] = padecoef (T, N)

Вычисляет приближение Паде N-го порядка для непрерывной задержки T в форме передаточной функции.

Приближение Паде для exp (-sT) определяется следующим уравнением

Pn(s)
exp (-sT) ~ -------
             Qn(s)

Где как Pn(s), так и Qn(s) являются рациональными функциями N-го порядка, определёнными следующими выражениями

N    (2N - k)!N!        k
Pn(s) = SUM --------------- (-sT)
        k=0 (2N)!k!(N - k)!

Qn(s) = Pn(-s)

Входы T и N должны быть неотрицательными числовыми скалярами. Если N не указан, он по умолчанию равен 1.

Выходные строковые векторы num и den содержат коэффициенты числителя и знаменателя в степенях s в порядке убывания. Оба являются полиномами N-го порядка.

Например:

t = 0.1;
n = 4;
[num, den] = padecoef (t, n)
⇒ num =

      1.0000e-04  -2.0000e-02   1.8000e+00  -8.4000e+01   1.6800e+03

⇒ den =

      1.0000e-04   2.0000e-02   1.8000e+00   8.4000e+01   1.6800e+03

Функция ppval, вычисляет значения кусочно-полиномиальных функций, созданных с помощью mkpp или другими способами, и unmkpp возвращает подробную информацию о кусочно-полиномиальной функции.

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

x = [-2, -1, 1, 2];
p = [ 0,  1, 0;
      1, -2, 1;
      0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);

pp = mkpp (breaks, coefs)

pp = mkpp (breaks, coefs, d)

Создать структуру кусково-полиномиальной функции (pp) из точек разрыва breaks и коэффициентов coefs.

breaks должен быть вектором строго возрастающих значений. Количество интервалов задаётся значением ni = length (breaks) - 1.

Когда m — порядок полинома, coefs должен иметь размерность: ni-by-(m + 1).

i-я строка coefs, coefs(i,:), содержит коэффициенты полинома на i-м интервале, упорядоченные от старшего (m) к младшему (0) степеням.

coefs также может быть многомерным массивом, задающим векторно-значный или многомерно-значный полином. В этом случае порядок полинома m определяется длиной последнего измерения coefs. Размерность первых измерений задаётся скалярным или векторным значением d. Если d не указан, он устанавливается в значение 1. В этом случае p(r, i, :) содержит коэффициенты r-го полинома, определённого на i-м интервале. В любом случае coefs преобразуется в двумерную матрицу размера [ni*prod(d) m].

Примечание по программированию: ppval вычисляет полиномы в точке xi - breaks(i), т.е. вычитает нижнюю границу текущего интервала из xi. Это необходимо учитывать при создании объектов кусково-полиномиальных функций с mkpp.

См. также: unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.

[x, p, n, k, d] = unmkpp (pp)

Извлечь компоненты структуры кусково-полиномиальной функции pp.

Эта функция является обратной к mkpp: она извлекает входные данные для mkpp необходимые для создания структуры кусково-полиномиальной функции PP. Следующий код делает эту связь явной:

[breaks, coefs, numinter, order, dim] = unmkpp (pp);
pp2  = mkpp (breaks, coefs, dim);

Структура кусково-полиномиальной функции pp2 полученная таким образом, идентична исходной pp. То же самое можно получить, напрямую обратившись к полям структуры pp.

Компоненты:

x

Точки разрыва.

p

Коэффициенты полинома для точек в интервале. p(i, :) содержит коэффициенты полинома на интервале i, упорядоченные от старшей к младшей степени. Если d > 1, тогда p является матрицей размера [n*prod(d) m], где i + (1:d) строки — это коэффициенты всех d полиномов на интервале i.

n

Количество кусков полинома или интервалов, n = length (x) - 1.

k

Порядок полинома плюс 1.

d

Количество полиномов, определённых для каждого интервала.

См. также: mkpp, ppval, spline, pchip.

yi = ppval (pp, xi)

Вычислить значения кусково-полиномиальной функции pp в точках xi.

Если pp описывает скалярную полиномиальную функцию, результат представляет собой массив с той же формой, что и xi. В противном случае размер результата равен [pp.dim, length(xi)] если xi является вектором, или [pp.dim, size(xi)] если это многомерный массив.

См. также: mkpp, unmkpp, spline, pchip.

ppd = ppder (pp)

ppd = ppder (pp, m)

Вычислить кусково m-ю производную кусково-полиномиальной структуры pp.

Если m опущено, вычисляется первая производная.

См. также: mkpp, ppval, ppint.

ppi = ppint (pp)

ppi = ppint (pp, c)

Вычислить интеграл кусково-полиномиальной структуры pp.

c, если задано, является константой интегрирования.

См. также: mkpp, ppval, ppder.

jumps = ppjumps (pp)

Вычислить скачки на границах кусково-полиномиальной функции.

Если имеется n интервалов и размерность pp равна d, размерность результирующего массива равна [d, n-1].

См. также: mkpp.