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

Spline Toolbox

Другие типы интерполяции в MATLAB

Функции, входящие в пакет MATLAB, реализуют следующие способы интерполяции:
  • интерполяция по ближайшим соседям;
  • кусочно-линейная интерполяция;
  • интерполяция кубическими сплайнами;
  • интерполяция кусочными полиномами Эрмита.
Эти возможности предоставляет функция interp1, обращение к которой выглядит следующим образом:
yy = interp1(ksi, y, xx, method)
где ksi - массив узлов, y - массив значений табличной функции, xx - значения абсцисс, для которых требуется вычислить значения интерполянта и записать в массив yy, а входной аргумент method определяет способ интерполяции: 'nearest' (по ближайшим соседям), 'linear' (кусочно-линейная), 'spline' (кубическими сплайнами), 'phcip' или 'cubic' (кусочными полиномами Эрмита). Указание 'spline' или 'phcip' приводит к вызову соответствующих функций pchip или spline, которыми можно пользоваться отдельно без interp1, т.е, следующие обращения приводят к одинаковым результатам:
yy = spline(ksi, y, xx)  yy = interp1(ksi, y ,xx, 'spline')
yy = pchip(ksi, y, xx)  yy = interp1(ksi, y, xx, 'pchip')
Другой способ обращения к этим функциям для построения кубических сплайнов и кусочных полиномов Эрмита
pp = spline(ksi, y)  pp = interp1(ksi, y, 'spline', 'pp')
pp = pchip(ksi, y)  pp = interp1(ksi, y, 'pchip', 'pp')
приводит к получению структуры pp с информацией о сплайне, с которой можно работать функциями Spline ToolBox (например интегрировать при помощи fnint, строить график, используя fnplt и т.д.).

Функция spline строит кубический интерполяционный сплайн с граничными условиями отсутствия узла (см. разд. Функция csapi для интерполяции с условиями отсутствия узла (not-a-knot)) и с точностью до ошибок округления приводит к тому же результату, что и spline. Для задания значений производной в граничных узлах следует дополнить вектор y значениями производной в левом и правом узлах, которые должны являться первым и последним элементами нового вектора, соответственно. Например:
pp = spline(ksi, [-3 y 4])
приводит к построению кубического сплайна f(x) такого, что и . Для функции pchip такого варианта вызова не предусмотрено потому, что при построении кусочного полинома Эрмита, значения его первой производной во всех узлах определяются автоматически.

Интерполяция кусочными полиномами Эрмита

При интерполяции кубическими сплайнами в узлах накладываются условия непрерывности на сплайн, его первую и вторую производные. Однако, часто имеет смысл отказаться от условий непрерывности второй производной и заменить их на другие: в узлах интерполяции первая производная принимает заданные значения. Такие кусочные полиномы называются полиномами Эрмита. Если значения производной интерполируемой функции известны, то имеет смысл их указать. Однако, при интерполяции табличных данных значения первых производных в узлах могут быть неизвестны, тогда необходимо их определить некоторым способом.

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


то в каждом k-ом внутреннем узле значения производной кусочного полинома f(x) определяют следующим образом:
  1. если знаки величин и различны, или одна из них равна нулю, то полагают

    ;

  2. если знаки и одинаковы и , то прибегают к усреднению

    ;

  3. если знаки и одинаковы и , то прибегают к усреднению

    .

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


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


где . Значение в правой граничной точке определяется по похожей формуле.

Интерполяция кусочными полиномами Эрмита может быть выполнена при помощи функции MATLAB, которая называется pchip (от piecewise cubic Hermite interpolating polynomial), функции interp1, либо средствами Curve Fitting Toolbox. Сравним интерполяцию кусочными полиномами Эрмита c интерполяцией кубическими сплайнами, выполненной функцией csapi, которую мы рассматривали выше:
% задание табличной функции
ksi = -1:0.2:1;
y = sin(15*ksi).*exp(-ksi.^2);
% построение данных маркерами
plot(ksi, y, 'Marker', 'o', 'MarkerSize', 10, 'LineStyle', 'none')
hold on
% задание абсцисс точек, в которых надо вычислить кусочный полином Эрмита
xx = -1:0.01:1;
% вычисление кусочного полинома Эрмита
yy = pchip(ksi, y, xx);
% построение графика кусочного полинома Эрмита
plot(xx, yy, 'LineWidth', 2, 'Color', 'r')
% конструирование кубического сплайна
pp = csapi(ksi,y)
% построение графика кубического сплайна
fnplt(pp, 'k')
% добавление легенды на график
legend('data', 'pchip', 'csapi', -1)
Результат интерполирования приведен ниже на рисунке, из которого видно, что кусочный полином Эрмита не превосходит значений табличной функции в отличие от кубического интерполяционного сплайна.


Вместо вычисления кусочного полинома Эрмита в точках с абсциссами xx можно было сразу получить структуру, содержащую информацию о нем в формате, принятом в Spline Toolbox. Для этого следовало обратиться к pchip так:
>> pp = pchip(ksi, y)
pp = 
      form: 'pp'
    breaks: [-1 -0.8000 -0.6000 -0.4000 -0.2000 0 0.2000 0.4000 0.6000 0.8000 1]
     coefs: [10x4 double]
    pieces: 10
     order: 4
       dim: 1
или к interp1:
pc=interp1(ksi, y, 'pchip', 'pp')
Проверим, что кусочный полином Эрмита имеет непрерывные первые и разрывные вторые производные. Для этого построим полином при помощи функции interp1 и запишем информацию о нем в структуру, затем вычислим его первую и вторую производные, обратившись к функции fnder Spline Toolbox и построим их графики, воспользовавшись функцией fnplt, так же входящей в состав Spline Toolbox (см. рис).
% задание табличной функции
ksi = [0 1 3 5 7 9 11];
y = [-1 2 1 4 3 0 -2];
% нахождение кусочного полинома Эрмита
pp = interp1(ksi, y, 'pchip', 'pp')
% построение графика кусочного полинома Эрмита
subplot(3, 1, 1)
fnplt(pp)
title('piecewise cubic Hermite interpolant')
% построение графика первой производной кусочного полинома Эрмита
subplot(3, 1, 2)
fnplt(fnder(pp))
title('its first derivative')
% построение графика второй производной кусочного полинома Эрмита
subplot(3, 1, 3)
fnplt(fnder(pp, 2))
title('its second derivative')


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

Система Orphus

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