介绍

由于曲线拟合工具箱可以处理样条与向量利用张量积样条对网格数据进行插值或逼近是很容易的。工具箱中的大多数样条构造命令都利用了这一点。

但是，您可能有兴趣了解如何对二元数据使用张量积近似网格数据的详细描述。当您需要一些工具箱中的命令未提供的张量积构造时，这也会很有用。下载188bet金宝搏

示例:网格数据的最小二乘近似

例如，考虑给定数据的最小二乘近似。

z(i,j) = f(x(i)，y(j)) for i = 1: i,j = 1: j。

这是一些网格数据，取自Franke的样本函数。请注意，网格在边界附近有些密集，这有助于确定那里的近似。

x=sort（[（0:10）/10.03.07.93.97]）；y=sort（[（0:6）/6.03.07.93.97]）；[xx，yy]=ndgrid（x，y）；%注意:ndgrid而不是meshgridz =因特网(xx和yy);网格(x, y, z。');包含(“x”）;ylabel (“y”）;视图(150年,50);标题(来自Franke函数的数据）;

关于NDGRID和MESHGRID的注意事项

注意这些语句

[xx, yy] = ndgrid (x, y);

z =因特网(xx和yy);

用在上面确定z（i，j）是在网格点处近似的函数值(x(我),y (j))．

不过，MATLAB®命令网格(x, y, z)期望z（j，i）(注意颠倒的顺序我和j)作为网格点上的值(x(我),y (j))．因此，上面的图是由语句生成的

网格(x, y, z。');

i、 例如，使用矩阵的转置z．

如果我们使用的话，这样的转换就没有必要了meshgrid而不是ndgrid．但由此产生的z不会遵循近似理论的标准。

样条空间在y方向的选择

接下来，我们选择样条线顺序肯塔基州还有一个结序列knotsy在y方向上

肯塔基州= 3;knotsy = augknt([0二十五分。5。1),肯塔基州);

然后获得

sp = spap2 (knotsy,肯塔基州,y, z);

样条曲线，其我-th分量是曲线的近似(y, z(我:))为我= 1:我．

评价

特别是,

yy =约:.05:1.1;瓦尔斯= fnval (sp、yy);

创建一个矩阵瓦尔斯谁的（i，j）-th元素可以作为该值的近似f (x(我),yy (j))基本功能的定义f在格点（x（i），yy（j））．这在我们绘图时很明显瓦尔斯．

网格(x, yy, vals。');包含(“x”）;ylabel (“y”）;视图(150年,50);标题(“同时逼近y方向上的所有曲线”）;

请注意，对于每个十（一），前两个和后两个值都是零，因为前两个和后两个站点都在yy是否在样条的基本区间之外sp．

还要注意沿着y方向的“脊”，最明显的是在表面的峰值附近。它们证实了我们只在一个方向上绘制平滑的曲线。

从曲线到曲面;在x方向上选择一个样条空间

为了得到一个实际的表面，我们现在必须进一步考虑系数。coefsy样条的sp，由

系数=fnbrk（sp，“c”）;

抽象地说，你可以考虑样条sp作为向量值函数

sum coefsy(:，r) B_{r,ky}(y)

r

与我th元素,系数（i，r），向量系数的coefsy (:, r)对应于十（一）为我= 1:我。这建议近似每个曲线（x，coefsy（：，r））通过样条，用同样的顺序kx和相同的结序列knotsx对于每一个r．

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

使用spap2这里的指挥官也许需要一个解释。

回想一下,spap2（节、k、x、fx）对待外汇(:,j)作为x（j），即每个列的外汇作为数据值。因为我们想要符合这个值coefsy(我,:)在十（一），全部我，我们必须提供spap2与转置的coefsy．

现在考虑转置所得到的样条“曲线”的系数矩阵sp2,获得

coefs=fnbrk（sp2，“c”)”。;

系数提供二元样条逼近

(x, y) | - >和系数之和(q, r) B_ {q, kx} (x) B_ {r,肯塔基州}(y)

q r

恢复原始数据

f(x(i)，y(j)) = z(i,j)

我们使用spcol提供值B_ {q, kx}(十五(i))和B{r，ky}（yv（j））需要计算样条曲面在某些网格点上的值（xv（i）、yv（j））然后画出这些值。

xv=0.025:1；yv=0.025:1；值=spcol（knotsx，kx，xv）*coefs*spcol（knotsy，ky，yv）。'；网格（xv，yv，值）。；xlabel(“x”）;ylabel (“y”）;视图(150年,50);标题(“样条近似式”）;

为什么这种评估是有效的?

该声明

= spcol(knotsx,kx,xv) * coefs * spcol(knotsy,ky,yv).'

上面的用法很有意义，因为，例如，spcol (knotsx kx,十五)矩阵是谁的(问我)-th条目等于该值B_ {q, kx}(十五(i))在十五(我)的问b样条的顺序kx对于结序列knotsx，当我们想对表达式求值时

f (q,r) B_{q,kx}(x) B_{r,ky}(y)

q r

= sum B_{q,kx}(x) coefs(q,r) B_{r,ky}(y)

q r

在(x, y) =(十五(我),青年志愿(j))．

更有效的选择

由于矩阵spcol (knotsx kx,十五)和spcol (knotsy、肯塔基州、青年志愿)如果是带状的，则对于“大型”可能更有效十五和伊夫（尽管可能更消耗内存）来利用fnval．

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

事实上,fnval和spmak可以直接处理多元样条，因此上面的表述可以用

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

更好的是，近似的构造可以用一个召唤spap2，因此我们可以通过语句直接从给定的数据中获得这些值

value4=fnval（spap2（{knotsx，knotsy}，[kx-ky]，{x，y}，z），{xv，yv}）；

检查一下

这里是一个检查，具体来说，是相对以这四种不同方式计算的值之间的差异。

差值=绝对值（值-值2）+绝对值（值-值3）+绝对值（值-值4）；最大值（最大值（差值））/最大值（最大值（绝对值）））

ans = 1.1206 e15汽油

这四个方法返回相同的值，直到舍入错误。

近似误差

这是一个误差图，即给定数据值与这些数据点的样条逼近值之间的差值。

Errors = z - spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).';网格(x, y,错误。');包含(“x”）;ylabel (“y”）;视图(150年,50);标题(“给定数据站点的错误”）;

的相对错误是

最大值（最大值（绝对值（误差））/最大值（绝对值（z）））

ans = 0.0539

这可能不是太令人印象深刻。另一方面，比率

(使用自由度)/(数据点数)

只是

努梅尔（coefs）/努梅尔（z）

ans=0.2909

这种方法的明显偏差

这里采用的方法似乎有偏见的:我们首先考虑给定的数据值z的向量值函数y然后，我们将由近似曲线的向量系数构成的矩阵视为描述的向量值函数x．

当我们以相反的顺序来处理事情，也就是，思考z的向量值函数x，然后将由近似曲线的向量系数构成的矩阵视为描述y？

也许令人惊讶的是，最终的近似值是相同的，直到舍入。下一节包含数值实验来证实这一点。

以另一种方式执行：从X中的样条线空间开始

首先，我们用样条曲线拟合数据，但这次用x作为自变量，因此它是行的z现在变成了数据值。相应地，我们必须供给z、 "(而不是z)spap2，并获得

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

所有曲线的样条近似（x，z（：，j））为j = 1: j.特别是，

十五,valsb = fnval (spb)。';

创建一个矩阵（i，j）-th元素可以作为该值的近似f（十五（i）、y（j））基本功能的定义f在格点(十五(i), y (j))．这在我们绘图时很明显瓦尔斯布．

网格(十五,y, valsb。');包含(“x”）;ylabel (“y”）;视图(150年,50);标题(“同时逼近x方向上的所有曲线”）;

再次注意这些脊线，这一次是沿着x方向运行的。它们再次确认，我们只在一个方向绘制平滑曲线。

从曲线到曲面:在y方向上使用样条空间

现在是第二步，得到实际的表面。

让coefsx是spb,也就是说,

coefsx=fnbrk（spb，“c”）;

抽象地说，你可以考虑样条spb作为向量值函数

B_{r,kx}(x)

r

与jth条目coefsx（r，j）向量系数coefsx (r,:)对应于y (j)，全部j.因此，我们现在拟合每条曲线（y，coefsx（r，：））通过样条，用同样的顺序肯塔基州和相同的结序列knotsy每人r．

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

在建设spb2，我们再次需要转置系数矩阵spb自从spap2将其最后一个输入参数的列作为数据值。

因此，现在不需要转置系数矩阵coefsb由此产生的“曲线”。

coefsb=fnbrk（spb2，“c”）;

索赔：coefsb等于前面的系数数组系数，直至四舍五入。有关这方面的证明，请参阅曲线拟合工具箱文档中有关张量积构造的讨论。在这里，我们只需进行以下检查。

最大值（最大值（绝对值（系数-系数）））

ans=1.9984e-15

因此,二元样条逼近

(x, y) | - >和总和coefsb (q, r) B_ {q, kx} (x) B_ {r,肯塔基州}(y)

q r

恢复原始数据

f(x(i)，y(j)) = z(i,j)

获得的与生成的前一个相同系数而不是coefsb．

正如前面所观察到的，您可以使用两个语句（一个用于构造最小二乘逼近，另一个用于在矩形网格上求值）执行我们刚刚经历的整个构造（以两种方式）。

tsp = spap2 ({knotsx, knotsy}, kx,肯塔基州,{x, y}, z);valuet = fnval (tsp,{十五,青年志愿});

这里，作为另一个检查，是早先计算的值和现在计算的值之间的相对差:

马克斯(max (abs (values-valuet))) / max (max (abs(值)))

ans = 3.7353 e-16

另一个例子:插值

由于数据来自平滑函数，我们应该对其进行插值，即使用spapi而不是spap2，或同等地使用spap2使用适当的结序列。为便于说明，以下是使用spapi．

回顾一下，(单变量)数据点为

x

x =1×150 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×110 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),肯塔基州);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,十五)* icoefs * spcol (knotsy、肯塔基州、青年志愿)。”;

然后画出插值图。

网格(十五,青年志愿,ivalues。');包含(“x”）;ylabel (“y”）;视图(150年,50);标题(“样条Interpolant”）;

同样，上面的步骤可以用两个语句来完成，一个用于插值的构造，另一个用于在矩形网格上求值。

tsp = spapi ({knotsx, knotsy}, {x, y}, z);valuet = fnval (tsp,{十五,青年志愿});

为了便于检查，我们还计算了先前计算的值与现在计算的值之间的相对差值。

马克斯(max (abs (ivalues-valuet))) / max (max (abs (ivalues)))

ans = 3.6712 e-16

近似误差

接下来，我们计算插值的误差作为Franke函数的近似。

fvalues =因特网(repmat(十五。',1,长度(青年志愿),repmat(青年志愿、长度(十五),1));Error = fvalues - ivalues;网格(十五、青年志愿错误。');包含(“x”）;ylabel (“y”）;视图(150年,50);标题(“内插误差”）;

的相对近似误差是

max（max（abs（误差））/max（max（abs（fvalue）））

ans = 0.0409