曲线拟合——最小二乘法算法

曲线拟合——最小二乘法算法
一、目的和要求
1）了解最小二乘法的基本原理，熟悉最小二乘算法；
2）掌握最小二乘进行曲线拟合的编程，通过程序解决实际问题。

二、实习内容
1）最小二乘进行多项式拟合的编程实现。

2）用完成的程序解决实际问题。

三、算法
1）输入数据节点数n ，拟合的多项式次数m ，循环输入各节点的数据x j , y j (j=0,1,…,n-1)
2）由x j 求S ；由x j ,y j 求T ：
S k =
∑-=10n j k j x ( k=0,1,2, … 2*m ) T k = ∑-=1
0n j k j j x y ( k=0,1,2,… m )
3）由S 形成系数矩阵数组c i,j ：c[i][j]=S[i+j] (i=0,1,2,…m, j=0,1,2,…,m)；由T 形成系数矩阵增广部分c i,m+1：c[i][m+1]=T[i] (i=0,1,2,…m)
4）对线性方程组CA=T[或A C ]，用列主元高斯消去法求解系数矩阵A=(a 0,a 1,…,a m )T
四、实验步骤
1）完成最小二乘法进行曲线拟合的程序设计及录入、编辑；
2）完成程序的编译和链接，并进行修改；
3）用书上P105例2的例子对程序进行验证，并进行修改；
4）用完成的程序求解下面的实际问题。

5）完成实验报告。

五、实验结果
1. 经编译、链接及例子验证结果正确的源程序：
#include<stdio.h>
#include<math.h>
#define Q 100
float CF(int,float);
main()
{
int i,j,n1,n,p,k,q;
float x[Q],y[Q],s[Q]={0},t[Q]={0},a[Q][Q]={0},l,sum=0;
/*以下是最小二乘的程序*/
printf("input 数据组数n");
scanf("%d",&n);
printf("input 拟合次数n1");
scanf("%d",&n1);
for(i=0;i<n;i++)
{
printf("x[%d]=",i);
scanf("%f",&x[i]);
printf("y[%d]=",i);
scanf("%f",&y[i]);
}
for(i=0;i<=2*n1;i++)
for(j=0;j<n;j++)
{
s[i]=s[i]+CF(i,x[j]);
if(i<=n1)
t[i]=t[i]+y[j]*CF(i,x[j]);
}
for(i=0;i<n1+1;i++)
for(j=0;j<n1+2;j++)
{
a[i][j]=s[i+j];
if(j==n1+1)
a[i][j]=t[i];
}
for(i=0;i<n1+1;i++)
for(j=0;j<n1+2;j++)
printf("a[%d][%d]=%f",i,j,a[i][j]); /*以下的是削去法的程序*/
for(j=0;j<=n1-1;j++)
{p=j;
for(i=j+1;i<=n1;i++)
{
if(fabs(a[j][j])<fabs(a[i][j]))
p=i;
}
if(p!=j)
for(i=j;i<=n1+1;i++)
{l=a[p][i];
a[p][i]=a[j][i];
a[j][i]=l;
}
for(k=j+1;k<=n1;k++)
{l=a[k][j]/a[j][j];
for(q=j;q<=n1+1;q++)
a[k][q]=a[k][q]-l*a[j][q];
}
}
for(i=0;i<n1+1;i++)
{for(j=0;j<n1+2;j++)
printf("a[%d][%d]=%f\n",i,j,a[i][j]);
printf("\n");}
x[n1]=a[n1][n1+1]/a[n1][n1];
for(i=n1-1;i>=0;i--)
{for(j=i+1;j<=n1;j++)
sum=a[i][j]*x[j]+sum;
x[i]=(a[i][n1+1]-sum)/a[i][i];
sum=0;
}
for(i=0;i<=n1;i++)
printf("x[%d]=%f\n",i,x[i]);
}
float CF(int i,float v)
{float a=1.0;
while(i--)a*=v;
return a;
}
2. 实例验证结果：
1）输入初始参数：
n=9，m=2
X：1 3 4 5 6 7 8 9 10
Y：10 5 4 2 1 1 2 3 4
2）结果输出：
1.实际应用
问题：
作物体运动的观测实验，得出以下实验测量数据，用最小二乘拟合求物体的运动方程。

时间t(秒) 0 0.9 1.9 3.0 3.9 5.0 距离s(cm) 0 10 30 50 80 110
解题步骤：
1）画草图
2）确定拟合方程次数为1：
用完成的程序输入数据，求取拟合方程中的未知数，得出方程：
S=-7.855.58+22.253763T
计算误差：
3）确定拟合方程次数为2：
用完成的程序输入数据，求取拟合方程中的未知数，得出方程：
S=-0.583370+11.081389T+2.248811T*T
计算误差：
六、分析和讨论
结合实际问题，进行拟合次数的分析和讨论：
七、心得（*可选）
①调试过程中遇到的问题和解决对策；②经验体会等。

以上都是我业余所写，程序肯定有很多不足之处希望大家多多指教。