Введение

Поскольку 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)).

Однако команда MATLAB® mesh(x,y,z) ожидает 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- B-сплайн th порядка 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

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

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

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

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 = 7.7716e-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