Приближение сплайнами продукта Tensor

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

Безусловно, большая часть приближения сплайна продукта тензора к данным с координатной сеткой может быть получена непосредственно с одной из команд конструкции сплайна, как spapi или csape, в этом тулбоксе, без беспокойства о деталях обсужден в этом примере. Скорее этот пример предназначается, чтобы проиллюстрировать теорию позади конструкции продукта тензора, и это поможет в ситуациях, не покрытых командами конструкции в этом тулбоксе.

В этом разделе рассматриваются эти аспекты проблемы сплайна продукта тензора:

Выбор сайтов и узлов

Считайте, например, приближение наименьших квадратов к определенным данным z (i, j) =f (x (i), y (j)), i =1:Nx, j =1:Ny. Вы берете данные из функции, используемой экстенсивно Франке для тестирования схем подбора кривой поверхности (см. Р. Франке, “Критическое сравнение некоторых методов для интерполяции данных, имеющий разброс”, Высшая школа ВМС США Techn. Член палаты представителей NPS 53-79-003, март 1979). Его область является модульным квадратом. Вы выбираете еще несколько сайтов данных в x - направлении, чем y - направление; также, для лучшего определения, вы используете более высокую плотность данных около контура.

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

Приближение Наименьших квадратов как Функция y

Обработайте эти данные как прибывающий из функции с векторным знаком, а именно, функции y, значением которого в y (j) является векторный z (: j), весь j. Ни по какой конкретной причине примите решение аппроксимировать эту функцию параболическим сплайном с векторным знаком с тремя однородно расположенными с интервалами внутренними узлами. Это означает, что вы выбираете порядок сплайна и последовательность узла для этого сплайна с векторным знаком как

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

и затем используйте spap2 предоставлять аппроксимирующую функцию наименьших квадратов данным:

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

В действительности вы находите одновременно дискретное приближение наименьших квадратов от S _ky,knotsy к каждому Nx наборы данных

$\begin{matrix} {(y (j), z (i, j))}_{j = 1}^{N y}, & i = 1 : N x \end{matrix}$

В частности, операторы

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

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

Это очевидно в следующем рисунке, полученное командой:

mesh(x,yy,vals.'), view(150,50)

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

Семейство плавных кривых, симулирующих быть поверхностью

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

Отметьте также гребни. Они подтверждают, что вы строите плавные кривые в одном направлении только.

Приближение к Коэффициентам как Функции x

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

coefsy = fnbrk(sp,'coefs');

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

$y | \to \sum_{r} c o e f s y (:, r) B_{r, k y} (y)$

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

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

Обратите внимание на то, что spap2(knots,k,x,fx) ожидает fxJ) быть данной величиной в x(j), т.е. ожидает каждый столбец fx быть значением функции. Соответствовать данной величине coefsy(q, :) в x (q), для всего q, существовать spap2 с транспонированием coefsy.

Двумерное приближение

Теперь рассмотрите транспонирование коэффициентов cxy из получившейся сплайновой кривой:

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

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

$(x, y) | \to \sum_{q} \sum_{r} c o e f s (q, r) B_{q, k x} (x) B_{r, k y} (y)$

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

$(x (i), y (j)) | \to z (x (i), y (j)), i = 1 : N x, j = 1 : N y$

Чтобы построить этот сплайн появляются по сетке, e.g., сетка

xv = 0:.025:1; yv = 0:.025:1;

можно сделать следующее:

values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv).'; 
mesh(xv,yv,values.'), view(150,50);

Это заканчивается в следующем рисунке.

Приближение сплайна к функции Франке

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

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

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

Это, на самом деле, что происходит внутренне когда fnval называется непосредственно со сплайном продукта тензора, как в

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

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

errors = z - spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).'; 
disp( max(max(abs(errors)))/max(max(abs(z))) )

Выходом является 0.0539, возможно, не слишком впечатляющий. Однако массив коэффициентов только имел размер 8 6

disp(size(coefs))

соответствовать массиву данных размера 15 11.

disp(size(z))

Переключитесь в порядок

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

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

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

Приближение Наименьших квадратов как Функция x

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

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

таким образом получая приближение сплайна ко всем кривым (x; z: J. В частности, оператор

valsb = fnval(spb,xv).';

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

mesh(xv,y,valsb.'), view(150,50)

Другое семейство плавных кривых, симулирующих быть поверхностью

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

Приближение к Коэффициентам как Функции y

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

coefsx = fnbrk(spb,'coefs');

Затем соответствуйте каждому вектору коэффициентов coefsx(r,:) сплайном того же порядка ky и с той же соответствующей последовательностью узла knotsy:

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

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

Соответственно, нет теперь никакой потребности транспонировать массив коэффициентов coefsb из получившейся кривой:

coefsb = fnbrk(spb2,'coefs');

Двумерное приближение

Требованием является тот coefsb равняется более раннему массиву коэффициентов coefs, до округления, и вот тест:

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

Выходом является 1.4433e-15.

Объяснение достаточно просто: коэффициенты c из сплайна s содержится в sp = spap2(knots,k,x,y) зависьте линейно от входных значений y. Это подразумевает, учитывая, что оба c и y матрицы с 1 строкой, что существует некоторый матричный A =Aknots_{, k, x} так, чтобы

$c = y A_{knots, k, x}$

для любых данных y. Этот оператор даже содержит когда y матрица, размера d-by-N, скажем, в этом случае каждая данная величина y (: j), взят, чтобы быть точкой в R^d, и получившимся сплайном является соответственно d - с векторным знаком, следовательно его массив коэффициентов c имеет размер d-by-n, с n = length(knots)-k.

В частности, операторы

sp = spap2(knotsy,ky,y,z); 
coefsy =fnbrk(sp,'coefs');

предоставьте нам матричный coefsy это удовлетворяет

$coefs y = z . A_{knotsy,ky,y}$

Последующие расчеты

 sp2 = spap2(knotsx,kx,x,coefsy.'); 
 coefs = fnbrk(sp2,'coefs').';

сгенерируйте массив коэффициентов coefs, которому, учитывая эти два перемещения, удовлетворяет

$\begin{matrix} coefs = ((z A_{knotsy,ky,y})' . A_{knotsx,kx,x})' \\ = (A_{knotsx,kx,x})' . z . A_{knotsy,ky,y} \end{matrix}$

Во втором, альтернативном, вычислении, вы сначала вычислили

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

следовательно coefsx=z _{'.Aknotsx, kx, x}. Последующее вычисление

 spb2 = spap2(knotsy,ky,y,coefsx.');
 coefsb = fnbrk(spb,'coefs');

затем, если

$coefsb = coefsx .' . A_{knotsy,ky,y} = (A_{knotsx,kx,x}) .' . z . A_{knotsy,ky, y}$

Следовательно, coefsb = coefs.

Сравнение и расширение

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

Если, например, определенные данные по 3D сетке содержатся в некотором 3D массиве v из размера [Nx,Ny,Nz], с v(i,j,k) содержа значение f (x (i), y (j), z (k)), затем вы начались бы с

coefs = reshape(v,Nx,Ny*Nz);

Принятие, что nj = knotsj - kj, для j = x,y,z, вы затем продолжили бы можно следующим образом:

sp = spap2(knotsx,kx,x,coefs.');
coefs = reshape(fnbrk(sp,'coefs'),Ny,Nz*nx);
sp = spap2(knotsy,ky,y,coefs.'); 
coefs = reshape(fnbrk(sp,'coefs'),Nz,nx*ny);
sp = spap2(knotsz,kz,z,coefs.'); 
coefs = reshape(fnbrk(sp,'coefs'),nx,ny*nz);

См. Главу 17 PGS или [К. де Бор, “Эффективная манипуляция с компьютерами продуктов тензора”, Математика Сделки ACM. Программное обеспечение 5 (1979), 173–182; Исправления, 525] для получения дополнительной информации. Те же ссылки также ясно дают понять, что нет ничего специального здесь об использовании приближения наименьших квадратов. Любой процесс приближения, включая интерполяцию сплайна, получившееся приближение которой имеет коэффициенты, которые зависят линейно от определенных данных, может быть расширен таким же образом к многомерному процессу приближения к данным с координатной сеткой.

Это точно, что используется в командах конструкции сплайна csapi, csape, spapi, spaps, и spap2, когда данные с координатной сеткой должны быть адаптированы. Это также используется в fnval, когда сплайн продукта тензора должен быть оценен на сетке.

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