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) = mean (x), а mu(2) = std (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)

Аппроксимирует зашумлённые данные x и y кусочным кубическим сплайном с узлами breaks.

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-на-(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.