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

Проблема

Вы хотите создать приближение наименьших квадратов к определенным данным (x, y) от пробела S “естественных” кубических сплайнов с данными пропусками b(1) <... <b(l+1).

Общее разрешение

Если вы знаете основание, (f1, f2..., из), для линейного пробела S всех “естественных” кубических сплайнов с последовательностью пропуска b, то вы учились находить приближение наименьших квадратов в c(1) f1 формы + c(2) f2 +... + c(m) из с векторным c решение методом наименьших квадратов к линейной системе A*c = y, матрицей коэффициентов которого дают

A(i,j) = fj(x(i)),   i=1:length(x),  j=1:m .

Другими словами, c = A\y.

Потребность в базисной карте

Общее решение, кажется, требует, чтобы вы знали основание. Однако в порядке создать содействующую последовательность c, только необходимо знать матричный A. Для этого достаточно иметь под рукой базисную карту, а именно, функциональный F говорит, так, чтобы F(c) возвратил сплайн, данный конкретной взвешенной суммой c(1) f1+c(2)f2+... +c(m) из. Поскольку, с этим можно получить, для j=1:m, j-th столбец A как fnval(F(ej),x), с ej j-th столбец eye(m), единичная матрица порядка m.

Еще лучше функции сплайна Curve Fitting Toolbox могут обработать функции с векторным знаком, таким образом, необходимо смочь создать базисную карту F, чтобы обработать коэффициенты с векторным знаком c(i) также. Однако по соглашению, в этом тулбоксе, коэффициент с векторным знаком является вектор-столбцом, следовательно последовательность c является обязательно вектором - строкой из вектор-столбцов, т.е. матрицей. С этим F(eye(m)) является сплайном с векторным знаком, i которого-th компонент является базисным элементом fi, i=1:m. Следовательно, принимая векторный x сайтов данных, чтобы быть вектором - строкой, fnval(F(eye(m)),x) является матрицей, (i,j) которой - запись является значением fi в x(j), т.е. транспонированием матричного A, который вы ищете. С другой стороны, как только указано, ваша базисная карта F ожидает, что содействующая последовательность c будет вектором - строкой, т.е. транспонированием векторного A\y. Следовательно, принятие, соответственно, векторный y значений данных, чтобы быть вектором - строкой, можно получить приближение наименьших квадратов от S до данных (x, y) как

F(y/fnval(F(eye(m)),x))

Безусловно, если бы вы хотели быть готовыми к x и y быть произвольными векторами (той же длины), вы использовали бы вместо этого

F(y(:).'/fnval(F(eye(m)),x(:).'))

Базисная карта для “естественных” кубических сплайнов

Что точно требуется базисной карты F для линейного пробела S “естественных” кубических сплайнов с последовательностью пропуска b(1) < ... < b(l+1)? Принятием размерности этого линейного пробела является m, карта, F должен настроить линейное взаимно-однозначное соответствие между m - векторы и элементами S. Но это точно, что делает csape(b, . ,'var').

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

function s = F(c)
s = csape(b,c,'var');

Для данного векторного c (той же длины как b), это предоставляет уникальному “естественному” кубическому сплайну последовательность пропуска b, который принимает значение c(i) в b(i), i=1:l+1. Уникальность является ключевой. Это гарантирует, что соответствие между векторным c и получившимся сплайном F(c) является непосредственным. В частности, m равняется length(b). Больше, чем это, потому что значение f (t) функционального f в точке t зависит линейно от f, эта уникальность, гарантируют, что F(c) зависит линейно от c (потому что c равняется fnval(F(c),b), и инверсия обратимой линейной карты является снова линейной картой).

Короткое решение

Соединяя все это, вы прибываете в следующий код

csape(b,y(:).'/fnval(csape(b,eye(length(b)),'var'),x(:).'),...
'var')

для приближения наименьших квадратов “естественными” кубическими сплайнами с последовательностью пропуска b.

Потребность в соответствующей экстраполяции

Давайте попробуем его на некоторых данных, данные о переписи, скажем, которые обеспечиваются в MATLAB^® командой

load census

и который предоставляет годы, 1790:10:1990, как cdate и значения как pop. Используйте последовательность пропуска 1810:40:1970.

b = 1810:40:1970;
s = csape(b, ...
pop(:)'/fnval(csape(b,eye(length(b)),'var'),cdate(:)'),'var');
fnplt(s, [1750,2050],2.2);
hold on
plot(cdate,pop,'or');
hold off

Взгляните на Приближение Наименьших квадратов “Естественными” Кубическими Сплайнами С Тремя Внутренними Пропусками, которое показывает, в толстом синем, получившемся приближении, наряду с определенными данными.

Это похоже на хорошее приближение, - за исключением того, что оно не похоже на “естественный” кубический сплайн. “Естественный” кубический сплайн, чтобы вспомнить, должен быть линейным слева от его первого пропуска и справа от его последнего пропуска, и это приближение не удовлетворяет никакое условие. Это происходит из-за следующих фактов.

“Естественный” кубический сплайн interpolant к определенным данным обеспечивается csape в ppform, с интервалом, заполненным по условию, располагает его основной интервал. С другой стороны, оценка ppform вне ее основного интервала сделана, в MATLAB ppval или функция сплайна Curve Fitting Toolbox fnval, при помощи соответствующей полиномиальной части конца ppform, т.е. экстраполяцией полного порядка. В случае “естественного” кубического сплайна вы хотите вместо этого экстраполяцию второго порядка. Это означает, что вы хотите, слева от первого пропуска, прямая линия, которая соглашается с кубическим сплайном в значении и наклоне в первом пропуске. Такая экстраполяция обеспечивается fnxtr. Поскольку “естественный” кубический сплайн имеет нулевую вторую производную в своем первом пропуске, такая экстраполяция является даже третьим порядком, т.е. это удовлетворяет три соответствующих условия. Таким же образом, вне последнего пропуска кубического сплайна, вы хотите прямую линию, которая соглашается со сплайном в значении и наклоне в последнем пропуске, и это также предоставляется fnxtr.

Правильное короткое решение

Следующий короткий код предоставляет правильное приближение наименьших квадратов данным (x, y) “естественными” кубическими сплайнами с последовательностью пропуска b:

fnxtr(csape(b,y(:).'/ ...
    fnval(fnxtr(csape(b,eye(length(b)),'var')),x(:).'),'var'))

Но это - по общему признанию, довольно длинная линия.

Следующий код использует эту правильную формулу и графики, в разбавителе, красной линии, получившемся приближении сверху более ранних графиков, как показано в Приближении Наименьших квадратов “Естественными” Кубическими Сплайнами С Тремя Внутренними Пропусками.

 ss = fnxtr(csape(b,pop(:)'/ ...
  fnval(fnxtr(csape(b,eye(length(b)),'var')),cdate(:)'),'var'));
hold on, fnplt(ss,[1750,2050],1.2,'r'),grid, hold off
legend('incorrect approximation','population', ...
'correct approximation')

Приближение наименьших квадратов кубическими сплайнами

Короткое решение работает отлично, если вы хотите аппроксимировать пробелом S всех кубических сплайнов с данной последовательностью пропуска b. Вы не должны даже использовать функции сплайна Curve Fitting Toolbox для этого, потому что можно полагаться на spline MATLAB. Вы знаете, что, с c последовательность, содержащая еще две записи, чем, делает b, spline(b,c) предоставляет уникальному кубическому сплайну последовательность пропуска b, который принимает значение c(i+1) в b(i), всем i, и берет наклонный c(1) в b(1) и наклонный c(end) в b(end). Следовательно, spline(b,.) является базисной картой для S.

Больше, чем это, вы знаете, что spline(b,c,xi) обеспечивает значение (значения) в xi этого сплайна интерполяции. Наконец, вы знаете, что spline может обработать данные с векторным знаком. Поэтому следующий короткий код создает приближение наименьших квадратов кубическими сплайнами с последовательностью пропуска b к данным (x, y):

spline(b,y(:)'/spline(b,eye(length(b)),x(:)'))