西北工业大学 计算方法课件 第五章 曲线拟合的最小二乘法 西工大 nwpu

ϕ
r
2.70 480
2.00 670
1.61 830
1.20 1080
1.02 1260
变形为: 解:变形为: 1 = 1 − e cosϕ, 则有如下数据
r
p
p
1 y= r
0.370370 0.669131
0.50000 0.390731
0.621118
0.83333
0.980392
t = cosϕ
其矩阵形式为 Ax = b 方程组的系数矩阵和增广矩阵的秩不相等时， 当方程组的系数矩阵和增广矩阵的秩不相等时， 方程组无解，此时方程组称为矛盾方程组 矛盾方程组。 方程组无解，此时方程组称为矛盾方程组。对于 的矛盾方程组（ ），我 rankA＝n（A的秩为n）的矛盾方程组（N>n），我 们寻求其最小二乘意义下的解。 们寻求其最小二乘意义下的解。
k
(2)矩阵
∂2 f 2 ∂x1 P 0 2 ∂ f M = ∂x2∂x1 P 0 ⋮ 2 ∂ f ∂x ∂x n 1 P0
P 0
∂2 f ∂x1∂x2
P 0
∂2 f 2 ∂x2 P 0 ⋮ ∂2 f ∂xn∂x2 P
0
P 0 2 ∂ f ⋯ ∂x2∂xn P 0 ⋱ ⋮ ∂2 f ⋯ 2 ∂xn P 0 ∂2 f ⋯ ∂x1∂xn
1.最小二乘原则 1.最小二乘原则
由于矛盾方程组的精确解不存在， 由于矛盾方程组的精确解不存在，我们转而 寻求其某种意义下，即最小二乘意义下的解。 寻求其某种意义下，即最小二乘意义下的解。 n 令 δi = ∑aij xj −bi (i =1,2,⋯, N) 偏差。 称 δi 为偏差。 工程实际中的许多问题都可以归结为矛盾方程组， 工程实际中的许多问题都可以归结为矛盾方程组， 实际中需要寻求矛盾方程组的一组解， 实际中需要寻求矛盾方程组的一组解，以使得偏差的 N 尽可能地小。 绝对值之和 ∑δi 尽可能地小。为了便于分析
曲线拟合问题的关键： 曲线拟合问题的关键：
选择合适的曲线类型： 选择合适的曲线类型：根据问题的物理规律或数 据特点， 据特点，选择合适函数空间 {φ0 ( x), φ1 ( x),⋯ , φm ( x)} ，则拟合曲线可以表示为
φ ( x) = c0φ0 ( x) + c1φ1 ( x) + ⋯ + cmφm ( x)(m < n)
2.最小二乘解的存在唯一性 2.最小二乘解的存在唯一性
引理1：设n元实函数f(x1,x2,…,xn)在点P0(a1,a2,…,an) 设 , , 的某个邻域内连续，且有一阶及二阶连续的偏导数， 的某个邻域内连续，且有一阶及二阶连续的偏导数，如 ∂f 果 =0 (k =1 2,⋯ n) , , (1) ∂x
是正（负）定矩阵，则f(a1,a2,…,an)是n元实函 是正（ 定矩阵， , 数f(x1,x2,…,xn)的极小（大）值。 , 的极小（
引理2 引理2：设非齐次线性方程组 Ax = b 的系数矩阵 A=(aij)N×n，若rankA=n，则 (1)矩阵 是对称正定矩阵； (1)矩阵ATA是对称正定矩阵； AT Ax = ATb 有唯一的解。 有唯一的解。 (2)n阶线性方程组 证明： 显然是对称矩阵。 证明：（1）矩阵ATA显然是对称矩阵。 设齐次线性方程组 Ax = 0 因为rankA=n，故齐次方程组有唯一零解。 因为 ，故齐次方程组有唯一零解。 因此， 因此，对于任意的 x ≠ 0，有Ax ≠ 0 从而 ， ( Ax)T ( Ax) = xT ( AT A)x > 0 是对称正定矩阵。 故矩阵ATA是对称正定矩阵。 (2)因为矩阵 是正定矩阵， (2)因为矩阵ATA是正定矩阵，故rank(ATA)=n，从 有唯一的解。 而线性方程组 AT Ax = ATb 有唯一的解。 证毕
3.最小二乘法解矛盾方程组 3.最小二乘法解矛盾方程组
计算步骤： 判断方程组的秩是否满足rank （1）判断方程组的秩是否满足rankA=n？ 写出正则方程组； （2）写出正则方程组； 求解正则方程组， （3）求解正则方程组，其解就是矛盾方程组 的最小二乘解。 的最小二乘解。
§5.3 曲线最小二乘拟合
(φk ,φ j ) = ∑φk ( xi )φ j ( xi ),
i =0
n
(φk , b) = ∑φk ( xi )yi
i =0
n
Remark
同一问题可以有不同的拟合曲线， 1.同一问题可以有不同的拟合曲线，通常根据均方误 差
ax 和最大偏差 m≤N ϕ(xi ) − yi 1≤i ∑[ϕ(xi ) − yi ]
Remark
中的待定常数是线性形式时， 3.当拟合曲线ϕ(x)中的待定常数是线性形式时，可直 接根据矛盾方程组得到正则方程组而求解。 接根据矛盾方程组得到正则方程组而求解。当待定常 数不是线性形式时，则应该先将待定常数线性化， 数不是线性形式时，则应该先将待定常数线性化，再 根据矛盾方程组写出正则方程组而求解。 根据矛盾方程组写出正则方程组而求解。
设线性方程组
a11x1 + a12x2 +⋯+ a1n xn = b1 a x + a x +⋯+ a x = b 21 1 22 2 2n n 2 ⋮ aN1x1 + aN2 x2 +⋯+ aNn xn = bN
或写为
∑a x
j =1 ij
n
j
= bi
( j =1,2,⋯, N)
0.121869 -0.309017 -0.587785
记
a=
1 e , b = − ，得拟合模型：a + bt 得拟合模型： p 1 1 1 0.669131 0.370370 0.390731 0.500000 a = 0.621118 0.121869 b − 0.309017 0.833333 0.980392 − 0.587785
明显该方程组无解，是矛盾方程组， 明显该方程组无解，是矛盾方程组，可以寻 求其在最小二乘意义下的解。 求其在最小二乘意义下的解。对应的正规方程 组为
(φ0 , φ0 ) (φ0 , φ1 ) ⋯ (φ0 , φm ) c0 (φ0 , b) (φ , φ ) (φ , φ ) ⋯ (φ , φ ) c (φ , b) 1 1 1 m 1 1 0 = 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ (φm , φ0 ) (φm , φ1 ) ⋯ (φm , φm ) cm (φm , b)
引理2说明,在条件RankA=n下,无论线性方程组Ax=b 引理2说明,在条件RankA=n下 无论线性方程组Ax=b RankA=n 是否有解,构造的n阶方程组ATAx=ATb一定有唯一解。 是否有解,构造的n阶方程组A 一定有唯一解。
定理： 定理：设矛盾方程组的系数矩阵的秩为n，则二次 函数 2 N n Q = f (x1, x2 ,⋯ xn ) = ∑ ∑aij xj −bi , i= j =1 1 一定存在最小值， 一定存在最小值，且最小值点为方程组 AT Ax = ATb 的解。 的解。 Remark1： 称为正则方程组。 Remark1：线性方程组 AT Ax = ATb称为正则方程组
i= 1
j =1
计算和应用， 计算和应用，常采用使偏差的平方和 达到最小值，这一条件称为最小二乘原则。 达到最小值，这一条件称为最小二乘原则 最小二乘原则
∑aij xj − bi Q = ∑δ = ∑ i=1 i= j =1 1
N N n 2 i 2
按照最小二乘原则来选择未知数x1,x2,…,xn的 , 一组取值的方法称为求解矛盾方程组的最小二乘法 最小二乘法。 一组取值的方法称为求解矛盾方程组的最小二乘法。 符合条件的一组取值称为矛盾方程组的最小二乘解 最小二乘解。 符合条件的一组取值称为矛盾方程组的最小二乘解。 把Q看成是n个自变量x1,x2,…,xn的二次函数， , 的二次函数， 因此， 记为Q＝f(x1,x2,…,xn)，因此，求矛盾方程组的 , 最小二乘解就是求二次函数Q＝f(x1,x2,…,xn)的 , 最小值点。 最小值点。 问题：二次函数 问题 二次函数Q＝f(x1,x2,…,xn)是否存在最小 , 若最小值存在，如何求出该最小值点？ 值？若最小值存在，如何求出该最小值点？
y = aebx 例1： ： ln y = ln a + bx u = ln y, A = ln a, B = b u = A+ Bx 1 1 y= 例2： ： = a + bx a + bx y 1 u = a + bx u= y
例1. 对彗星1968Tentax的移动在某极坐标系下有如 对彗星1968Tentax的移动在某极坐标系下有如 1968Tentax 下表所示的观察数据，假设忽略来自行星的干扰, 下表所示的观察数据，假设忽略来自行星的干扰,坐 p r 标应满足： 其中： 为参数， 为偏心率, 标应满足： = 1− ecosϕ 其中：p为参数，e为偏心率, 试用最小二乘法拟合p和e。
设曲线拟合模型为
φ ( x) = c0φ0 ( x) + c1φ1 ( x) + ⋯ + cmφm ( x)(m < n)
如果试图插值，即 φ ( xi ) = yi (i = 0,1, 2,⋯ , n) 如果试图插值， 则得到关于未知数 c0 , c1 ,⋯ , cm的线性方程组
c0φ0 ( x0 ) + c1φ1 ( x0 ) + ⋯ + cmφm ( x0 ) = y0 c φ (x ) + c φ (x ) +⋯ + c φ (x ) = y 0 0 1 1 1 1 m m 1 1 ⋮ ⋮ ⋮ ⋮ c0φ0 ( xn ) + c1φ1 ( xn ) + ⋯ + cmφm ( xn ) = yn