曲线拟合最小二乘法

Code:

1. #include <stdio.h>   
2. #include <stdlib.h>   
3. #include <math.h>   
4. #include <malloc.h>   
5.   
6. int main()   
7. {   
8.     int i;   
9.     float *a;   
10.     float x[16] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};   
11.     float y[16] = {4.00,6.40,8.00,8.80,9.22,9.50,9.70,   
12.                     9.86,10.00,10.20,10.32,10.42,10.50,10.55,   
13.                     10.58,10.60};   
14.     float  *Approx(float *,float *,int ,int );   
15.     a =Approx(x,y,16,2);   
16.   
17.     for(i = 0; i<=2; i++)   
18.     {   
19.         printf("a[%d]=%f/n",i,a[i]);   
20.     }   
21.     return 0;   
22. }   
23.   
24. float * Approx(float * x,float * y,int m,int n)   
25. {   
26.     float * c,*a;   
27.     int i,j,t;   
28.     float power(int, float);   
29.     float *ColPivot(float *,int);   
30.     c = (float *)malloc((n+1)*(n+2)*sizeof(float));   
31.     for(i = 0; i<=n; i++)   
32.     {   
33.         for(j =0; j<=n ; j++)   
34.         {   
35.             *(c+i*(n+2)+j)=0.0;   
36.             for(t=0; t<=m-1;t++)   
37.                *(c+i*(n+2)+j)+= power(i+j,x[t]);   
38.         }   
39.         *(c+i*(n+2)+n+1)=0.0;   
40.         for(j=0 ; j<m-1; j++)   
41.                *(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);   
42.     }   
43.     a= ColPivot((float *)c,n+1);   
44.     return a;   
45. }   
46.   
47. float *ColPivot(float *a,int n)   
48. {   
49.     int i,j,t,k;   
50.     float *x,*c,p;   
51.     x = (float *)malloc(n*sizeof(float));   
52.     c = (float *)malloc(n*(n+1)*sizeof(float));   
53.     for(i = 0; i<=n-1; i++)   
54.     for(j=0;j<=n;j++)   
55.     {   
56.         *(c+i*(n+1)+j)=(*(a+i*(n+1)+j));   
57.     }   
58.     for(i=0; i<=n-2; i++)   
59.     {   
60.   
61.         k=i;   
62.         for(j=i+1;j=<n-1; j++)    //添加等号<=n-1
63.         {   
64.             if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))   
65.               k=j;   
66.         }   
67.         if(k != i)   
68.            for(j =i; j<=n ; j++)   
69.            {   
70.                p = *(c+i*(n+1)+j);   
71.                *(c+i*(n+1)+j) = *(c+k*(n+1)+j);   
72.                *(c+k*(n+1)+j) = p;   
73.            }   
74.         for(j=i+1;j<=n-1;j++)   
75.         {   
76.             p=(*(c+j*(n+1)+i))   
77.             /(*(c+i*(n+1)+i));   
78.             for(t=i; t<=n-1; t++)   
79.             {   
80.                 *(c+j*(n+1)+t)   
81.                     = *(c+j*(n+1)+t)   
82.                     -p*(*(c+i*(n+1)+t));   
83.             }   
84.           *(c+j*(n+1)+n) -= *(c+i*(n+1)+n)*p;   
85.         }   
86.     }   
87.     for(i=n-1; i>=0;i--)   
88.     {   
89.         for(j=n-1;j>=i+1;j--)   
90.         {   
91.             (*(c+i*(n+1)+n)) -=   
92.                x[j]*(*(c+i*(n+1)+j));   
93.         }   
94.         x[i] =*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));   
95.     }   
96.     free(c);   
97.     return x;   
98. }   
99.   
100. float power(int i,float v)   
101. {   
102.     float a=1.0;   
103.     while(i--) a*=v;   
104.     return a;   
105. }