Интерполяция кубических сплайнов

Открыть сценарий в реальном времени

В этом примере показано, как использовать csapi и csape команды из Toolbox™ Фитинг кривой (Curve Fitting) для построения интерполяций кубических сплайнов.

Команда CSAPI

Команда

values = csapi(x,y,xx)

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

Задание только двух точек данных приводит к прямолинейному интерполятору.

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

Figure contains an axes. The axes with title Interpolant to Two Points contains 2 objects of type line.

Указание трех точек данных дает параболу.

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

Figure contains an axes. The axes with title Interpolant to Three Points contains 2 objects of type line.

В более общем случае четыре или более точек данных дают кубический сплайн.

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

Figure contains an axes. The axes with title Cubic Spline Interpolant to Five Points contains 2 objects of type line.

Проверка выходных данных CSAPI

Они выглядят как хорошие интерполяторы, но как мы можем это проверить? csapi выступает в рекламе?

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

Пример: Усеченная функция мощности

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

$f (x) = ((_{} {x-xi}^{)}$ +) 3,

где xi является одним из разрывов, а индекс «+» указывает функцию усечения, обеспечиваемую командой subplus:

help subplus

 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

    Documentation for subplus
       doc subplus

Усеченная 3-я мощность нанесена ниже для конкретного выбора xi = 2. Как и ожидалось, она равна нулю слева от 2 и поднимается как (x-2) ^ 3 справа от 2.

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

Figure contains an axes. The axes contains an object of type line.

Теперь мы интерполируем данный кубический сплайн на участках данных 0:6 и строим график интерполятора поверх сплайна черным цветом.

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

Figure contains an axes. The axes with title Interpolant to ((x-2)_+)^3 contains 3 objects of type line.

Ошибка интерполяции

При сравнении двух функций обычно гораздо информативнее строить их разницу.

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

Figure contains an axes. The axes with title Error in Cubic Spline Interpolation to ((x-2)_+)^3 contains an object of type line.

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

max_y = max(abs(y))

max_y = 64

Усеченная мощность, которую невозможно воспроизвести

В качестве дополнительного теста мы интерполируем усеченную мощность, csapi-производимый интерполятор на участках 0:6 не может совпадать с ним. Например, первый внутренний разрыв интерполирующего сплайна на самом деле не является узлом, поскольку csapi использует условие «не а-узел», следовательно, интерполятор имеет три непрерывных производных в этом месте. Это означает, что мы не должны быть в состоянии воспроизвести усеченную 3-ю мощность, центрированную на этом сайте, так как его третья производная прерывистая на этом сайте.

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

Figure contains an axes. The axes with title Error in Not-a-Knot Interpolant to ((x-1)_+)^3 contains an object of type line.

Поскольку 1 является первым внутренним узлом, он не активен для этого интерполятора.

Разница равна 0,18, но быстро распадается, когда мы уходим от 1. Это показывает, что интерполяция кубических сплайнов является по существу локальной.

Использование ppform вместо значений

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

pp = csapi(x,y)

который возвращает ppform интерполятора. Вы можете оценить эту форму в некоторых новых точках xx командой

values = fnval(pp,xx)

Интерполант можно отличить с помощью команды

dpp = fnder(pp)

или интегрировать его с помощью команды

ipp = fnint(pp)

которые возвращают ppform производной или интеграл соответственно.

Пример: Дифференцирование и интеграция интерполятора

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

$_{f2}^{}' (x) = 3 (_{(}^{} x-2$ ) +) 2,

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

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

Figure contains an axes. The axes with title Derivative of Interpolant to ((x-2)_+)^3 contains 2 objects of type line.

Опять же, более информативное сравнение заключается в построении их различий, и, как и раньше, это не больше, чем ошибка округления.

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

Figure contains an axes. The axes with title Error in Derivative of interpolant to ((x-2)_+)^3 contains an object of type line.

Второй производной усеченной мощности является

$_{f2}^{'}' (x) = 6_{}$ (x-2) +

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

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

Figure contains an axes. The axes with title Error in Second Derivative of Interpolant to ((x-2)_+)^3 contains an object of type line.

Интеграл усеченной мощности равен

$_{F2} (x) = ((_{} {x-2}^{)} +)$ 4/4.

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

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

Figure contains an axes. The axes with title Error in Integral of Interpolant to ((x-2)_+)^3 contains an object of type line.

Команда CSAPE

Как csapi, csape обеспечивает интерполяцию кубического сплайна для заданных данных. Однако он допускает различные дополнительные конечные условия. Его простейшая версия,

pp = csape(x,y)

использует конечное условие Лагранжа, которое является общей альтернативой условию без узла, используемому csapi. csape непосредственно возвращает не значения интерполятора, а только его ppform.

Например, снова рассмотрим интерполяцию к функции

$_{f1} (x) = ((_{} {x-1}^{)}$ +) 3,

который csapi плохо воспроизводится. Строим график ошибки неузлового интерполятора, возвращаемого csapi (черным цветом) вместе с ошибкой интерполятора, полученной из csape (красным цветом).

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

Figure contains an axes. The axes with title Error in Not-a-Knot vs. Lagrange End Conditions contains 2 objects of type line. These objects represent Not-a-Knot, Lagrange.

В этом случае между двумя интерполяторами нет большой разницы.

Другие конечные условия: «естественный» сплайн-интерполятор

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

pp = csape(x,y,'variational')

использует так называемые «естественные» конечные условия. Это означает, что вторая производная равна нулю в двух крайних разрывах.

На этом шаге показано, как применить интерполяцию «натурального» кубического сплайна к функции

$_{f2} (x) = ((_{} {x-2}^{)}$ +) 3,

и постройте график ошибки. Приведенный ниже код вычисляет «естественный» сплайн-интерполятор с альтернативным синтаксисом аргументов, эквивалентным «вариационному» строковому аргументу: использование строки «second» указывает, что csape следует установить для второй производной в крайних узлах данных значение по умолчанию 0.

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

Figure contains an axes. The axes with title Error in 'Natural' Spline Interpolation to ((x-2)_+)^3 contains an object of type line.

Обратите внимание на большую ошибку около правого конца. Это связано с тем, что «естественные» конечные условия неявно настаивают на наличии там производной от нуля секунды.

Другие конечные условия: Назначение вторых производных

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

endcond = 6*subplus(x([1 end])-2);

Затем мы создаем интерполятор, указывая, что вторые производные в конечных точках должны соответствовать только что вычисленным значениям второй производной. Мы делаем это, предоставляя endcond(1) для условия левой конечной точки, и endcond(2) справа вместе со значениями данных.

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

Figure contains an axes. The axes with title Error in Spline Interpolation to ((x-1)_+)^3 When Matching the 2nd Derivative at Ends contains an object of type line.

Другие конечные условия: назначение откосов

csape также позволяет специфицировать откосы конечных точек. Это усеченный (или полный) кубический сплайн-интерполятор. Заявление

pp = csape(x,[sl,y,sr],'clamped')

создает кубический сплайн, интерполированный с данными (x, y), который также имеет наклон sl на крайнем левом участке данных и уклоне sr на самом правом сайте данных.

Другие конечные условия: смешанные конечные условия

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

$_{f1} (x) = ((_{} {x-1}^{)}$ +) 3

имеет уклон 0 на x= 0 и вторая производная 30 приx= 6 (последний узел данных).

Поэтому, сопоставляя наклон на левом конце и кривизну на правом, мы ожидаем отсутствие ошибок в полученном интерполяторе.

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

Figure contains an axes. The axes with title Error in Spline Interpolation to ((x-1)_+)^3 with Mixed End Conditions. contains an object of type line.

Другие конечные условия: Периодические условия

Также можно назначать периодические конечные условия. Например, синусоидальная функция является 2 * pi-периодической и имеет значения[0 -1 0 1 0] на площадках (pi/2)*(-2:2). Разница между синусоидальной функцией и ее периодическим кубическим сплайновым интерполятором в этих местах составляет всего 2 процента. Неплохо.

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

Figure contains an axes. The axes with title Error in Periodic Cubic Spline Interpolation to sin(x) contains an object of type line.

Конечные условия, явно не охватываемые CSAPI или CSAPE

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

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

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

и обеспечить выполнение условия

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

по первой и второй производным от s в точке e.

Данные были сгенерированы из квартального полинома, который удовлетворяет этому боковому условию с определенными параметрами

e = x(1);
a = 2; b = -3; c = 4;

Чтобы создать интерполятор, удовлетворяющий этому конкретному условию, сначала построим интерполятор с конечными условиями по умолчанию.

pp1 = csape(x,y);

и первую производную его первой полиномиальной части.

dp1 = fnder(fnbrk(pp1,1));

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

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

а также построение первой производной ее первой полиномиальной части.

dp0 = fnder(fnbrk(pp0,1));

Затем вычисляем lambda для обоих pp1 и pp0,

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

и построить правильную линейную комбинацию pp1 и pp0 для получения кубического сплайна

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

это соответствует требуемому условию, а также конечному условию по умолчанию в правой конечной точке. Мы формируем эту линейную комбинацию с помощью fncmb.

s = fncmb(pp0,(c-lam1)/lam0,pp1);

График ошибки интерполяции показывает, что s подходит к квартирному многочлену немного лучше вблизи e чем интерполятор pp1 с условиями по умолчанию.

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Default conditions, Nonstandard conditions.

Если мы хотим выполнить условие

mu(s) := (D^3 s)(3) = 14.6

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

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

Затем мы находим коэффициенты d0 и d2 в линейной комбинации

s := pp1 + d0*pp0 + d2*pp2

которая решает линейную систему

lambda(s) = c

mu(s) = 14.6

Обратите внимание, что оба pp0 и pp2 исчезают на всех участках интерполяции, следовательно s будет соответствовать данным для любого выбора d0 и d2.

Для развлечения мы используем средство кодирования MATLAB ® для записи цикла для вычисленияlambda(pp_j) и mu(pp_j), для j=0:2.

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

Учитывая значения lambda и mu для pp0, pp1, и pp2, затем решаем для коэффициентов, которые определяют правильную линейную комбинацию.

d = dd(:,[1,3])\([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

Figure contains an axes. The axes with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains an object of type line.

Для подтверждения сравним эту ошибку с ошибкой, полученной при полной интерполяции кубического сплайна, с этой функцией:

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

Figure contains an axes. The axes with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains 2 objects of type line. These objects represent Nonstandard conditions, Endslope conditions.

Ошибки различаются (и не сильно) только вблизи конечных точек, свидетельствуя о том, что оба pp0 и pp2 значимы только вблизи их соответствующих конечных точек.

В качестве окончательной проверки мы проверяем, что s удовлетворяет требуемому состоянию третьей производной при 3.

fnval(fnder(s,3),3)

ans = 14.6000

Документация