最小二乘法程序

最小二乘法程序

#include

#include

#include

#include

#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{

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;ifor(j=0;j{

res=0;

for(k=0;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{

for(j=0;j{

b[j][i]=a[i][j];

}

}

}

void init(float x_y[][2],int n)

{

int i;

printf(\"请输入%d个已知点：\\n\

for(i=0;i{

printf(\"(x%d y%d):\

scanf(\"%f %f\

}

}

void get_A(float matrix_A[][T+1],float x_y[][2],int n)

{

int i,j;

for(i=0;i{

for(j=0;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{

for(j=0;j{

printf(\"%-g\j]);

}

printf(\"\\n\");

}

}

void convert(float argu[][T+2],int n)

{

int i,j,k,p,t;

float rate,temp;

for(i=1;i{

for(j=i;j{

if(argu[i-1][i-1]==0)

{

for(p=i;p{

if(argu[p][i-1]!=0)

break;

}

if(p==n)

{

printf(\"方程组无解!\\n\");

exit(0);

}

for(t=0;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{

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{

temp=0;

for(j=0;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{

for(j=0;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\

for(i=0;i{

printf(\"a[%d]=%g,\

}

printf(\"\\n\");

printf(\"拟合曲线方程为:\\ny(x)=%g\

for(i=1;i{

printf(\" + %g\

for(j=0;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();

}

]]>

2007-4-19 19:23:57

一级(初级)

user1

100

403872

78010

1526752

jiangxc2004

0

你可以改一下

不从终端输入,直接在程序中给出参数

请输入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 1 -1 1 0 1 1 1 2 法矩阵为：

5 0 0 10

4 -8

1 -1

0 0

1 1

4 8

10 0

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

]]>

2007-4-19 19:26:11

一级(初级)

user1

100

40390406

78010

1526752

jiangxc2004

0

#include

void main ()

{

int num,i;

float x,y,l,m,n,p,a,b;

i=1;

l=0.0;

m=0.0;

n=0.0;

p=0.0;

printf (\"请输入你想计算的x，y的个数：\");

scanf(\"%d\

if (num>=1)

{

while (i<=num);

{

printf(\"请输入x的值\");

scanf (\"%lf\

printf(\"请输入y的值\");

scanf (\"%lf\

l+=x;

m+=y;

n+=x*y;

p+=x*x;

i++;

}

a=(num*n-l*m)/(num*p-l*l);

b=(p*m-n*l)/(num*p-l*l);

printf(\"最小二乘法所算得的斜率和截距分别为%f和%f\\n\

}

else printf(\"mun\"输入有误！);

}

#include

#include

//////////////////////////////////////////////////////////////////////////////////////////

//矩阵结构体

struct Matrix

{

int m,n;//m为行数,n为列数

double **pm;//指向矩阵二维数组的指针

};

//初始化矩阵mt，并置矩阵的行为m，列为n

void InitMatrix(struct Matrix *mt,int m,int n)

{

int i;

//指定矩阵的行和列

mt->m=m;

mt->n=n;

//为矩阵分配内存

mt->pm=new double *[m];

for(i=0;i{

mt->pm[i]=new double[n];

}

}

//用0初始化矩阵mt，并置矩阵的行为m，列为n

void InitMatrix0(struct Matrix *mt,int m,int n)

{

int i,j;

InitMatrix(mt,m,n);

for(i=0;ifor(j=0;jpm[i][j]=0.0;

}

//销毁矩阵mt

void DestroyMatrix(struct Matrix *mt)

{

int i;

//释放矩阵内存

for(i=0;im;i++)

{

delete []mt->pm[i];

}

delete []mt->pm;

}

//矩阵相乘mt3=mt1*mt2

//成功返回1，失败返回0

int MatrixMul(Matrix *mt1,Matrix *mt2,Matrix *mt3)

{

int i,j,k;

double s;

//检查行列号是否匹配

if(mt1->n!=mt2->m||mt1->m!=mt3->m||mt2->n!=mt3->n) return 0;

for(i=0;im;i++)

for(j=0;jn;j++)

{

s=0.0;

for(k=0;kn;k++) s=s+mt1->pm[i][k]*mt2->pm[k][j];

mt3->pm[i][j]=s;

}

return 1;

}

//矩阵转置mt2=T(mt1)

//成功返回1，失败返回0

int MatrixTran(Matrix *mt1,Matrix *mt2)

{

int i,j;

//检查行列号是否匹配

if(mt1->m!=mt2->n||mt1->n!=mt2->m) return 0;

for(i=0;im;i++)

for(j=0;jn;j++) mt2->pm[j][i]=mt1->pm[i][j];

return 1;

}

//控制台显示矩阵mt

void ConsoleShowMatrix(Matrix *mt)

{

int i,j;

for(i=0;im;i++)

{

printf(\"\\n\");

for(j=0;jn;j++) printf(\"%2.4f \

}

printf(\"\\n\");

}

//////////////////////////////////////////////////////////////////////////////////////////

//Guass列主元素消元法求解方程Ax=b,a=(A,b)

int Guassl(double **a,double *x,int n)

{

int i,j,k,numl,*h,t;

double *l,max;

l=new double[n];

h=new int[n];

for(i=0;ifor(i=1;i{

max=fabs(a[h[i-1]][i-1]);

numl=i-1;

//列元的最大值

for(j=i;j{

if(fabs(a[h[j]][i-1])>max)

{

numl=h[j];

max=fabs(a[h[j]][i-1]);

}

}

if(max<0.00000000001) return 0;

//交换行号

if(numl>i-1)

{

t=h[i];

h[i]=h[numl];

h[numl]=t;

}

for(j=i;jfor(j=i;jfor(k=i;k}

for(i=n-1;i>=0;i--)

{

x[i]=a[h[i]][n];

for(j=i+1;jx[i]=x[i]/a[h[i]][i];

}

//清除临时数组内存

delete []l;

delete []h;

return 1;

}

///////////////////////////////////////////////////////////////////////////////////////////

//最小二乘法求解矩阵Ax=b

int MinMul(Matrix A,Matrix b,double *x)

{

int i,j;

if(b.n!=1) return 0;

if(A.m!=b.m) return 0;

Matrix TranA;//定义A的转置

InitMatrix0(&TranA,A.n,A.m);

MatrixTran(&A,&TranA);

Matrix TranA_A;//定义A的转置与A的乘积矩阵

InitMatrix0(&TranA_A,A.n,A.n);

MatrixMul(&TranA,&A,&TranA_A);//A的转置与A的乘积

Matrix TranA_b;//定义A的转置与b的乘积矩阵

InitMatrix0(&TranA_b,A.n,1);

MatrixMul(&TranA,&b,&TranA_b);//A的转置与b的乘积

DestroyMatrix(&TranA);//释放A的转置的内存

Matrix TranA_A_b;//定义增广矩阵

InitMatrix0(&TranA_A_b,TranA_A.m,TranA_A.m+1);

//增广矩阵赋值

for(i=0;i{

for(j=0;jTranA_A_b.pm[i][TranA_A_b.m]=TranA_b.pm[i][0];

}

DestroyMatrix(&TranA_A);

DestroyMatrix(&TranA_b);

//Guass列主消元法求解

Guassl(TranA_A_b.pm,x,TranA_A_b.m);

DestroyMatrix(&TranA_A_b);

return 1;

}

int MinMul(double **A,double *b,int m,int n,double *x)

{

int r,i;

Matrix Al,bl;

Al.pm=new double *[m];

Al.m=m;

Al.n=n;

InitMatrix(&bl,m,1);

for(i=0;i{

Al.pm[i]=A[i];

bl.pm[i][0]=b[i];

}

r=MinMul(Al,bl,x);

delete []Al.pm;

DestroyMatrix(&bl);

return r;

}