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

Введение

Поскольку Curve Fitting Toolbox может обработать сплайны с векторными коэффициентами, легко реализовать интерполяцию или приближение к данным с координатной сеткой сплайнами продукта тензора. Большинство команд конструкции сплайна в тулбоксе использует в своих интересах это.

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

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

Рассмотрите, например, приближение наименьших квадратов к определенным данным

z(i,j) = f(x(i),y(j)) for i = 1:I, j = 1:J.

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

x = sort([(0:10)/10,.03 .07, .93 .97]);
y = sort([(0:6)/6,.03 .07, .93 .97]);
[xx,yy] = ndgrid(x,y); % note: ndgrid rather than meshgrid
z = franke(xx,yy);

mesh(x,y,z.');
xlabel('x'); ylabel('y');
view(150,50);
title('Data from the Franke Function');

Примечание о NDGRID по сравнению с MESHGRID

Обратите внимание на то, что операторы

[xx,yy] = ndgrid(x,y);

z = franke(xx,yy);

используемый выше удостоверяются, что z(i,j) является значением функции, аппроксимируемой в узле решетки (x(i),y(j)).

Однако mesh(x,y,z) команды MATLAB® ожидает z(j,i) (отметьте обратный порядок i и j) как значение в узле решетки (x(i),y(j)). По этой причине вышеупомянутый график был сгенерирован оператором

mesh(x,y,z.');

т.е. использование транспонирования матричного z.

Такое перемещение не было бы необходимо, имел, мы использовали meshgrid вместо ndgrid. Но получившийся z не следовал бы стандартам теории приближения.

Выбор пробела сплайна в направлении Y

Затем, мы выбираем порядок сплайна ky и последовательность узла knotsy для направления Y

ky = 3;
knotsy = augknt([0,.25,.5,.75,1],ky);

и затем получите

sp = spap2(knotsy,ky,y,z);

сплайновая кривая, i которой-th компонент является приближением к кривой (y,z(i,:)) для i=1:I.

Оценка

В частности,

yy = -.1:.05:1.1;
vals = fnval(sp,yy);

создает матричный vals, (i,j) которого-th элемент может быть взят в качестве приближения к значению f(x(i),yy(j)) базового функционального f в узле решетки (x(i),yy(j)). Это очевидно, когда мы строим vals.

mesh(x,yy,vals.');
xlabel('x'); ylabel('y');
view(150,50);
title('Simultaneous Approximation to All Curves in the Y-Direction');

Обратите внимание на то, что, для каждого x(i), и первые два и последние два значения являются нулем, поскольку и первыми двумя и последними двумя сайтами в yy является вне основного интервала для сплайна sp.

Также отметьте "гребни", которые запускаются вдоль направления Y, самого примечательного около peaks поверхности. Они подтверждают, что мы строим плавные кривые в одном направлении только.

От кривых до поверхности; выбор пробела сплайна в направлении X

Чтобы получить фактическую поверхность, мы теперь должны пойти один шаг вперед. Считайте коэффициенты coefsy сплайна sp, как получено

coefsy = fnbrk(sp,'c');

Абстрактно, можно думать о сплайне sp как функция с векторным знаком

y |--> sum coefsy(:,r) B_{r,ky}(y)

r

с i-th элемент, coefsy(i,r), векторного коэффициента coefsy(:,r), соответствующий x(i) для i=1:I. Это предлагает аппроксимировать каждую кривую (x,coefsy(:,r)) сплайном, с помощью того же порядка kx и той же соответствующей последовательности узла knotsx для каждого r.

kx = 4;
knotsx = augknt(0:.2:1,kx);
sp2 = spap2(knotsx,kx,x,coefsy.');

Для использования команды spap2 здесь нужно, возможно, объяснение.

Вспомните, что spap2(knots,k,x,fx) обрабатывает fx(:,j) как значение в x(j), т.е. берет каждый столбец fx как значение данных. Поскольку мы хотели соответствовать значению coefsy(i,:) в x(i) для всего i, мы должны предоставить spap2 транспонирование coefsy.

Теперь полагайте, что транспонирование матрицы коэффициентов получившегося сплайна "изгибает" sp2, полученный как

coefs = fnbrk(sp2,'c').';

coefs обеспечивает приближение двумерного сплайна

(x,y) |--> sum sum coefs(q,r) B_{q,kx}(x) B_{r,ky}(y)

q r

к исходным данным

(x(i),y(j)) |--> f(x(i),y(j)) = z(i,j).

Мы используем spcol, чтобы обеспечить значения, B_{q,kx}(xv(i)) и B_{r,ky}(yv(j)) должны были оценить эту поверхность сплайна в некоторых узлах решетки (xv(i),yv(j)) и затем построить значения.

xv = 0:.025:1; yv = 0:.025:1;
values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv).';
mesh(xv,yv,values.');
xlabel('x'); ylabel('y');
view(150,50);
title('The Spline Approximant');

Почему эта оценка работает?

Оператор

values = spcol(knotsx,kx,xv) * coefs * spcol(knotsy,ky,yv).'

используемый выше проявляет здравый смысл с тех пор, например, spcol(knotsx,kx,xv) является матрицей, (i,q) которой-th запись равняется значению B_{q,kx}(xv(i)) в xv(i) q-th B-сплайн порядка kx для последовательности узла knotsx, в то время как мы хотим выполнить выражение

sum sum coefs(q,r) B_{q,kx}(x) B_{r,ky}(y)

q r

= sum sum B_{q,kx}(x) coefs(q,r) B_{r,ky}(y)

q r

в (x,y) = (xv(i),yv(j)).

Более эффективные альтернативы

Начиная с матриц соединены spcol(knotsx,kx,xv) и spcol(knotsy,ky,yv), может быть более эффективно для "большого" xv и yv (хотя, возможно, больше потребления памяти) использовать fnval.

value2 = fnval(spmak(knotsx,fnval(spmak(knotsy,coefs),yv).'),xv).';

На самом деле fnval и spmak могут иметь дело непосредственно с многомерными сплайнами, следовательно вышеупомянутый оператор может быть заменен

value3 = fnval( spmak({knotsx,knotsy},coefs), {xv,yv} );

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

value4 = fnval( spap2({knotsx,knotsy},[kx ky],{x,y},z), {xv,yv} );

Проверьте его

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

diffs = abs(values-value2) + abs(values-value3) + abs(values-value4);
max(max(diffs)) / max(max(abs(values)))

ans = 1.1206e-15

Эти четыре метода возвращают те же значения до ошибки округления.

Ошибка приближения

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

errors = z - spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).';
mesh(x,y,errors.');
xlabel('x'); ylabel('y');
view(150,50);
title('Error at the Given Data Sites');

Относительная погрешность

max(max(abs(errors))) / max(max(abs(z)))

ans = 0.0539

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

(degrees of freedom used) / (number of data points)

только

numel(coefs)/numel(z)

ans = 0.2909

Очевидное смещение этого подхода

Подход, сопровождаемый здесь, кажется biased: Мы сначала думаем о значениях определенных данных z как описание функции с векторным знаком y, и затем мы обрабатываем матрицу, сформированную векторными коэффициентами кривой приближения как описание функции с векторным знаком x.

Что происходит, когда мы берем вещи в противоположном порядке, т.е. думаем о z как об описании функции с векторным знаком x, и затем обрабатываем матрицу, составленную от векторных коэффициентов кривой приближения как описание функции с векторным знаком y?

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

Выполнение его наоборот: запустите с пробела сплайна в X

Во-первых, мы соответствуем сплайновой кривой к данным, но на этот раз с x как независимая переменная, следовательно это - строки z, которые теперь становятся значениями данных. Соответственно, мы должны предоставить z.' (а не z) к spap2 и получить

spb = spap2(knotsx,kx,x,z.');

приближение сплайна ко всем кривым (x,z(:,j)) для j=1:J. В частности,

valsb = fnval(spb,xv).';

создает матрицу, (i,j) которой-th элемент может быть взят в качестве приближения к значению f(xv(i),y(j)) базового функционального f в узле решетки (xv(i),y(j)). Это очевидно, когда мы строим valsb.

mesh(xv,y,valsb.');
xlabel('x'); ylabel('y');
view(150,50);
title('Simultaneous Approximation to All Curves in the X-Direction');

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

От кривых до поверхности: использование пространства сплайна в направлении Y

Теперь прибывает второй шаг, чтобы получить фактическую поверхность.

Позвольте coefsx быть коэффициентами для spb, т.е.

coefsx = fnbrk(spb,'c');

Абстрактно, можно думать о сплайне spb как функция с векторным знаком

x |--> sum coefsx(r,:) B_{r,kx}(x)

r

с j-th запись coefsx(r,j) векторного коэффициента coefsx(r,:), соответствующий y(j), для всего j. Таким образом мы теперь соответствуем каждой кривой (y,coefsx(r,:)) сплайном, с помощью того же порядка ky и той же соответствующей последовательности узла knotsy для каждого r.

spb2 = spap2(knotsy,ky,y,coefsx.');

В конструкции spb2 мы снова должны транспонировать матрицу коэффициентов от spb, поскольку spap2 берет столбцы своего последнего входного параметра как значения данных.

Поэтому нет теперь никакой потребности транспонировать матрицу коэффициентов coefsb получившейся "кривой".

coefsb = fnbrk(spb2,'c');

Требование: coefsb равняется более раннему массиву коэффициентов coefs до округления. Для доказательства этого смотрите обсуждение построения продукта тензора в документации Curve Fitting Toolbox. Здесь, мы просто осуществляем следующую проверку.

max(max(abs(coefs - coefsb)))

ans = 9.9920e-16

Таким образом, приближение двумерного сплайна

(x,y) |--> sum sum coefsb(q,r) B_{q,kx}(x) B_{r,ky}(y)

q r

к исходным данным

(x(i),y(j)) |--> f(x(i),y(j)) = z(i,j)

полученный совпадает с более ранним, который сгенерировал coefs, а не coefsb.

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

tsp = spap2({knotsx,knotsy},[kx,ky],{x,y},z);
valuet = fnval(tsp,{xv,yv});

Здесь, как другая проверка, относительная разница между значениями, вычисленными ранее и вычисленные теперь:

max(max(abs(values-valuet))) / max(max(abs(values)))

ans = 3.7353e-16

Другой пример: Интерполяция

Поскольку данные прибывают из сглаженной функции, мы должны интерполировать ее, т.е. использовать spapi вместо spap2, или, эквивалентно, использовать spap2 с соответствующими последовательностями узла. Для рисунка вот тот же процесс, сделанный с spapi.

Чтобы вспомнить, (одномерные) сайты данных были

x

x = 1×15

         0    0.0300    0.0700    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    0.9300    0.9700    1.0000

y

y = 1×11

         0    0.0300    0.0700    0.1667    0.3333    0.5000    0.6667    0.8333    0.9300    0.9700    1.0000

Мы используем снова квадратичные сплайны в y, следовательно используем узлы на полпути между сайтами данных.

knotsy = augknt( [0 1 (y(2:(end-2))+y(3:(end-1)))/2 ], ky);
spi = spapi(knotsy,y,z);
coefsy = fnbrk(spi,'c');

Интерполяция получившихся коэффициентов

Мы используем снова кубические сплайны в x и используем условие не-узла. Мы поэтому используем все кроме второго и предпоследних точек данных как узлы.

knotsx = augknt(x([1,3:(end-2),end]), kx);
spi2 = spapi(knotsx,x,coefsy.');
icoefs = fnbrk(spi2,'c').';

Оценка

Мы теперь готовы оценить interpolant

ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv).';

и постройте interpolant в мелкой сетке.

mesh(xv,yv,ivalues.');
xlabel('x'); ylabel('y');
view(150,50);
title('The Spline Interpolant');

Снова, шаги выше могут быть выполнены с помощью всего два оператора, один для конструкции interpolant, другого для его оценки в прямоугольной mesh.

tsp = spapi({knotsx,knotsy},{x,y},z);
valuet = fnval(tsp,{xv,yv});

Для проверки мы также вычисляем относительную разницу между значениями, вычисленными ранее и вычисленные теперь.

max(max(abs(ivalues-valuet))) / max(max(abs(ivalues)))

ans = 3.6712e-16

Ошибка приближения

Затем, мы вычисляем ошибку interpolant как приближение к функции Франке.

fvalues = franke(repmat(xv.',1,length(yv)),repmat(yv,length(xv),1));
error = fvalues - ivalues;
mesh(xv,yv,error.');
xlabel('x'); ylabel('y');
view(150,50);
title('Interpolation error');

Относительная ошибка приближения

max(max(abs(error))) / max(max(abs(fvalues)))

ans = 0.0409