最小二乘法用于多项式的拟合及程序实现

改写自：https://blog.csdn.net/piaoxuezhong/article/details/54973750

  1 #include <stdio.h>
  2 #include "stdlib.h"
  3 #include "math.h"
  4 //#include <vector>
  5 //using namespace std;
  6 
  7 #define n 4 /* x最高次幂 */
  8 #define N (n + 1)
  9 #define SMPNUM 5    /* 采样个数 */
 10  
 11 struct point
 12 {
 13     double x;
 14     double y;
 15 };
 16 
 17 float sampleX[SMPNUM] = {0.0};
 18 float sampleY[SMPNUM] = {0.0};
 19 float coeff[N] = {0};
 20 
 21 void getFile(char *File);  //获取文件数据
 22 void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM]);   //矩阵方程
 23  
 24 int main()
 25 {
 26     int i;
 27     char *File = "XY.txt";
 28     //vector<point> sample;
 29     //doubleVector  coefficient;
 30     //sample = getFile(File);
 31     getFile(File);
 32     printf("拟合多项式阶数n=");
 33     //scanf("%d", &n);
 34     // coefficient = getCoeff(sample, n);
 35     //n = 3;
 36     getCoeff(sampleX, sampleY);
 37     printf("\n拟合矩阵的系数为：\n");
 38     for (i = 0; i < N; i++)
 39         printf("a%d = %lf\n", i, coeff[i]);
 40 
 41     return 0;
 42 }
 43 //矩阵方程
 44 void getCoeff(float sampleX[SMPNUM], float sampleY[SMPNUM])
 45 {
 46     double sum;
 47     int i, j, k;
 48 
 49     float matX[N][N];  //公式3左矩阵
 50     float matY[N][N];  //公式3右矩阵
 51     float temp2[N];
 52     //公式3左矩阵
 53     for (i = 0; i <= n; i++)
 54     {
 55         for (j = 0; j <= n; j++)
 56         {
 57             sum = 0.0;
 58             for (k = 0; k < SMPNUM; k++)
 59             {
 60                 sum += pow(sampleX[k], j + i);
 61             }
 62             matX[i][j] = sum;
 63         }
 64     }
 65  
 66     //打印matFunX矩阵
 67     printf("matFunX矩阵：\n");
 68     for (i = 0;i <= n;i++) {
 69         for (j = 0;j <= n;j++)
 70             printf("%f\t", matX[i][j]);
 71         printf("\n");
 72     }
 73  
 74     //公式3右矩阵
 75     for (i = 0; i <= n; i++)
 76     {
 77         //temp.clear();
 78         sum = 0;
 79         for (k = 0; k < SMPNUM; k++)
 80         {
 81             sum += sampleY[k]*pow(sampleX[k], i);
 82         }
 83         matY[i][0] = sum;
 84     }
 85     //printf("matFunY.size=%d\n", matFunY.size());
 86     //打印matFunY的矩阵
 87     printf("matFunY的矩阵：\n");
 88     for (i = 0;i <=n;i++) 
 89     {
 90         printf("%f\n", matY[i][0]); 
 91     }
 92     //矩阵行列式变换，将matFunX矩阵变为下三角矩阵,将matFunY矩阵做怎样的变换呢？
 93     //AX=Y在将X变换为下三角矩阵X'时，是左右两边同乘ratio
 94     double num1, num2, ratio;
 95     for (i = 0; i <= n; i++)
 96     {
 97         num1 = matX[i][i];
 98         for (j = i + 1; j <= n; j++)
 99         {
100             num2 = matX[j][i];
101             ratio = num2 / num1;
102             for (k = 0; k <= n; k++)
103                 matX[j][k] = matX[j][k] - matX[i][k] * ratio;
104             matY[j][0] = matY[j][0] - matY[i][0] * ratio;
105         }
106     }
107     //打印matFunX行列变化之后的矩阵
108     printf("matFunX行列变换之后的矩阵：\n");
109     for (i = 0;i <=n;i++) 
110     {
111         for (j = 0;j <= n;j++)
112         {
113             printf("%f\t", matX[i][j]);
114         }
115         printf("\n");
116     }
117     //打印matFunY行列变换之后的矩阵
118     printf("matFunY行列变换之后的矩阵：\n");
119     for (i = 0;i <=n;i++) 
120     {
121         printf("%f\n", matY[i][0]);
122     }
123     //计算拟合曲线的系数
124     //doubleVector coeff(n + 1, 0);
125     for (i = n; i >= 0; i--)
126     {
127         if (i == n)
128         {
129             coeff[i] = matY[i][0] / matX[i][i];
130         }
131         else
132         {
133             for (j = i + 1; j <= n; j++)
134             {
135                 matY[i][0] = matY[i][0] - coeff[j] * matX[i][j];
136             }
137             coeff[i] = matY[i][0] / matX[i][i];
138         }
139     }
140     return;
141 }
142  
143 //获取文件数据
144 void getFile(char *File)
145 {
146     int i = 0, j = 0, k = 0;
147     //vector<point> dst;
148 
149     FILE *fp = fopen(File, "r");
150  
151     if (fp == NULL)
152     {
153         printf("Open file error!!!\n");
154         exit(0);
155     }
156  
157     point temp;
158     double num;
159  
160     while (fscanf(fp, "%lf", &num) != EOF)
161     {
162         if (i % 2 == 0)
163         {
164             temp.x = num;
165             sampleX[j++] = num;
166 
167             
168         }
169         else
170         {
171             temp.y = num;
172             sampleY[k++] = num;
173         }
174         i++;
175     }
176     fclose(fp);
177     return;
178     //return dst;
179 }