В этом примере показано, как использовать команды сплайна в Curve Fitting Toolbox™, чтобы соответствовать, продукт тензора шлицует к двумерным данным с координатной сеткой.
Поскольку 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- компонент 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 поверхности. Они подтверждают, что мы строим плавные кривые в одном направлении только.
Чтобы получить фактическую поверхность, мы теперь должны пойти один шаг вперед. Считайте коэффициенты 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 как независимая переменная, следовательно это - строки 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. Они подтверждают, что еще раз мы строим плавные кривые в одном направлении только.
Теперь прибывает второй шаг, чтобы получить фактическую поверхность.
Позвольте 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 = 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').';Мы теперь готовы оценить 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