第 五 章 插 值 法 与曲线拟合 插值法

上一讲的简单回顾 ● 插值多项式的存在惟一性: 满足插值条件 Pn(xi)=f(xi), ( i=0,1,2,…,n) n次插值多项式Pn(x)=a0+a1x+a2x2+……+anxn 存在而且惟一. ● 插值余项: Rn(x)= f(x)- Pn(x)= , ● Lagrange插值多项式 ,i =0,1,…,n 其中 称为Lagrange插值基函数

§1.3 牛顿途径 对于n+1个不同的节点 ，考虑n次多项式 （6） 如果满足： ，那么它就是n+1个点上的n次插值多项式，对于这样的 ，有

本讲主要内容: 缺点: 不具有承袭性,即每当增加一个节点时,不仅要增加求 和的项数,而且以前的各项也必须重新计算. 优点: 具有严格的规律性,便于记忆. 缺点: 不具有承袭性,即每当增加一个节点时,不仅要增加求 和的项数,而且以前的各项也必须重新计算. 为了克服这一缺点,本讲将建立具有承袭性的插值公式— Newton插值公式. 本讲主要内容: ● Newton插值多项式的构造 ● 差商的定义及性质 ● 差分的定义及性质 ● 等距节点Newton插值公式

一 基函数 问题1 求作n次多项式 使满足 (1) (2) 为了使 的形式得到简化,引入如下记号 (3)

称为Newton插值的以x0,,x1,…,xn为节点的基函数,即 可以证明,这样选取的基函数是线性无关的,由此得 出的问题4.1的解便于求值,而且新增加一个节点 xn+1时 只需加一个新项 即 而 (4)

依据条件(2),可以依次确定系数c0,c1,…,cn..例如, 取x=x0,,得 取x=x1 ,得

. 为了得到计算系数ci的一般方法,下面引进差商的概念.

二 差商的定义 给定[a,b]中互不相同的点x0,x1,x2,…,以及f(x)在这些点处相 二 差商的定义 给定[a,b]中互不相同的点x0,x1,x2,…,以及f(x)在这些点处相 应的函数值f(x0),f(x1),f(x2),…,用记号 表示f(x)在x0及x1两点的一阶差商. 用记号 表示f(x)在x0,x1,x2三点的二阶差商. 一般地,有了k-1阶差商之后, 可以定义f(x)在x0,x1,..,xk的k阶差商

三 Newton插值公式 由差商定义,有 f(x)= f[x0]+(x-x0)f[x,x0] f[x,x0]= f[x0,x1]+(x-x1)f[x,x0,x1] f[x,x0,x1]= f[x0,x1,x2]+(x-x2)f[x,x0,x1,x2] ……….. f[x,x0,…xn-1]= f[x0,…,xn]+(x-xn)f[x,x0,….,xn] 将以上各式，由下而上逐步代入，得到 　f(x)= f(x0)+(x-x0) f[x0,x1]+(x-x0)(x-x1) f[x0,x1,x2] +…+(x-x0)…(x-xn-1) f[x0,…,xn] (5) +(x-x0)…(x-xn-1)(x-xn)f[x,x0,…xn]

记 Nn(x)= f(x0)+(x-x0) f[x0,x1]+(x-x0)(x-x1) f[x0,x1,x2] +…+(x-x0)…(x-xn-1) f[x0,…,xn] (6) Rn(x)= (x-x0)…(x-xn) f[x,x0,…,xn]= f[x,x0,…,xn]wn+1(x) (7) 则(5)可表示为 f(x)= Nn(x)+ Rn(x) (8) 显然, Nn(x)是次数不超过n的多项式,且有 Rn(xi)= f[x,x0,…,xn]wn+1(xi)=0, i=0,1,…,n 即 Nn(xi)= f(xi) , i=0,1,…,n 由此可知,如此构造的函数Nn(x)就是问题1的解,且 ci =f[x0,…,xi] , i=0,1,…,n (9)

公式(6)称为函数f(x)在节点x0,…,xn上的n阶Newton 插值公式, (7)式称为Newton插值公式余项,即截断误差.注 意到,余项表达式(7)只要求被插值函数f(x)在插值区间[a,b] 上连续. 由函数f(x)的插值多项式的惟一性,函数f(x) 的Newton 插值多项式与Lagrange插值多项式实为同一个多项式,即 Nn(x) ≡ Ln(x) 两者不过是表现形式不同而已.由此有: 若f(x)∈Cn+1[a,b], 则有 Rn(x)= f[x,x0,…,xn]wn+1(x)= (10) ,

四 差商的性质 性质1(差商与函数值的关系) 证明： 性质2(对称性) 其中i0,…,ik是0,1,…k的任排列. 证明: 由性质1可知.

性质3(差商与导数的关系) f(x)∈Ck[a,b], 证明: 由式(10)即得. 性质4(多项式的差商)设f(x)为n次多项式,则其一阶差商 是x的n-1次多项式 证明: 推论 n次多项式pn(x)的k阶差商pn[x0,…xk],当k≤n时 是n-k次多项式,当k>n时恒为0

运用Newton插值公式(6)进行计算时,先计算f(x)关于 节点x0,…,xn的各阶差商.计算过程如下表所示 . xk f(xk) 一阶差商 二阶差商 三阶差商 n 阶差商

计算Nn(x)时,常采用秦九韶算法,即 . 下面给出Newton插值法的计算机算法.开始时,f(k)存放函数值f(xk),运算完毕f(k)存放k阶差商f[x0,…,xk]

Newton插值算法 (1) 输入: xi, fi; di=fi (i=0,1,…,n); 计算差商 对i=1,2,…,n 做 (3.1) 对j=i,i+1,…,n 做 fj=(dj-dj-1)/(xj-xj-i); (3.2) 对j=i,i+1,…,n 做 dj=fj; (4) 计算插值N(u) (4.1)输入插值点u; (4.2) v=0; (4.3) 对i=n,n-1,…,1,0 做 v=v(u-xi)+fi; (5) 输出u,v.

五 等距节点的Newton插值公式与差分 (6)可以简化. 当插值节点x0,…,xn为等距分布时,Newton插值公式 (6)可以简化. 设插值节点xj=x0+jh, j=0,1,…,n; h=(b-a)/n称为步长,且 x0=a, xn=b. 令x=x0+th,则当x0≤x≤xn时,0≤t≤n. 基函数 此时差商也可进一步化简,为此引进差分的概念. 定义 称 △f(xi)=f(xi+h)-f(xi) 和 ▽f(xi)= f(xi)-f(xi-h) 分别为函数f(x)在点xi处的一阶向前差分和一阶向后差分.

一般地,称k阶差分的差分为k+1阶差分,如二阶向前和向 后差分分别为 计算各阶差分可按如下差分表进行.

向前差分表

性质3(多项式的差分)若f(x)∈Pn(n次多项式类), 则 差分具有如下性质: . 性质1(差分与函数值的关系) 各阶差分均可表示为函值 fj=f(xj)的线性组合: 性质2(前差与后差的关系): 性质3(多项式的差分)若f(x)∈Pn(n次多项式类), 则 其中

性质4(差分与差商的关系): 性质5(差分与导数的关系 利用这些性质,可将Newton公式 Nn(x)= f(x0)+(x-x0) f[x0,x1]+(x-x0)(x-x1) f[x0,x1,x2] +…+(x-x0)…(x-xn-1) f[x0,…,xn] 进一步简化

(11) 称公式(11)为Newton向前差分插值公式,其余项为 (12) 如果将式(6)改为按节点xn,xn-1,…,x0的次序排列的Newton插 值公式,即

令x=xn-th, 则当x0≤x≤xn时,0≤t≤n.利用差商与向后 差分的关系, 式(13)可简化为 (13) (14) 称式(14)为Newton向后差分插值公式,其余项为

若要计算的插值点x较靠近点x0,则用向前插值公式(4. 8),这时t=(x-x0)/n的值较小,数值稳定性较好 若要计算的插值点x较靠近点x0,则用向前插值公式(4.8),这时t=(x-x0)/n的值较小,数值稳定性较好.反之,若x靠近xn,,,则用向后插值公式(14). 利用向前与向后差分的关系(差分性质2): 式(14)可表示成 (15) 这样,计算靠近x0或xn的点的值时,都只需构造向前差分表

例 给定f(x)在等距节点上的函数值表如下: xi 0.4 0.6 0.8 1.0 f(xi) 1.5 1.8 2.2 2.8 分别用Newton向前和向后差分公式,求f(0.5)及f(0.9)的近似值. 解 先构造向前差分表如下: xi fi △fi △2fi △3fi 0.4 1.5 0.6 1.8 0.3 0.8 2.2 0.4 0.1 1.0 2.8 0.6 0.2 0.1 x0=0.4, h=0.2, x3=1.0. 由(4.8)和(4.12),分别用差分表中对角线上的值和最后一行的值,得Newton向前和向后插值公式如 下:

(1) (2) 当 x=0.5时,用公式(1),这时t=(x-x0)/h=0.5. 将t=0.5代入(1),得 f (0.5)≈N3(0.5)=1.64375. 当x=0.9时, 用公式(2), 这时t=(x3-x)/h=0.5. 将t=0.5代入(2), 得 f (0.9)≈N3(0.9)=2.46875.

课后练习题 P217: 第2题, 第4题, P219: 第11题.

§6 曲线拟合 6.1.2 曲线拟合问题 仍然是已知 x1 … xm ; y1 … ym, 求一个简单易算的近似函数 f(x) 来拟合这些数据。 但是① m 很大； ② yi 本身是测量值，不准确，即 yi  f (xi) 这时没必要取 f(xi) = yi , 而要使 i=f(xi)  yi 总体上 尽可能地小。 这种构造近似函数 的方法称为曲线拟合，f(x) 称为拟合函数 称为“残差”

“使 i=P(xi)  yi 尽可能地小”有不同的准则 常见做法： 使 最小 不可导，求解困难，P283 使 最小 使 最小

6.2 线性拟合问题 6.2.1 ||.||2 意义下的线性拟合（线性最小二乘问题） 确定拟合函数 ，对于一组数据(xi, yi) (i = 1, 2, …, m) 使得 达到极小，这里 n <＝ m。 Denote:

称方程组Ax=b为超定方程组

记 E 实际上是 c0, c1, …, cn 的多元函数，在 E 的极值点应有

得到关于c1,c2,…,cn的方程组 法方程组(或正规方程组)

例1 数据 ti 0 20 40 60 80 100 fi 81.4 77.7 74.2 72.4 70.3 68.8

6.3 线性最小二乘问题 设A是m×n阶矩阵(m>n), 称线性方程组 Ax=b (1) 为超定方程组; 这里x∈Rn,b∈Rm. 6.3 线性最小二乘问题 设A是m×n阶矩阵(m>n), 称线性方程组 Ax=b (1) 为超定方程组; 这里x∈Rn,b∈Rm. 如果A的秩r(A)=n, 称A为列满秩矩阵. 记残向量r=b-Ax，考虑确定一个向量x， 使‖r‖2 2＝‖b-Ax‖2 2, 达到最小的问题称为线性最小二乘问题， 这样的x称为方程组(1)的最小二乘解.

6.3.4 最小二乘解的存在惟一性 结论1 ：设A是m×n阶矩阵，x∈Rn, b∈Rm. 由线性方程组理论可知，线性方程组 6.3.4 最小二乘解的存在惟一性 结论1 ：设A是m×n阶矩阵，x∈Rn, b∈Rm. 由线性方程组理论可知，线性方程组 Ax=b (24) 有解的充分必要条件是  r (A)= r (A|b). (25)

定理6. 3. 7假设方程组(24)有解，令x是其一个解. 那么，方程组(24)的所有解的集合为{x}+N(A) 定理6.3.7假设方程组(24)有解，令x是其一个解. 那么，方程组(24)的所有解的集合为{x}+N(A). 方程组(24)有 惟一解的充分必要条件是null(A)=0. 这里， null (A)表示A的核子空间的维数.

证明: 首先证明任意的向量y∈{x}+N(A)都是方程组(24)的解. 事实上，将y记为y=x+z, 其中z∈N(A)， 即Az=0，x∈{x}. 因此， Ay=Ax+Az=b,即y满足方程组(24). 反过来， 若y满足方程组(24)， 有 Ay-Ax =A(y-x)= 0， 即y-x∈N(A). 记y=x+(y-x)，从而有y∈｛x｝+N(A). 惟一性. 因为齐次方程组Ax=0有惟一零解的 充 分必要条件是A为满秩矩阵，即null (A)=0.

定理6.3.8当m>n时， 超定方程组(1)的最小二乘解总是存在的. 最小二乘解惟一的充分必要条件是null (A)=0.

证: 记b=b１＋b２， 其中b1∈R(A)，b2∈N(AＴ）. 对任意x∈Rn, Ax∈R(A), b1-Ax∈R(A). 因此， ‖r‖22＝‖b-Ax‖22＝‖(b１-Ax)+b２‖22. 由定理6.3.3的推论1和定理6.3.2， ‖r‖22＝‖b１-Ax‖22+‖b２‖22． 要使‖r‖22达到最小等价于确定x，使‖b１-Ax‖22 为0， 即求方程组Ax=b１的解x. 因为b１,Ax, b１－Ax都是R(A)中的向量，因此，可以 把b１看成由A的列向量线性表示， 即b１＝Ax. 换句话说，方程组Ax=b１的解总是存在的，从而方程 组(1)的最小二乘解也总是存在的. 惟一性的证明可直接由定理6.3.7得到.

6.3.1 正交性的有关性质 在线性代数欧氏空间理论中， 将R3中两个向量x,y之间的夹角φ满足的关系式 6.3.1 正交性的有关性质 在线性代数欧氏空间理论中， 将R3中两个向量x,y之间的夹角φ满足的关系式 xTy=‖x‖2‖y‖2cosφ (2) 推广到Rn. 设x,y∈Rn， 由Cauchy不等式  -1≤ ≤1 从而得到Rn中两个向量之间的夹角为 φ=arccos (3)

定理6.3.1 设x, y是Rn中的向量， x与y正交 的充分必要条件为xＴy=0.

定理6.3.2：设x， y∈Rn， 且x⊥y，那 么： ‖x+y‖22＝‖x‖22＋‖y‖22． = xＴx+2yＴx+yＴy 而xＴy=yＴx=0, 因此 ‖x+y‖22＝‖x‖22+‖y‖22 注：推广到Rn中的向量组α１，α２，…，αk， 如果αiＴαj=0 (i≠j), 称α1，α2，…，αk是 正交向量组. 特别地: 如果‖αi‖2=1(i=1,2,…，k)， 即 αiTαj=δij，称α1，α2，…，αk 为标准正交向量组.

设U是Rn中的子空间， x∈Rn. 如果x与U中任意向量正交， 称向量x与子空间U正交， 记为x⊥U. 设U，V是Rn中两个子空间， 如果任意x∈U和任意y∈V是正交的， 称子空间U与子空间V正交， 记为U⊥V. 设U，V是Rn中互补的子空间. 如果U⊥V， 那么称U，V互为正交补子空间， 记U=V⊥或V=U⊥. 可以证明， 一个子空间的正交补子空间是惟一的.

定理6.3.3 设A是n×k阶矩阵，x∈Rn, 那么下列三种情况是 等价的： ①x⊥R(A); ②AＴx=0; ③x∈N(AＴ). 这里，N(AＴ)={AＴx=0, x∈Rn}称为AＴ的核子空间. 证：由N(AＴ)的定义， ②与③显然等价. 下面证明①与②等价. 记A=(α1,α2,…，αk)， 那么，αi∈R(A) (i=1,2,…,k). 假设x⊥R(A), 即αiＴx=0 (i=1,2,…，k). 从而AＴx=0. 另一方面，如果AＴx=0, 那么有z∈Rk， 使Az=y∈R(A). 这时，yＴx=zＴAＴx=0,即x⊥y. 由z的任意性， 得Az是任意的， 因此x⊥R(A). 由这个定理， 容易得到： 推论1设A是n×k阶矩阵， 那么R(A)有惟一的正交补子 空间N(AＴ).

6.3 线性最小二乘问题 设A是m×n阶矩阵(m>n), 称线性方程组 Ax=b (1) 6.3 线性最小二乘问题 设A是m×n阶矩阵(m>n), 称线性方程组 Ax=b (1) 为超定方程组; 这里x∈Rn, b∈Rm. 如果A的秩r (A) =n, 称A为列满秩矩阵. 记残向量r=b-Ax， 考虑确定一个向量x， 使‖r‖2 2＝‖b-Ax‖2 2, 达到最小的问题称为线性最小二乘问题， 这样的x称为方程组(1)的最小二乘解.

6.3.1 正交性的有关性质 在线性代数欧氏空间理论中， 将R3中两个向量x,y之间的夹角φ满足的关系式 6.3.1 正交性的有关性质 在线性代数欧氏空间理论中， 将R3中两个向量x,y之间的夹角φ满足的关系式 xTy=‖x‖2‖y‖2cosφ (2) 推广到Rn. 设x,y∈Rn， 由Cauchy不等式  -1≤ ≤1 从而得到Rn中两个向量之间的夹角为 φ=arccos (3)

定理6.3.1 设x, y是Rn中的向量， x与y正交 的充分必要条件为xＴy=0.

定理6.3.2：设x， y∈Rn， 且x⊥y，那 么： ‖x+y‖22＝‖x‖22＋‖y‖22． = xＴx+2yＴx+yＴy 而xＴy=yＴx=0, 因此 ‖x+y‖22＝‖x‖22+‖y‖22 注：推广到Rn中的向量组α１，α２，…，αk， 如果αiＴαj=0 (i≠j), 称α1，α2，…，αk是 正交向量组. 特别地: 如果‖αi‖2=1(i=1,2,…，k)， 即 αiTαj=δij，称α1，α2，…，αk 为标准正交向量组.

设U是Rn中的子空间， x∈Rn. 如果x与U中任意向量正交， 称向量x与子空间U正交， 记为x⊥U. 设U，V是Rn中两个子空间， 如果任意x∈U和任意y∈V是正交的， 称子空间U与子空间V正交， 记为U⊥V. 设U，V是Rn中互补的子空间. 如果U⊥V， 那么称U，V互为正交补子空间， 记U=V⊥或V=U⊥. 可以证明， 一个子空间的正交补子空间是惟一的.

定理6.3.3 设A是n×k阶矩阵，x∈Rn, 那么下列三种情况是 等价的： ①x⊥R(A); ②AＴx=0; ③x∈N(AＴ). 这里，N(AＴ)={AＴx=0, x∈Rn}称为AＴ的核子空间. 证：由N(AＴ)的定义， ②与③显然等价. 下面证明①与②等价. 记A=(α1,α2,…，αk)， 那么，αi∈R(A) (i=1,2,…,k). 假设x⊥R(A), 即αiＴx=0 (i=1,2,…，k). 从而AＴx=0. 另一方面，如果AＴx=0, 那么有z∈Rk， 使Az=y∈R(A). 这时，yＴx=zＴAＴx=0,即x⊥y. 由z的任意性， 得Az是任意的， 因此x⊥R(A). 由这个定理， 容易得到： 推论1设A是n×k阶矩阵， 那么R(A)有惟一的正交补子 空间N(AＴ).

6.3.2 矩阵的QR分解 定理6.3.4 设A=(α1,α2,…，αn)是列满秩矩阵，αi∈Rm (i=1,2,…，n)且m≥n. 那么， A有惟一的QR分解，记为 A=QR， (4) 这里，Q是有n个标准正交列的m×n阶矩阵，R是有正对角元的n阶上三角矩阵. 证：由A是列满秩矩阵可知， AＴA是n阶正定矩阵， 因此有 惟一的Cholesky分解： AＴA=RＴR, (5) 这里R是有正对角元的上三角矩阵， R－１存在. 令Q=AR－１, (6) 那么， QＴQ=R－ＴAＴAR－１＝In, 即Q是有标准正交列的m×n阶矩阵. 由(6)式， (4)式成立， 且由(5)式的惟一性， 分解式(4)也是惟一的.

6.3.3 Householder矩阵与矩阵的正交三角化 定义1设w是欧氏空间Rn中的单位向量， 形如H=I-2wwＴ (10) 的n阶矩阵称为Householder矩阵， 也称 为反射(镜像)矩阵或称为Householder变 换(反射变换、镜像变换).

定理6.3.5 设H是Rn中的Householder 矩阵， 那么， ①HＴ=H (对称性)； ②H－１=HＴ (正交性)； ③H2=I (对合性). 证：直接验证， 性质①成立. 由H=I-2wwＴ, wＴw=1和性质①， H2＝HＴH=(I-2wwＴ)Ｔ (I-2wwＴ) =I-4wwＴwwＴ+4wwＴ=I, 即HＴ=H-1，性质②和③成立.

定理6.3.6设Rn中有非零向量x≠y 且‖x‖2=‖y‖2，那么，存在 Householder矩阵H，使Hx=y 证：不妨令w=(x-y)/ ‖x-y‖2 ,有wTw=1.设H=I- 2wwT ,那么 Hx=(I-2wwT) 由已知， xＴx=yＴy. 所以 (x-y)Ｔ(x-y)=xＴx- 2yＴx+yＴy=2(xＴx-yＴx). 从而， Hx=x-(x-y)=y. 证明中如果令u=x-y， 那么Householder矩阵形为(11)式，即H=I-β－１uuＴ,β=uＴu. 同样可验证Hx=y成立.

例1 设R3中向量x=(0,4,3)Ｔ和向量σe1，其中σ2=‖x‖２２，即σ= =±5 例1 设R3中向量x=(0,4,3)Ｔ和向量σe1，其中σ2=‖x‖２２，即σ= =±5. 按(11)式构造一个3阶Householder矩阵H，使Hx=σe1. 解:令u=x-σe1.设σ=-5, 于是 u=(-5, 4, 3)Ｔ 又由uＴu=25+16+9=50, 所以β=25. 由H=I-β-１uuＴ, 得

容易验证Hx=σe1. 在Rn中，如果已知非零向量x=(x1,x2,…，xn)Ｔ和向量σe1≠x, 其中，σ=±‖x‖2, 即‖σe1‖2=‖x‖2 由定理6.3.6， 存在一个形如(11)式或(10)式的 Householder矩阵H，使Hx=σe1. 由定理6.3.6， 令u=x-σe1 (13) 又由σ2=‖x‖22＝xＴx和(12)式， β=1/2× (x-σe1)Ｔ(x-σe1)= 1/2× (xＴx-2σx１+σ2) =σ(σ-x1). (14) 记H=I-β-１uuＴ，容易验证Hx=σe1.