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

Команда CSAPI

Команда

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')

Как проверять Вывод CSAPI

Они похожи на хороший interpolants, но как мы проверяем, что csapi выполняет, как рекламируется?

Мы уже видели, что csapi интерполирует, потому что мы построили точки данных, и interpolant пошел прямо через те точки. Но быть уверенным, что мы получаем кубический сплайн, лучше запускаться с данных из кубического сплайна ожидаемого вида и проверять, воспроизводит ли 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.

Усеченная 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. Это иллюстрирует, что интерполяция кубическим сплайном чрезвычайно локальна.

Используя ppform Вместо Значений

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

pp = csapi(x,y)

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

values = fnval(pp,xx)

Можно дифференцировать interpolant командой

dpp = fnder(pp)

или интегрируйте его командой

ipp = fnint(pp)

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

Пример: дифференциация и интеграция Interpolant

Чтобы показать дифференцирование 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')

Команда CSAPE

Как 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 в этом случае.

Другие граничные условия: 'естественный' сплайн Interpolant

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

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

использует так называемые 'естественные' граничные условия. Это означает, что вторая производная является нулем в двух экстремальных пропусках.

Этот шаг показывает, как применить 'естественную' интерполяцию кубическим сплайном к функции

и постройте ошибку. Код ниже вычисляет 'естественный' сплайн interpolant с альтернативным синтаксисом аргумента, который эквивалентен 'вариационному' аргументу строки: использование 'второй' строки указывает, что 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 к данным (x, y), который также имеет наклонный 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

Любое граничное условие, не покрытое явным образом 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