В этом примере показано, как использовать команды сплайна в Toolbox™ Фитинг кривой (Curve Fitting) для подгонки тензорных сплайнов изделия для двумерного отображения данных с сеткой.
Поскольку панель инструментов «Фитинг кривой» (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');
Обратите внимание, что инструкции
[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 не следовал бы стандартам теории аппроксимации.
Далее выберите порядок сплайнов. ky и последовательность узлов knotsy для направления y
ky = 3; knotsy = augknt([0,.25,.5,.75,1],ky);
а затем получить
sp = spap2(knotsy,ky,y,z);
сплайновая кривая, i-й компонент - аппроксимация кривой (y,z(i,:)) для i=1:I.
В частности,
yy = -.1:.05:1.1; vals = fnval(sp,yy);
создает матрицу vals чей (i,j)-й элемент может быть взят в качестве аппроксимации к значению 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-направления, наиболее заметного вблизи вершин поверхности. Они подтверждают, что мы строим гладкие кривые только в одном направлении.
Чтобы получить реальную поверхность, мы теперь должны пойти на шаг дальше. Рассмотрим коэффициенты coefsy сплайна sp, как получено
coefsy = fnbrk(sp,'c');Абстрактно можно подумать о сплайне sp как функция векторного значения
y |--> sum coefsy(:,r) B_{r,ky}(y)
r
с i-й элемент, 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)-й элемент равен значению B_{q,kx}(xv(i)) в xv(i) из q-й 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) banded, это может быть более эффективным для «больших» 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 в качестве независимой переменной, следовательно, это строки z которые теперь становятся значениями данных. Соответственно, мы должны поставлять z.' (а не zКому spap2, и получить
spb = spap2(knotsx,kx,x,z.');
сплайновое приближение ко всем кривым (x,z(:,j)) для j=1:J. В частности,
valsb = fnval(spb,xv).';
создает матрицу, (i,j)-й элемент может быть взят в качестве аппроксимации к значению 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-направления. Они подтверждают, что мы снова строим гладкие кривые только в одном направлении.
Теперь идет второй шаг, чтобы получить фактическую поверхность.
Давайте coefsx быть коэффициентами для spb, т.е.
coefsx = fnbrk(spb,'c');Абстрактно можно подумать о сплайне spb как функция векторного значения
x |--> sum coefsx(r,:) B_{r,kx}(x)
r
с j-я запись 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 = 1.9984e-15
Таким образом, двухмерное сплайновое приближение
(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.
Как отмечалось ранее, можно осуществить всю только что пройденную конструкцию (двумя способами), используя всего два утверждения, одно для построения аппроксиманта наименьших квадратов, другое для его оценки на прямоугольной сетке.
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').';Сейчас мы готовы оценить интерполятор
ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv).';
и постройте график интерполятора на тонкой сетке.
mesh(xv,yv,ivalues.'); xlabel('x'); ylabel('y'); view(150,50); title('The Spline Interpolant');
Опять же, вышеописанные шаги могут быть выполнены с использованием только двух утверждений, одно для построения интерполятора, другое для его оценки в прямоугольной сетке.
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
Далее вычисляем ошибку интерполятора как приближение к функции Франке.
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