В этом примере показано, как использовать csapi и csape команды из Curve Fitting Toolbox™ для создания кубических сплайн интерполяций.
Команда
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')
Установка трех точек данных дает параболу.
x = [2 3 5]; y = [1 0 4]; plot(xx,csapi(x,y,xx),'k-',x,y,'ro') title('Interpolant to Three Points')
В более общем случае четыре или более точек данных дают кубический сплайн.
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')
Они выглядят как хорошие интерполяции, но как мы проверяем это csapi выступает как объявлено?
Мы уже видели это csapi интерполирует, потому что мы построили графики точек данных, и интерполянт прошел прямо через эти точки. Но чтобы быть уверенными, что мы получаем кубический сплайн, лучше всего начать с данных кубического сплайна ожидаемой сортировки и проверить, csapi ли воспроизводит тот кубический сплайн, т.е. возвращает тот кубический сплайн, из которого были взяты данные.
Одним из простых примеров кубической сплайн для проверки является усеченная третья степень, т.е. функция
где 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])
Теперь мы интерполируем этот конкретный кубический сплайн в местах данных 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')
При сравнении двух функций обычно гораздо информативнее построить график их различий.
plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')
Чтобы поместить размер их различий в контекст, можно также вычислить максимальное значение данных. Это показывает, что ошибка не хуже неизбежной ошибки округления.
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')
Поскольку 1 является первым внутренним узлом, он не является активным для этой интерполяции.
Размер различия равен 18, но быстро затухает, когда мы отходим от 1. Это иллюстрирует, что кубическая сплайн интерполяция по существу локальна.
Возможно сохранить интерполирующий кубический сплайн в форме, подходящей для последующей оценки, или для вычисления его производных, или для других манипуляций. Это делается по телефону csapi в форме
pp = csapi(x,y)
который возвращает ppform интерполяции. Вы можете оценить эту форму в некоторых новых точках xx по команде
values = fnval(pp,xx)
Вы можете дифференцировать интерполяцию по команде
dpp = fnder(pp)
или интегрировать его по команде
ipp = fnint(pp)
которые возвращают ppform производной или интеграл, соответственно.
Чтобы показать дифференциацию интерполяции, мы построим производную этой усеченной степени
(снова в желтом), а затем, поверх него, производная нашей интерполяции от исходной усеченной третьей степени (снова в черном).
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')
Снова, более информативным сравнением является построение графика их различий, и, как и прежде, это не больше, чем ошибка округления.
plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')
Вторая производная усеченной степени
График различия между этой функцией и второй производной интерполяции исходной функции показывает, что теперь есть переходы, но они все еще находятся в пределах округления.
ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')
Интеграл усеченной степени
График различия между этой функцией и интегралом интерполяции в исходную функцию снова показывает, что ошибки находятся в пределах ошибки округления.
ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')
Как csapi, а csape команда предоставляет кубическую сплайн интерполяцию к данным. Однако он допускает различные дополнительные граничные условия. Его простейший вариант,
pp = csape(x,y)
использует граничное условие Лагранжа, которое является общей альтернативой условию не узла, используемому csapi. csape не возвращает непосредственно значения интерполяции, а только ее ppform.
Например, рассмотрите снова интерполяцию в функцию
который 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
В этом случае нет большого различия между двумя интерполянтами.
The csape команда также предоставляет способы задать несколько других типов граничных условий для интерполирующего кубического сплайна. Для примера используйте команду
pp = csape(x,y,'variational')
использует так называемые «естественные» граничные условия. Это означает, что вторая производная равна нулю на двух крайних пропусках.
Этот шаг показывает, как применить 'естественную' кубическую сплайн интерполяцию к функции
и постройте график ошибки. Приведенный ниже код вычисляет 'natural' сплайн интерполяцию с альтернативным синтаксисом аргумента, который эквивалентен 'variational' string аргументу: использование строки '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')
Обратите внимание на большую ошибку около правого конца. Это связано с тем, что 'естественные' граничные условия неявно настаивают на наличии там нулевой второй производной.
Мы также можем явным образом использовать правильные вторые производные, чтобы получить небольшую ошибку. Во-первых, мы вычисляем правильные вторые значения производной усеченной степени в конечных точках.
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 '])
csape также допускает спецификацию уклонов оконечных точек. Это зажатая (или полная) кубическая сплайн интерполяция. Оператор
pp = csape(x,[sl,y,sr],'clamped')
создает кубическую сплайн интерполяцию к данным (x, y), который также имеет уклон sl на крайнем левом узле данных и уклоне sr на самом правом сайте данных.
Можно даже смешать эти условия. Для примера, наша широко реализованная усеченная функция степени
имеет уклон 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. '])
Также возможно назначать периодические граничные условия. Например, функция синуса является 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)')
Любое граничное условие, не покрытое явно 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')
Если мы хотим применить условие
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))')
Для заверения, мы сравниваем эту ошибку с той, которая получена в полной кубической сплайн интерполяции к этой функции:
hold on plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o') hold off legend({'Nonstandard conditions' 'Endslope conditions'})
Ошибки различаются (и не сильно) только вблизи конечных точек, свидетельствуя о том, что обе pp0 и pp2 являются значительными только вблизи их соответствующих конечных точек.
В качестве окончательной проверки мы проверяем, что s удовлетворяет желаемому третьему условию производной при 3.
fnval(fnder(s,3),3)
ans = 14.6000