最小二乘法-cpp.txt

最小二乘法-cpp.txt

#include<iostream.h> 
#include<math.h> 

void main() 
{ 
int i; 
float *a; 
float x[16]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16}; 
float y[16]={4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86, 
10.00,10.20,10.32,10.42,10.50,10.55,10.58,10.60}; 
float *Approx(float *,float *,int ,int ); 
a=Approx(x,y,16,2); 
for(i=0;i<=2;i++) 
cout<<"a["<<i<<"]="<<a[i]<<endl; 
} 
float *Approx(float *x,float *y,int m,int n) 
{ 
float *c,*a; 
int i,j,t; 
float power(int ,float); 
float *ColPivot(float *,int ); 
c= new float[((n+1)*(n+2)*sizeof(float))]; 
for(i=0;i<=n;i++) 
{ 
for (j=0;j<=n;j++) 
{ 
*(c+i*(n+2)+j)=0.0; 
for(t=0;t<=m-1;t++) 
*(c+i*(n+2)+j)+=power(i+j,x[t]); 
} 
*(c+i*(n+2)+n+1)=0.0; 
for(j=0;j<=m-1;j++) 
*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]); 
} 
a=ColPivot((float *)c,n+1); 
return a; 
} 
float *ColPivot(float *a,int n) 
{ 
int i,j,t,k; 
float *x,*c,p; 
x=new float[(n*sizeof(float))]; 
c=new float[(n*(n+1)*sizeof(float))]; 
for(i=0;i<=n-1;i++) 
for (j=0;j<=n;j++) 
*(c+i*(n+1)+j)=(*(a+i*(n+1)+j)); 
for(i=0;i<=n-2;i++) 
{ 
k=i; 
for(j=i+1;j<=n-1;j++) 
if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i)))) 
k=j; 
if(k!=j) 
for(j=i;j<=n;j++) 
{ 
p=*(c+i*(n+1)+j); 
*(c+i*(n+1)+j)=*(c+k*(n+1)+j); 
*(c+k*(n+1)+j)=p; 
} 
for(j=i+1;j<=n-1;j++) 
{ 
p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i)); 
for(t=i;t<=n-1;t++) 
*(c+j*(n+1)+t)=*(c+j*(n+1)+t)-p*(*(c+i*(n+1)+t)); 
*(c+j*(n+1)+n)-=*(c+i*(n+1)+n)*p; 
} 
} 
for( i=n-1;i>=0;i-- ) 
{ 
for( j=n-1;j>=i+1;j--) 
(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j)); 
x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i)); 
} 
delete[] c; 
return x; 
} 
float power(int i,float v) 
{ 
float a=1.0; 
while(i--) 
a*=v; 
return a; 
}