80分急求：用c++写成的最小二乘法的源代码
80分急求：用c++写成的最小二乘法的源代码 问题描述：本人要做一个最小2乘法的曲线拟合的vc程序，现在已经得到了在假定坐标轴下面的5个点的坐标，分别存放在x[],y[]数组里面，现在需要根据已知点的坐标求出曲线拟合的2次多项式的3个系数，本来想用matlab去拟合的，可是在vc中如何使用matlab的polyfit函数不太会用，于是自己想找一个c++的最小二乘法的曲线拟合的函数自己调用，希望有过此类编程经验的高手能够帮忙救急，能够提供c语言的最小2乘法的函数，最好是纯c++（不是mfc代码），决不吝啬分数，解决即给分，谢谢了！！！！
[解决办法]
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <process.h>
#define N 5//N个点
#define T 3 //T次拟合
#define W 1//权函数
#define PRECISION 0.00001
float pow_n(float a,int n)
{
int i;
if(n==0)
return(1);
float res=a;
for(i=1;i <n;i++)
{
res*=a;
}
return(res);
}
void mutiple(float a[][N],float b[][T+1],float c[][T+1])
{
float res=0;
int i,j,k;
for(i=0;i <T+1;i++)
for(j=0;j <T+1;j++)
{
res=0;
for(k=0;k <N;k++)
{
res+=a[i][k]*b[k][j];
c[i][j]=res;
}
}
}
void matrix_trans(float a[][T+1],float b[][N])
{
int i,j;
for(i=0;i <N;i++)
{
for(j=0;j <T+1;j++)
{
b[j][i]=a[i][j];
}
}
}
void init(float x_y[][2],int n)
{
int i;
printf( "请输入%d个已知点：\n ",N);
for(i=0;i <n;i++)
{
printf( "(x%d y%d): ",i,i);
scanf( "%f %f ",&x_y[i][0],&x_y[i][1]);
}
}
void get_A(float matrix_A[][T+1],float x_y[][2],int n)
{
int i,j;
for(i=0;i <N;i++)
{
for(j=0;j <T+1;j++)
{
matrix_A[i][j]=W*pow_n(x_y[i][0],j);
}
}
}
void print_array(float array[][T+1],int n)
{
int i,j;
for(i=0;i <n;i++)
{
for(j=0;j <T+1;j++)
{
printf( "%-g ",array[i][j]);
}
printf( "\n ");
}
}
void convert(float argu[][T+2],int n)
{
int i,j,k,p,t;
float rate,temp;
for(i=1;i <n;i++)
{
for(j=i;j <n;j++)
{
if(argu[i-1][i-1]==0)
{
for(p=i;p <n;p++)
{
if(argu[p][i-1]!=0)
break;
}
if(p==n)
{
printf( "方程组无解!\n ");
exit(0);
}
for(t=0;t <n+1;t++)
{
temp=argu[i-1][t];
argu[i-1][t]=argu[p][t];
argu[p][t]=temp;
}
}
rate=argu[j][i-1]/argu[i-1][i-1];
for(k=i-1;k <n+1;k++)
{
argu[j][k]-=argu[i-1][k]*rate;
if(fabs(argu[j][k]) <=PRECISION)
argu[j][k]=0;
}
}
}
}
void compute(float argu[][T+2],int n,float root[])
{
int i,j;
float temp;
for(i=n-1;i> =0;i--)
{
temp=argu[i][n];
for(j=n-1;j> i;j--)
{
temp-=argu[i][j]*root[j];
}
root[i]=temp/argu[i][i];
}
}
void get_y(float trans_A[][N],float x_y[][2],float y[],int n)
{
int i,j;
float temp;
for(i=0;i <n;i++)
{
temp=0;
for(j=0;j <N;j++)
{
temp+=trans_A[i][j]*x_y[j][1];
}
y[i]=temp;
}
}
void cons_formula(float coef_A[][T+1],float y[],float coef_form[][T+2])
{
int i,j;
for(i=0;i <T+1;i++)
{
for(j=0;j <T+2;j++)
{
if(j==T+1)
coef_form[i][j]=y[i];
else
coef_form[i][j]=coef_A[i][j];
}
}
}
void print_root(float a[],int n)
{
int i,j;
printf( "%d个点的%d次拟合的多项式系数为:\n ",N,T);
for(i=0;i <n;i++)
{
printf( "a[%d]=%g, ",i+1,a[i]);
}
printf( "\n ");
printf( "拟合曲线方程为:\ny(x)=%g ",a[0]);
for(i=1;i <n;i++)
{
printf( " + %g ",a[i]);
for(j=0;j <i;j++)
{
printf( "*X ");
}
}
printf( "\n ");
}
void process()
{
float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_formu[T+1][T+2],y[T+1],a[T+1];
init(x_y,N);
get_A(matrix_A,x_y,N);
printf( "矩阵A为：\n ");
print_array(matrix_A,N);
matrix_trans(matrix_A,trans_A);
mutiple(trans_A,matrix_A,coef_A);
printf( "法矩阵为：\n ");
print_array(coef_A,T+1);
get_y(trans_A,x_y,y,T+1);
cons_formula(coef_A,y,coef_formu);
convert(coef_formu,T+1);
compute(coef_formu,T+1,a);
print_root(a,T+1);
}
void main()
{
process();
}
[解决办法]
你可以改一下
不从终端输入,直接在程序中给出参数
请输入5个已知点：
(x0 y0):-2 -0.1
(x1 y1):-1 0.1
(x2 y2):0 0.4
(x3 y3):1 0.9
(x4 y4):2 1.6
矩阵A为：
1 -2 4 -8
1 -1 1 -1
1 0 0 0
1 1 1 1
1 2 4 8
法矩阵为：
5 0 10 0
0 10 0 34
10 0 34 0
0 34 0 130
5个点的3次拟合的多项式系数为:
a[1]=0.408571, a[2]=0.391667, a[3]=0.0857143, a[4]=0.00833333,
拟合曲线方程为:
y(x)=0.408571 + 0.391667*X + 0.0857143*X*X + 0.00833333*X*X*X
[解决办法]
double CFunctionSimulate::CurveSimulate(double n[],double T[],int M,int N,int xi)
//参数说明//
//n[] 自变量数据数组
//T[] 因变量数据数组
//M 拟和公式的最高阶数
//N 数据组数
//xi 所返回系数值的对应参数
{
double b[10][20];
double A[10][10],B[10],t,x[10],y[10];
int i,j,k,l;
for(i=0;i <M;i++)
for(j=0;j <N;j++)
{
t=1;
for(l=0;l <i;l++)
{t=t*n[j];}
b[i][j]=t;
}
for(i=0;i <M;i++)
for(k=0;k <M;k++)
{
t=0;
for(j=0;j <N;j++)
t+=b[i][j]*b[k][j];
A[i][k]=t;
}
for(i=0;i <M;i++)
{
t=0;
for(j=0;j <N;j++)
t+=T[j]*b[i][j];
B[i]=t;
}
for(i=1;i <M;i++)
A[i][0]=A[i][0]/A[0][0];
for(i=1;i <M;i++)
for(j=i;j <M;j++)
{
t=0;
for(k=0;k <i;k++)
t+=A[k][j]*A[i][k];
A[i][j]=A[i][j]-t;
if(j+1!=M)
{
t=0;
for(k=0;k <i;k++)
t+=A[k][i]*A[j+1][k];
A[j+1][i]=(A[j+1][i]-t)/A[i][i];
}
}
y[0]=B[0];
for(i=1;i <M;i++)
{
t=0;
for(j=0;j <i;j++)
t+=A[i][j]*y[j];
y[i]=B[i]-t;
}
x[M-1]=y[M-1]/A[M-1][M-1];
for(i=M-2;i> =0;i--)
{
t=0;
for(j=i+1;j <M;j++)
t+=A[i][j]*x[j];
x[i]=(y[i]-t)/A[i][i];
}
return x[xi];
}
作毕业设计时搞的一个最小二乘法拟合曲线的函数，希望对楼主有用
[解决办法]
矩阵论里有公式的，但忘了，寒一个