张量积样条逼近- MATLAB & Simulink - MathWorks España万博1manbetx - 万博1manbetx,s manbetx 845,万博尤文图斯

张量积样条逼近

因为工具箱可以处理样条向量系数，很容易实现插值或逼近网格化数据张量积样条，如下图所示。你也可以运行例子“二元张量积样条”。

可以肯定的是，大多数张量积样条近似网格化数据可以直接获得样条构造命令之一，如spapi或csape在这个工具箱中，无需考虑本例中讨论的细节。相反，这个例子是为了说明张量积构造背后的理论，在这个工具箱中的构造命令没有涵盖的情况下，这将有所帮助。

本节讨论张量积样条问题的这些方面:

地点和结的选择

例如，考虑给定数据的最小二乘近似z（我，j) =f（x（我)，y（j))，我= 1:Nx，j= 1:纽约．从广泛使用的函数中获取数据Franke对表面拟合方案的测试(参见R. Franke，“分散数据插值的一些方法的关键比较”，海军研究生院技术学院。众议员nps - 53 - 79 - 0031979年3月)。定义域是单位正方形。中选择更多的数据站点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);

实际上，你同时发现离散最小二乘近似年代_{肯塔基州,knotsy}给每一个Nx数据集

$\begin{matrix} {（ y （ j ）， z （我， j ））}_{j ＝ 1}^{N y} ， & 我＝ 1 ： N x \end{matrix}$

特别是这些表述

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

提供数组瓦尔斯，其条目瓦尔斯(i, j)可以作为一个近似值吗f（x（我)，yy（j))的基本功能f在网格点x（我)，yy（j),因为瓦尔斯(:,j)的值是yy（j)的近似样条曲线sp．

命令获取的如下图所示:

网格(x, yy,瓦尔。”),视图(150年,50)

注意使用瓦尔斯。”，在网命令，因为MATLAB需要^®绘制数组时的面向矩阵视图。在二元近似中，这可能是一个严重的问题，因为习惯上是这样考虑的z（我，j)作为函数在点(x（我)，y（j))，而MATLAB认为z（我，j)作为函数在点(x（j)，y（我))。

一组假装成曲面的平滑曲线

请注意，每条平滑曲线的前两个和后两个值实际上都是零，因为前两个和后两个位置都在yy都在外面中样条的基本区间sp．

还要注意脊线。它们确认你只是在一个方向上绘制平滑曲线。

作为x函数的系数逼近

为了得到一个实际的表面，你现在必须更进一步。看看系数coefsy样条的sp：

Coefsy = fnbrk(sp，'coefs');

抽象地说，你可以把样条放在sp作为函数

$y | \to \sum_{r} c o e f 年代 y （：， r ） B_{r ， k y} （ y ）$

与我th条目coefsy (ir)矢量系数的coefsy (:, r)对应于x（我)，适用于所有人我．这建议近似每个系数向量coefsy (q,:)用相同顺序的样条kx用同样合适的绳结顺序knotsx．没有特别的原因，这次使用立方样条函数与四个均匀分布的内部结:

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

请注意,spap2（节,k, x,外汇)预计外汇（j:,)作为基准x (j)，即期望每一个列的外汇是一个函数值。拟合基准coefsy (q,:)在x（问)，适用于所有人问、现在spap2与转置的coefsy．

二元逼近

现在考虑系数的转置cxy的结果样条曲线：

Coefs = fnbrk(sp2，' Coefs ').';

它提供二元样条逼近

$（ x ， y ） | \to \sum_{问} \sum_{r} c o e f 年代（问， r ） B_{问， k x} （ x ） B_{r ， k y} （ y ）$

到原始数据

$（ x （我）， y （ j ）） | \to z （ x （我）， y （ j ）），我＝ 1 ： N x ， j ＝ 1 ： N y$

在网格(例如网格)上绘制此样条曲面

Xv = 0:.025:1;Yv = 0:.025:1;

您可以执行以下操作:

值=spcol (knotsx kx,十五)*系数* spcol (knotsy、肯塔基州、青年志愿)。”;网(十五、青年志愿值。”),视图(150年,50);

结果如下图所示。

Franke函数的样条近似

这很有道理，因为spcol (knotsx kx,十五)矩阵的(我，问)第一个条目等于值B_问，kx（十五（我))十五（我)的问第th条b样条kx对于结序列knotsx．

因为矩阵spcol (knotsx kx,十五)而且spcol (knotsy、肯塔基州、青年志愿)是带状，它可能更有效，尽管可能更多的内存消耗，为大十五而且青年志愿利用，利用fnval，详情如下:

Value2 =…fnval (spmak (knotsx fnval (spmak (knotsy系数),青年志愿)。“),十五)。”;

事实上，这是在内部发生的fnval是用张量积样条直接调用的，如在

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

以下是相对误差的计算，即给定数据与这些数据点的近似值与给定数据的大小之间的差值:

Errors = z - spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).';Disp (max(max(abs(错误)))/max(max(abs(z))))

输出为0.0539，也许不太令人印象深刻。然而，系数数组只是大小8 6

disp(大小(系数)

以适应大小的数据数组15日11．

disp(大小(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（十五（我)，y（j))的基本功能f在网格点处(十五（我)，y（j))。这在绘图时很明显valsb使用网：

网格(十五,y, valsb。”),视图(150年,50)

另一类平滑曲线假装成曲面

注意山脊。它们再次证明，你只是在一个方向上绘制平滑曲线。但这一次，曲线走向了另一个方向。

作为y函数的系数近似值

现在是第二步，得到实际的表面。首先，提取系数:

Coefsx = fnbrk(spb，'coefs');

然后拟合每个系数向量coefsx (r,:)用相同顺序的样条肯塔基州用同样合适的绳结顺序knotsy：

spb2 =spap2 (knotsy,肯塔基州,y, coefsx。');

注意，再一次，你需要转置系数数组spb,因为spap2将其最后一个输入参数的列作为数据值。

相应地，现在不需要转置系数数组coefsb的结果曲线：

Coefsb = fnbrk(spb2，'coefs');

二元逼近

结论是coefsb等于前面的系数数组系数，直到四舍五入，下面是测试:

Disp (max(max(abs(coefs - coefsb))))

输出为1.4433 e15汽油．

解释很简单:系数c样条的曲线年代包含在Sp = spap2(节，k,x,y)依赖线性在输入值上y．这意味着，考虑到两者c而且y是单行矩阵吗一个＝一个_节,k，x这

$c ＝ y {一个}_{结， k ， x}$

对于任何数据y．这句话甚至在什么时候也成立y是一个矩阵，大小d——- - - - - -N，在这种情况下，每个数据y(:,j)被认为是……中的一个点R^d，得到的样条为d-vector-value，因此它的系数数组c大小d——- - - - - -n， n =长度(节)- k．

特别是这些表述

Sp = spap2(knotsy,ky,y,z);coefsy = fnbrk (sp,“系数”);

给我们矩阵coefsy满足

$系数 y ＝ z ． {一个}_{knotsy,肯塔基州,y}$

后续计算

Sp2 = spap2(knotsx,kx,x,coefsy.');Coefs = fnbrk(sp2，' Coefs ').';

生成系数数组系数，考虑两次转置，则满足

$\begin{matrix} 系数＝（（ z {一个}_{knotsy,肯塔基州,y} ）＇． {一个}_{knotsx, kx, x} ）＇ \\ ＝（ {一个}_{knotsx, kx, x} ）＇． z ． {一个}_{knotsy,肯塔基州,y} \end{matrix}$

在第二种，替代，计算中，你首先计算

SPB = spap2(knotsx,kx,x,z.');Coefsx = fnbrk(spb，'coefs');

因此coefsx= z”。一个_{knotsx, kx, x}．后续计算

Spb2 = spap2(knotsy,ky,y,coefsx.');Coefsb = fnbrk(spb，'coefs');

然后提供

$coefsb ＝ coefsx ．＇． {一个}_{knotsy,肯塔基州,y} ＝（ {一个}_{knotsx, kx, x} ）．＇． z ． {一个}_{knotsy,肯塔基州, y}$

因此,coefsb＝系数．

比较与推广

第二种方法比第一种方法更对称，因为换位发生在每次调用中spap2没有其他地方。这种方法可以用于在任意数量的变量中近似网格数据。

例如，给定a上的数据三个-维网格包含在一些三维数组中v的大小(Nx、纽约、新西兰),v (i, j, k)包含值f（x（我)，y（j)，z（k))，然后开始

coefs =重塑(v,Nx,Ny*Nz);

假设nj =结j -kj，表示j =x, y, z，然后按以下步骤进行:

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

见第17章动力或(C。de Boor，“张量积的有效计算机操作”，s manbetx 845ACM反式。数学。软件5(1979), 173 - 182;更正，525]以获取更多细节。同样的参考文献也清楚地表明，这里使用最小二乘近似没有什么特别之处。任何近似过程，包括样条插值，其结果近似具有线性依赖于给定数据的系数，可以以同样的方式扩展到网格数据的多元近似过程。

这正是样条构造命令中所使用的csapi，csape，spapi，spaps,spap2，当拟合网格化数据时。它也用于fnval，当张量积样条要在网格上求值时。