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 обычно является целым числом ≥ 0, определяющим степень аппроксимирующего полинома. Если n является логическим вектором, он используется в качестве маски для выборочного принудительного использования или игнорирования соответствующих коэффициентов полинома.

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

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

‘yf’

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

‘X’

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

‘R’

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

‘C’

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

‘df’

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

‘normr’

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

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

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

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

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

Пример 1: логическое n и целое n

f = @(x) x.^2 + 5;   # data-generating function
x = 0:5;
y = f (x);
## Fit data to polynomial A*x^3 + B*x^1
p = polyfit (x, y, logical ([1, 0, 1, 0]))
⇒ p = [ 0.0680, 0, 4.2444, 0 ]
## Fit data to polynomial using all terms up to x^3
p = polyfit (x, y, 3)
⇒ p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]

См. также: 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.

Кусочно-полиномиальная аппроксимация, предоставляемая 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.

Когда подлежащая функция для подгонки является периодической, 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.

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

## 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.

Функция 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.

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

: [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 преобразуется в матрицу 2-D размера [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.