В этом примере показано, как использовать команды сплайна в 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-направлению, наиболее заметные вблизи peaks поверхности. Они подтверждают, что мы строим графики только в одном направлении.
Чтобы получить фактическую поверхность, мы теперь должны пойти на шаг дальше. Рассмотрим коэффициенты 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) ленточные, это может быть более эффективным для «больших» 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.');
a сплайна приближения ко всем кривым (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');
Снова обратите внимание на гребни, на этот раз идущие по направлению Х. Они подтверждают, что, еще раз, мы строим графики только в одном направлении.
Теперь наступает второй шаг, чтобы получить фактическую поверхность.
Позвольте 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.
Как наблюдалось ранее, можно провести всю конструкцию, которую мы только что прошли (двумя способами), используя всего два операторов, одно для построения аппроксимации методом наименьших квадратов, другое для ее оценки в прямоугольном 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').';Теперь мы готовы оценить интерполяцию
ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv).';
и постройте график интерполяции на мелком mesh.
mesh(xv,yv,ivalues.'); xlabel('x'); ylabel('y'); view(150,50); title('The Spline 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
Далее вычисляем ошибку интерполяции как приближение функции Франке.
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