离散试验数据点的正交多项式最小二乘拟合
数值分析--32正交多项式与最小二乘拟合
j 0 n
{ j(x) = x j } 对应代数多项式 /* algebraic polynomial */
{ j(x) = cos jx }、{ j(x) = sin jx } { j(x), j(x) }对应三 角多项式 /* trigonometric polynomial */ { j(x) = e kj x , ki kj } 对应指数多项式 /* exponential
n
定理 Ba = c 存在唯一解 0(x), 1(x), … , n(x) 线性无关。
证明：若存在一组系数 {i } 使得 0 0 + 1 1 + ... + n n 0 则等式两边分别与0, 1, … , n作内积，得到：
0 ( 0 , 0 ) + 1 ( 1 , 0 ) + ... + n ( n , 0 ) 0 ( , ) + ( , ) + ... + ( , ) 0 1 1 1 n n 1 0 0 1 . . . ( , ) + ( , ) + ... + ( , ) 0 1 1 n n n n 0 0 n
1 2 49 3 y P( x) x + x 2 10 2
cond ( B ) 7623
|| B || 484,
|| B
1
||
63 4
§6 Orthogonal Polynomials & L-S Approximation
j 例：连续型拟合中，取 j ( x) x , ( x) 1, y( x) C[0, 1]
3.4 离散数据的曲线拟合——数值分析课件PPT
4
(P0 , P0 ) i P02 (xi ) 5 i0
4
(xP0 , P0 ) i xi P02 (xi ) 2.5 i0
a0
(xP0 , P0 ) (P0 , P0 )
0.5
P1(x) x a0 x 0.5
由此得 从而有
4
(P1, P1) i P12 (xi ) 0.625 i0
aj j (x)存在唯一；
j0
(b) p *(x)
n
aj j (x)的系数
a
j
n 可由法方程组
j0
j0
(0 ,0 ) (1 ,0 )
(n ,0 )
(0 ,1 ) (0 ,n ) a0 ( f ,0 )
(1 ,1 )
(1
,n
)
a1
( f
,1
)
(n ,1 )
(n ,n )an
i1
m
xi
i 1
m
xi2
i 1
m
xi3
i 1
m
xi
2
m
yi
i1 m
xi
3
a0
a1
i 1
m
xi yi
i1 m
xi 4
a2
i1
i1
m
xi
2
yi
i1
例 3.4.1 用多项式拟合表3-4中的离散数据。
表3-4
i
1
2
3
45
xi 0.00 0.25 0.50 0.75 1.00 yi 0.10 0.35 0.81 1.09 1.96
(
f
,
n
)
或Ga
d
最小二乘法多项式拟合实验报告.docx
连续系统仿真实验报告实验数据拟合建模姓名：专业：学号：时间：2013年5月1日实验单元二实验数据拟合建模一、实验目的1、 用C 语言实现最小二乘的多项式拟合和LU 分解法；2、 熟练掌握最小二乘拟合和LU 分解法的基本原理。
3、 体会用计算机编程解决计算问题的方法。
二、需求说明（一） 、需求阐述本次实验是要求根据己知的自变量和函数值，通过多项式拟合來分别计 算2、3、4阶拟合多项式，并根据拟合结果分別计算出待求点的函数值。
其中解 拟合系数方程组时采用LU 分解的方法计算拟合多项式的系数。
（二） 、实验公式m 次拟合函数公式为：（p （x ）=ao 七1対~・・・可点"计算系数4的方程组为:Sg a 0 +S] a 】 +...4-s ni a ni =t 0 < S]a ()+s?a]+...+s mF ]a m =t]k Sm a 0 +S mH a i +• •丹加^冃 其中 》= 士疋E ，i-0所以，在编程计算时，先计算出方程组①，再用LU 分解法计算求出耳的 值，即可得到拟合多项式。
LU 分解法的公式为：其中L 矩阵和U 矩阵的计算公式如下: 第一步，当k 二1,有：「1 0 0・・・01〔21 1 0-0^31 彳32 1 ••::::0 厶L ……1-i=0n-l最后求 u nn : U nn =a nn -^l m u m r=l三、设计说明（一） 、数据结构程序采用一维数组的形式来读取文件中给出的己知点处的值和要计算的未 知点处的H 变量值，最终的拟合计算结果也是采用一维数组的形式输出到文件中。
拟合多项式的系数a 和拟合系数方程组的参数t 都是采用一维数组來存储的，而 拟合系数方程组中的参数s 和L 、U 矩阵都是用二维数组來表示的。
由于要分别 计算2、3、4阶拟合结果，所以数组的规模取为5,矩阵的规模取为5*5.（二） 、算法设计及效率分析在进行LU 分解函数中，在计算L 矩阵和U 矩阵时，因为当k=2,3.-,n-l 时， 计算丈M 和土皿的循环条件不允许k=l 时进入,而正好k=l 时，计算1“和i 町不 x-1 r-1k-1 k ・l需要工1丿匕和工1以崎，因而对k=l 和k=2,3,-,n-l,就可以和在一起计算，这样就减少了 r=lr=l程序的长度。
实验题目用正交多项式做小二乘曲线拟合
实验题目：用正交多项式做小二乘曲线拟合实验题目： 用正交多项式做最小二乘的曲线拟合 学生组号： 6 完成日期： 2011/11/27 1 实验目的针对给定数据的煤自燃监测数据中煤温与NO 22,之间的非线性关系，用正交多项式做最小二乘曲线拟合。
2 实验步骤2.1 算法原理设给定n+1个数据点：（yx kk,），k=0,1,···,n ，则根据这些节点作一个m 次的最小二乘拟合多项式pm（x ）=a+x a x a a mm x +++ (2)21=x a jmj j ∑=0①其中，m ≤n,一般远小于n.。
若要构造一组次数不超过ｍ的在给定点上正交的多项式函数系｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝，则可以首先利用｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝作为基函数作最小二乘曲线的拟合，即pm（x ）＝)(...)()(11x x x Qq Q q Q q mm+++ ②根据②式，其中的系数qj（ｊ＝０，１，．．．，ｍ）为∑∑===nk kjnk kjkjx Q x Q y q2)()(，ｊ＝０，１，．．．，ｍ ③将④代入③后展开就成一般的多项式。
构造给定点上的正交多项式)(x Qj（ｊ＝０，１，．．．，ｍ）的递推公式如下：⎪⎪⎩⎪⎪⎨⎧-=--=-==-+1,...,2,1),()()()()()(1)(11010m j x x x x x x x QQ Q Q Q j jj j j βαα ④其中αj＝dx x jk j=0，ｊ＝０，１，．．．，ｍ－１ ⑤βj＝dd j j1-，ｊ＝１，２，．．．，ｍ－１ ⑥∑==nk k jjx Q d2)(，ｊ＝０，１，．．．，ｍ－１ ⑦则实际计算过程中，根据⑤式逐步求出个正交多项式)(x Qj，并用公式④计算出q j，并将每次计算展开后累加到拟合多项式①中。
2.2 算法步骤用三个向量Ｂ,Ｔ,Ｓ,存放多项式)(1x Qj -，)(x Q j，)(1x Qj +的系数。
数值分析(22)离散数据的最小二乘拟合
xi
m i0
(
m
n)上至多只有n个不同零点，则称0
(
x
),
1
(
x),
...,
n
(
x
)在点集
xi
m i0
上满足Haar条件。
可以证明，如果0
(
x
),
1
(
x),
...,
n
(
x
)在
xi
m i0
上满足
Haar条件，则法方程的系数矩阵非奇异。
数值分析
数值分析
法方程的另一种形式
记 A (0 , 1 , ..., n ), 由矩阵乘法可知 G ATWA，F ATWY
则
1
A
1
x0 x1
x02 x12
x0n x1n
1
xm
xm2
xmn
m
实际计算中要计算x0n，x1n，...，xmn , xik ,(k 0,1, ..., 2n)
i 1
当数据点数值很大时，可能会在计算中溢出，可采用
P240的方法处理。
数值分析
数值分析
三、法方程的求解
求解法方程AT AC ATY或GC F (取W I )
1.确定n+1个线性无关函数
j( x)
n 的具体形式
j=0
常用的代数多项式拟合是取
j( x)＝x j (j 0,1, , n) 最佳平方逼近函数的形式是
s( x)＝c0 c1 x c2 x2 cn xn
G
(1
最小二乘拟合多项式
则法方程系数矩阵为： T 常数项为：
T a y
T
n y i i 1 n x y i i y i 1 n xm y i i i 1
其他类型的拟合问题
最小二乘法并不只限于多项式，也可用于任何具体给出的 函数形式。特别重要的是有些非线性最小二乘拟合问题通 过适当的变换可以转化为线性最小二乘问题求解。
因此，我们需要一种新的逼近原函数的办法
解决方案：
1. 2. 不要求过所有数据点（可以消除误差影响）； 尽可能地刻画数据点的趋势，靠近这些数据点。
插值与拟合的关系： 问题：给定一组数据点，构造一个函数作为近似（或逼近）。 解决方案：
1. 若要求所求曲线通过给定的所有数据点，就是插值问题； 2. 若不要求曲线通过所有数据点，而是要求它反映数据点的整体变化趋势， 这就是数据拟合，又称曲线拟合，所求出的曲线称为拟合曲线。
a0 sn a1sn 1 an s2 n un
称为正规方程组（或法方程组）。
可以证明：当 x0 , x1, 是最小值问题的解。
法方程组 可写成以 下形式:
m m xk k 1 m n xk k 1
xn , 互异时，该方程组有唯一解，并
k
x
k 1 m k 1
m
2 x k
令
x
k 1
m
n 1 k
n x y i k 1 a0 i 1 m n xy n 1 x a k 1 i 1 i i k 1 an n m 2n xin yi xk k 1 i 1
曲线拟合
可做变换
Y ln y ,
X
1 x
,
A ln a ,
B b
Y A BX 就是个线性问题
将( xi , yi ) 化为( X i ,Yi ) 后易解 A 和B
a eA , b B , P(x) a eb/x
二、 一般的最小二乘法
2
2
m
2 i
m
m
[S*(
xi
)
yi
]2
min
S ( x )
则 勒 让 德 多 项 式 为:
P~n ( x)
n! (2n)!
dn dx n
[(x2
1)n ]
勒 让 德 多 项 式 有 以 下 几个 重 要 性 质 ：
性 质1： 正 交 性
1 1
Pn
(
x)Pm
(
x)d(
x)
0, 2
2n
1
,
m n; m n.
性 质2： 奇 偶 性 Pn ( x) (1)n Pn ( x) 性 质3： 递 推 关 系 (n 1)Pn1( x) (2n 1)xPn ( x) nPn1( x) 性 质4： 在 所 有 最 高 项 系 数 为1的n次 多 项 式 中 ， 勒 让 德 多项 式
条件.
可以证明,如果0(x), 1(x),… n(x)C[a,b]在{xi}0m上满
足哈尔(Haar)条件,则法方程(5.6)的系数矩阵G非奇异.
用最小二乘法得到的法方程组(3. 6)，其系数矩阵G是
病态的,但如果0(x), 1(x),… n(x)是关于点集 {xi}(i=0,
l, ..., m)带权(xi) (i=0,l,...,m) 正交的函数族，即
课堂练习
数据拟合的最小二乘方法的实现
数据拟合的最小二乘方法的实现班级学号姓名榴莲一、实验任务给定离散样本点，采用最小二乘方法拟合样本数据，涉及的线性方程组请用高斯列主元消去法求解。
求函数f(x)=cosπx,x∈[0,1]的一次和二次最佳平方逼近多项式。
二、编程环境Windows7，Codeblock.三、算法步骤1，采用目标函数对多项式系数求偏导，得到最优值条件，组成一个方程组；2，方程组的解法采用行列式变换（两次变换：普通行列式——三角行列式——对角行列式——求解），行列式的求解算法上优化过一次了，目前还没有更好的思路再优化运算方法，限幅和精度准备再修改修改。
四、程序流程图数据结构：五、程序#include<stdio.h>#include<math.h>double qwer(double a[],double b[],int n);void guess(double a[][10],double b[],double x[],int n); int main() {double a[10][10],b[10],asd[10];int i,j,k,n;int num = 6;double y[10] = {1.0,1.004,1.031,1.117,1.223,1.422};double x1[10] = {1.0,1.0,1.0,1.0,1.0,1.0};double x2[10] = {0,0.15,0.31,0.5,0.6,0.75};n = 2;a[0][0] = qwer(x1,x1,num);a[0][1] = qwer(x1,x2,num);a[1][0] = qwer(x2,x1,num);a[1][1] = qwer(x2,x2,num);b[0] = qwer(x1,y,num);b[1] = qwer(x2,y,num);printf("矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}guess(a,b,asd,n);printf("\n三角化矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}printf("\n系数值asd：");for(i=0; i<n; i++) printf("%8.4lf",asd[i]);printf("\n一次多项式为：f(x) = %.4lf+%.4lfx\n",asd[0],asd[1]);double x3[10] = {0*0,0.15*0.15,0.31*0.31,0.5*0.5,0.6*0.6,0.75*0.75};n = 3;a[0][0] = qwer(x1,x1,num);a[0][1] = qwer(x1,x2,num);a[0][2] = qwer(x1,x3,num);a[1][0] = qwer(x2,x1,num);a[1][1] = qwer(x2,x2,num);a[1][2] = qwer(x2,x3,num);a[2][0] = qwer(x3,x1,num);a[2][1] = qwer(x3,x2,num);a[2][2] = qwer(x3,x3,num);b[0] = qwer(x1,y,num);b[1] = qwer(x2,y,num);b[2] = qwer(x3,y,num);printf("\n\n矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}guess(a,b,asd,n);printf("\n三角化矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}printf("\n系数值asd：");for(i=0; i<n; i++) printf("%8.4lf",asd[i]);printf("\n二次多项式为：f(x) = %.4lf%.4lfx+%.4lfx^2\n",asd[0],asd[1],asd[2]);/* while(1) {double cc;scanf("%lf",&cc);printf("1=%lf\n",asd[0]+asd[1]*cc+asd[2]*cc*cc);printf("2=%lf\n",0.8115+1.3418*cc-0.9262*cc*cc);}*/}double qwer(double a[],double b[],int n) {int i,j;double s = 0;for(i=0; i<n; i++) {s += a[i]*b[i];}return s;}void guess(double a[][10],double b[],double x[],int n) { int k,i,j;//guess法化下三角矩阵for(k=0; k<n-1; k++) {//找主元double ma = a[k][k];int tab = k;for(i=k+1; i<n; i++) {if(fabs(ma)<fabs(a[i][k])) {ma = a[i][k];tab = i;}}double mid;mid = b[k];b[k] = b[tab];b[tab] = mid;for(i=k; i<n; i++) {mid = a[k][i];a[k][i] = a[tab][i];a[tab][i] = mid;}//三角化for(i=k+1; i<n; i++) {b[i]=b[i]-a[i][k]/a[k][k]* b[k];for(j=k+1; j<n; j++) {a[i][j]=a[i][j]-a[i][k]/a[k][k]*a[k][j];}a[i][k]=0;}}//回代求解x值for(k=n-1; k>=0; k--) {double s =0;for(j=k+1; j<n; j++) s += a[k][j]*x[j];x[k]=(b[k]-s)/a[k][k];}}六、实验结果及分析通过数据可以看到asd: 0,9295 0.5281。