В этом примере показано, как использовать csapi и csape команды от Curve Fitting Toolbox™, чтобы создать кубический сплайн interpolants.
Команда
values = csapi(x,y,xx)
возвращает значения в xx из кубического сплайна interpolant к определенным данным (x,y), с помощью граничного условия не-узла. Этот interpolant является кусочной кубической функцией с последовательностью пропуска x, чьи кубические части объединяются, чтобы сформировать функцию с двумя непрерывными производными. Граничное условие "не-узла" означает, что в первом и последнем внутреннем пропуске даже третья производная непрерывна (до ошибки округления).
Определение только двух точек данных приводит к прямой линии interpolant.
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')
Они похожи на хороший interpolants, но как мы проверяем тот csapi выполняет, как рекламируется?
Мы уже видели тот csapi интерполирует, потому что мы построили точки данных, и interpolant пошел прямо через те точки. Но быть уверенным, что мы получаем кубический сплайн, лучше запускаться с данных из кубического сплайна ожидаемого вида и проверять ли csapi воспроизводит тот кубический сплайн, i.e., отдает тот кубический сплайн, из которого были взяты данные.
Одним простым примером кубической функции сплайна, чтобы проверять по является усеченная третья степень, i.e., функция
где 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 и строим interpolant сверху сплайна в черном цвете.
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- произведенный interpolant на сайтах 0:6 не может совпасть с ним. Например, первый внутренний пропуск сплайна интерполяции не является действительно узлом начиная с csapi использует условие "не-узла", следовательно interpolant имеет три непрерывных производные на том сайте. Это подразумевает, что мы не должны мочь воспроизвести усеченную 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 первый внутренний узел, это не активно для этого interpolant.
Различие является столь же большим как.18, но затухает быстро, когда мы переезжаем от 1. Это иллюстрирует, что интерполяция кубическим сплайном чрезвычайно локальна.
Возможно сохранить интерполирующий кубический сплайн в форме, подходящей для последующей оценки, или для вычисления ее производных, или для других манипуляций. Это сделано путем вызова csapi в форме
pp = csapi(x,y)
который возвращает ppform interpolant. Можно оценить эту форму в некоторых новых точках xx командой
values = fnval(pp,xx)
Можно дифференцировать interpolant командой
dpp = fnder(pp)
или интегрируйте его командой
ipp = fnint(pp)
которые возвращают ppform производной или интеграла, соответственно.
Чтобы показать дифференцирование interpolant, мы строим производную этой усеченной степени
(снова в желтом) и затем, сверху его, производной нашего interpolant к исходной усеченной функции третьей степени (снова в черном цвете).
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')
Вторая производная усеченной степени
График различия между этой функцией и второй производной interpolant к исходной функции показывает, что существуют теперь скачки, но они все еще в ошибке округления.
ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')
Интеграл усеченной степени
График различия между этой функцией и интегралом interpolant к исходной функции снова показывает, что ошибки в ошибке округления.
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 команда предоставляет кубический сплайн interpolant определенным данным. Однако это разрешает различные дополнительные граничные условия. Его самая простая версия,
pp = csape(x,y)
использует Лагранжево граничное условие, которое является общей альтернативой условию не-узла, используемому csapi. csape не делает непосредственно возвращаемых значений interpolant, но только его ppform.
Например, рассмотрите снова интерполяцию к функции
который csapi сбои, чтобы воспроизвести хорошо. Мы строим ошибку не-узла interpolant возвращенный csapi (в черном цвете), наряду с ошибкой interpolant получен из 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
Нет большого различия между двумя interpolants в этом случае.
csape команда также обеспечивает способы задать несколько других типов граничных условий для интерполирующего кубического сплайна. Например, команда
pp = csape(x,y,'variational')
использует так называемые 'естественные' граничные условия. Это означает, что вторая производная является нулем в двух экстремальных пропусках.
Этот шаг показывает, как применить 'естественную' интерполяцию кубическим сплайном к функции
и постройте ошибку. Код ниже вычисляет 'естественный' сплайн interpolant с альтернативным синтаксисом аргумента, который эквивалентен 'variational' аргумент: использование '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);
Затем мы создаем interpolant, указывая, что вторые производные в конечных точках должны быть соответствующими к вторым производным значениям, которые мы только вычислили. Мы делаем это путем обеспечения 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 также спецификация разрешений наклонов конечной точки. Это - зафиксированный (или, завершенное) кубический сплайн interpolant. Оператор
pp = csape(x,[sl,y,sr],'clamped')
создает кубический сплайн interpolant к данным (xY) это также имеет наклонный sl на крайнем левом сайте данных и наклонном sr на самом правом сайте данных.
Это является четным возможное смешать эти условия. Например, наша очень осуществленная усеченная функция степени
имеет наклон 0 в x=0 и вторые производные 30 в x=6 (последний сайт данных).
Поэтому путем соответствия с наклоном в левом конце и искривлении справа, мы не ожидаем ошибки в получившемся interpolant.
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-periodic и имеет значения [0 -1 0 1 0] на сайтах (pi/2)*(-2:2). Различие, между синусоидальной функцией и ее периодическим кубическим сплайном interpolant на этих сайтах, составляет только 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 может быть обработан путем построения interpolant с csape условия стороны по умолчанию, и затем добавляющий к нему соответствующее скалярное кратное interpolant к нулевым значениям и некоторым условиям стороны. Если существует два 'нестандартных' условия стороны, которым удовлетворят, вам, вероятно, придется решить линейную систему 2 на 2 сначала.
Например, предположите, что вы хотите вычислить кубический сплайн interpolant 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;
Чтобы создать interpolant, который удовлетворяет этому особому условию, мы сначала создаем interpolant с граничными условиями по умолчанию
pp1 = csape(x,y);
и первая производная его первой полиномиальной части.
dp1 = fnder(fnbrk(pp1,1));
Кроме того, мы создаем кубический сплайн interpolant, чтобы обнулить значения данных, указывая, что он имеет наклон 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 чем interpolant 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
на третьей производной interpolant (биквадратное удовлетворяет этому условию), затем мы создаем дополнительную кубическую интерполяцию сплайна к нулевым значениям, и с нулевой первой производной в левой конечной точке, следовательно убеждающейся быть независимой от 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