专题三 二维数据插值拟合与最小二乘

易验证，上的双三次多项式可以唯一地写成
⎤ ⎡ ⎤ S ( x, y ) = ⎡ ⎣( x − xi −1 ) ⎤ ⎦⎡ ⎣ A ( hi ) ⎤ ⎦ [C ]ij ⎡ ⎣ A ( g j ) ⎦ ⎣( y − y j −1 ) ⎦ , ( x, y ) ∈ Rij , （3.5）
T T
其中符号 T 表示矩阵的转置， ⎡ ⎣ A ( hi ) ⎤ ⎦ 为 4 × 4 矩阵。
pi , j (i = 0,1, ⋅⋅⋅, n; j = 0,1, ⋅⋅⋅, m) 。分别对 j = 0,1, ⋅⋅⋅, m ，求解上 m 连续性方程就可求
出 R 的所有节点 ( xi , y j ) 处沿 x 方向的一阶偏导数 pij (i = 0,1, ⋅⋅⋅, n; j = 0,1, ⋅⋅⋅, m) 。 其次固定 x = xi ( i = 0,1, ⋅⋅⋅, n ) 则函数 S ( xi , y ) 也是一维三次样条函数。由边界 条件 qi 0和qim 及 m 连续性方程有
⎡ 2 / hi3 − 2 / hi3 1/ hi2 ⎢ 3 3 / hi3 − 2 / hi ⎢ −3 / hi A h = ⎡ ⎤ ( ) i ⎦ ⎣ ⎢0 0 1 ⎢ 0 0 ⎢ ⎣1
这是因为从可导出
1/ hi2 ⎤ ⎥ − 1/ hi ⎥ ⎥ 0 ⎥ 0 ⎥ ⎦
（3.6）
' '' si , j = S ( xi , y j ) , pij = S x ( xi , y j ) , qij = S y' ( xi , y j ) , ri, j = S xy ( xi , y j ) ,
λi ri−1,β + 2ri,β + μi ri+1,β = 3 ⎢λi
r −r ⎤ ⎡ ri,β − ri−1,β + μi i+1,β i,β ⎥ (i = 1,2, ⋅⋅⋅, n −1), β = 0, m. （3.9） hi hi+1 ⎦ ⎣
而边界条件 r0 β 和 rnβ ( β = 0, m) 已知，由此即可解出边界上的 ri 0 和
(3.1)
由此导出 R 上的一个矩形网格分割 Δ = Δx × Δy ，将 R 分成 mn 个子矩形
Rij : [ xi −1 , xi ] × [ yi −1 , yi ] 。
它的两条邻边长分别是
hi = xi − xi −1 ( i = 1, ⋅⋅⋅, n ) ， gi = yi − yi −1 ( j = 1, ⋅⋅⋅, m ) 。
i, பைடு நூலகம் i j ij
为插值双三次样条函数。 实际上，双三次样条函数是由两个一维三次样条函数作直积产生的。对任意 固定的 y0 ∈ [ c, d ] ， S ( x, y0 ) 是关于 x 的三次样条函数，同理，对任意固定的
x0 ∈ [ a, b ] ， S ( x0 , y ) 是关于 y 的三次样条函数。 n + 1 个节点的一维三次样条函数
(y− y ) ;
l j −1
（2） 在整 R 个上， 函数 S ( x, y ) 的偏导数 ∂α + β S ( x, y ) / ∂xα ⋅ ∂y β (α , β = 0,1, 2) 都 是连续的。则称 S ( x, y ) 为双三次样条函数。如果给定一数组
{s } ( i = 0,1, ⋅⋅⋅, n; j = 0,1, ⋅⋅⋅, m ) 则称为双三次样条。如果给定一数组，双三次样条 函数 S ( x, y ) 还满足插值条件 S ( x , y ) = s ( i = 0,1, ⋅⋅⋅, n; j = 0,1, ⋅⋅⋅, m ) ， 就称 S ( x, y )
专题三 二维数据的插值与拟合模型
如果我们要描述一个方桌或一座建筑物的形状， 我就要测量出其角点的位置 与高度，用直线连接角点就行了。但如果要我们描绘出一个不规则的曲面，如一 座山峰，一片海床面等等，问题就不会那么简单了，测量出足够密的样本点的高 度值，尤其是有特殊特征，如山脊、峡谷等位置的高度值，用直线连接相邻样本 点， 在图纸上或计算机上就能恢复曲线的形状， 但无疑这个工作量是相当巨大的。 在这一章里，我们讨论一个更一般的情形：在一个平面区域里，随机地测量一些 样本点的高度值， 如何描绘出原有的空间曲面？这类问题是典型的二维数据的插 值与拟合建模问题， 它有两个基本特点： 一是样本点是很有限的， 也就是并不 “足 够密” ，因此无法直接通过样本点恢复出曲线的形状，需要设计出合理的模型， 通过已知样本点，计算出其他未知点上的高度值；二是样本点是随机分布的，并 不能现成地套用一维的建模方法，需要克服如何从“随机样本点”到规则样本点 的过渡。同时也正是因为这两个特点，建模的方法具有很大的灵活性。 在这一讲里，我们仍然从现实问题入手，讨论 AMCM − 86 的 A 题的建模方 法，从中可以领会到二维数据的拟合这类问题的建模方法。内容是这样安排的， 第 1 节先给出一些必备的数学理论知识，第 2 节提出问题并作出合理的假设，同 时对问题进行分析，最后在第 3 节，我们给出一种解法。 §3.1 双三次样条函数 设 R:[ a, b ] × [ c, d ] 是 xy 平面上的一个矩形区域，在 x 轴和 y 轴上分别取定分 割 Δx : a = x0 < x1 < ⋅⋅⋅ < xn = b, Δy : c = y0 < y1 < ⋅⋅⋅ < yn = d .
λi pi −1, j + 2 pij + μi pi +1, j = 3 ⎢λi
⎡
sij − si −1, j hi
+ μi
si +1, j − si , j ⎤ ⎥ (i = 1, 2, ⋅⋅⋅, n − 1) hi +1 ⎦ （3.7）
在边界上，我们已知切向 p0 j 和 pnj ，由此即可解出
∗ ∗ λ∗ j qi , j −1 + 2qij + μ j qi , j +1 = 3 ⎢ λ j
⎡ ⎢ ⎣
sij − si , j −1 gj
+ μ∗ j
其中λ ∗ j =
gi +1 gi , μ∗ j = gi + gi +1 gi + gi +1
si , j +1 − si , j ⎤ ⎥ ( j = 1, 2, ⋅⋅⋅, m − 1) gj ⎥ ⎦ （3.8）
'' rαβ = S xy ( xα , yβ ) (α = 0, n; β = 0, m ) 。
（3.3）
给出了上面的条件后，我们有 定理 给 定 xy 平 面 上 的 矩 形 区 域 R ， 分 割 Δ 和 ( m + 3)( n + 3) 个 常 数
si , j , pα j , qiβ , rαβ ( i = 0,1, ⋅⋅⋅, m; α = 0, n; β = 0, m ) ， 则存在唯一的关于分割 Δ 的插值三
' '' 然记 S ( xi , y j ) = sij , S x ( xi , y j ) = pij , S y' ( xi , y j ) = qij , S xy ( xi , y j ) = rij ，则由 m 连续性方
程有
⎣ hi +1 hi 其中λi = , μi = hi + hi +1 hi + hi +1
直线 x = xi ( i = 0,1, ⋅⋅⋅, n ) 和直线 y = y j ( j = 0,1, ⋅⋅⋅, m ) 称为 Δ 分割的两族网格 线。网格线的交点 ( xi , y j ) ( i = 0,1, ⋅⋅⋅, n, j = 0,1, ⋅⋅⋅m ) 称为 Δ 分割的节点，总共有
( m + 1)( n + 1) 个节点，见图 3-1。
的全体构成 n + 3 维线性空间，所以要求求出插值三次样条函数除节点插值外， 还要加上两个边界条件。对双三次样条函数，关于分割 Δ 的双三次样条函数的全 体组成了 ( m + 3)( n + 3) 维线性空间，现在插值条件只有 ( m + 1)( n + 1) 个，还差
2 ( m + 1) + 2 ( n + 1) + 4 个边界条件。常用的边界条件是：
满足前面的条件。这样， S ( x, y ) 在子矩形 Rij 上的表示式（3.5）完全决定于矩阵
[C ]ij ，对 [C ]ij ，左上角 4 个元素 sij 由插值条件直接给出，剩下 12 个元素就要通
过边界条件解连续性方程求得。 对于一元三次样条插值，我们可导出 m 连续性方程和 M 连续性方程，并讨 论了其具体求解算法，在这里从略，只是简单的介绍以下思路。 首先，固定 y = y j ( j = 0,1, ⋅⋅⋅, m ) 则函数 S ( x, y j ) 是一维三次样条函数。如仍
即可解出 R 的所有节点 ( xi , y j ) 处沿 y 方向的一阶偏导数 qij (i = 0,1, ⋅⋅⋅, n; j = 0,1, ⋅⋅⋅, m) 。
' 最后，我们来求解 rij ，在 R 的两条边界 y = y0 和 y = ym 上，函数 S y ( x , y0 ) 和 ' 它们在各节点 xi 处的函数值 （实则为 S ( x, y0 ) Sy ( x, ym ) 是关于 x 的三次样条函数， ' 和 S ( x, ym ) 的关于 y 的一阶偏导数值） S y ( xi , y0 ) = qi 0 , S y' ( xi , ym ) = qim 已由（3.8） ' 式求出，同样可导出 S y ( x, y0 ) 和 S y' ( x, ym ) 关于 x 的 m 连续性方程
y
d
yi
Rij
yi −1
c
R
O
a
xi −1
xi
xi +1
x
图 3-1 矩形分割 Δ
设函数 S ( x, y ) 定义在矩形域 R 上且满足下列条件： （1）在每个子矩形 Rij ( i = 1, ⋅⋅⋅, n, j = 1, ⋅⋅⋅m ) 上， S ( x, y ) 是一个关于 x 和 y 都
ij 是三次的多项式函数，即 S ( x, y ) = ∑∑ bkl ( x − xi −1 ) k =0 l =0 3 3 k
。 次样条函数 S ( x, y ) 满足插值条件及边界条件（3.2）和（3.3） 这个定理的证明并不难，但属纯数学范畴，在这我们不再引述。我们关心的 是如何求出插值双三次样条函数的具体表达式。 由定义， S ( x, y ) 限制在单个子矩形 Rij 上时为双三次多项式。记 Rij 四个顶点 上的函数值为 sij , si −1, j , si , j −1 , si −1, j −1 ，两个方向的一阶偏导数为 pi −1, j −1 , pi −1, j , pi , j −1 , pij
（1）矩形 R 的四条边界条件上所有节点处的一阶法向偏导数 ' ⎧ ⎪ pα j = S x ( xα , y j ) ( j = 0,1, ⋅⋅⋅, m; α = 0, n ) （3.2） ⎨ ' = = ⋅⋅⋅ = q S x , y i 0,1, , n ; β 0, m ( ) ( ) ⎪ i β y i β ⎩ （2）矩形 R 的四个角点处的二阶混合偏导数
[C ]ij
（3.4）
又记
3 2 ⎡ ⎤ ⎡ ⎣( x − xi −1 ) ⎤ ⎦ = ⎣( x − xi −1 ) , ( x − xi −1 ) , ( x − xi −1 ) ,1⎦ 3 2 ⎡( y − y j −1 ) ⎤ = ⎡( y − y j −1 ) , ( y − y j −1 ) , ( y − y j −1 ) ,1⎤ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦
和 qi −1, j −1 , qi −1, j , qi , j −1 , qij 以及二阶混合偏导数为 ri −1, j −1 , ri −1, j , ri , j −1 , rij ，也就是给定了四
阶方阵
⎡ si −1, j −1 si −1, j qi −1, j −1 qi −1, j ⎤ ⎢ ⎥ si , j qi , j −1 qi , j ⎥ ⎢ si , j −1 =⎢ p pi −1, j ri −1, j −1 ri −1, j ⎥ ⎢ i −1, j −1 ⎥ ⎢ ⎥ p p r r i, j i , j −1 i, j ⎦ ⎣ i , j −1