【C语言】最小二乘法拟合一组数据：对于y=a+bx+cx^2拟合

注：本代码仅针对本公式：y=a bx cx^2拟合

#include<stdio.h> #include<math.h> //显示输入x和y void x_y_print(float x[],float y[],int n){ 
            system("cls");     printf("您输入的x为：\n");     printf("|\t");     for(int i=0;i<n;i ){ 
                printf("%f ",x[i]);     }     printf("\t|\n");     printf("您输入的y为：\n");     printf("|\t");     for(int i=0;i<n;i ){ 
                printf("%f ",y[i]);     }     printf("\t|\n"); } void main(){ 
            printf(/span>"请输入已知点个数：\n");//6 int n; scanf("%d",&n); printf("请输入未知数个数：\n");//3 int wn; scanf("%d",&wn); float x[n],y[n]; printf("请输入x的值:\n"); for(int i=0;i<n;i++){ 
         scanf("%f",&x[i]); } printf("请输入y的值:\n"); for(int i=0;i<n;i++){ 
         scanf("%f",&y[i]); } //x_y_print(x,y,n); float A[n][wn];//6行3列 for(int i=0;i<n;i++){ 
        //6 for(int j=0;j<wn;j++){ 
        //3 //printf("此时的值为：%lf\n",pow(x[i],j)); A[i][j]=pow(x[i],j); } } double AT[wn][n];//3行6列 for(int i=0;i<wn;i++){ 
        //3 for(int j=0;j<n;j++){ 
        //6 AT[i][j]=pow(x[j],i); } } printf("此时的A：\n"); for(int i=0;i<n;i++){ 
         for(int j=0;j<wn;j++){ 
         printf("%f ",A[i][j]); } printf("\n"); } printf("此时的AT：\n"); for(int i=0;i<wn;i++){ 
         for(int j=0;j<n;j++){ 
         printf("%f ",AT[i][j]); } printf("\n"); } float ATtimesA[wn][wn]; float sum=0; int q=0,p=0; //计算AT*A int k=0; for(int i=0;i<wn;i++){ 
        //AT行遍历 sum=0; for(int j=0;j<n;j++){ 
         sum=sum+AT[i][j]*A[j][k]; } ATtimesA[q][p++]=sum; k++; if(k==wn){ 
         k=0; }else{ 
         i--; } if(p==wn){ 
         p=0; q++; } } for(int i=0;i<wn;i++){ 
         printf("| "); for(int j=0;j<wn;j++){ 
         printf("%f ",ATtimesA[i][j]); } printf("\t|\n"); } float ATtimesY[wn];//AT*Y sum=0; for(int i=0;i<wn;i++){ 
         sum=0; for(int j=0;j<n;j++){ 
         sum=sum+AT[i][j]*y[j]; } ATtimesY[i]=sum; } for(int i=0;i<wn;i++){ 
         printf("%f ",ATtimesY[i]); } printf("\n"); float b[wn+1];//保存三个未知数 float Matrix[wn+1][wn+2]; int o=0; //合并ATtimesA和ATtimesY，得到的矩阵然后使用高斯消元法 for (int i = 1; i <= wn; i++) { 
         for (int j = 1; j <= wn+1; j++) { 
         if(j==wn+1){ 
         Matrix[i][j]=ATtimesY[o++]; }else{ 
         Matrix[i][j]=ATtimesA[i-1][j-1]; } } } printf("输出合并矩阵:\n"); for (int i = 1; i <= wn; i++) { 
         for (int j = 1; j <= wn+1; j++) { 
         printf("%f ",Matrix[i][j]); } printf("\n"); } float m[100][100]; //以下为高斯消元部分 //高斯消元得到最简矩阵 for (int i = 1; i <= wn-1; i++) { 
         for (int j = i+1; j <= wn; j++) { 
         m[j][i] = Matrix[j][i] / Matrix[i][i]; for (int k = i; k <= wn+1; k++) { 
         Matrix[j][k] = Matrix[j][k] - m[j][i] * Matrix[i][k]; } } } printf("输出简化矩阵：\n"); for (int i = 1; i <= wn; i++) { 
         for (int j = 1; j <= wn+1; j++) { 
         printf("%f ",Matrix[i][j]); } printf("\n"); } //输出解 b[wn] = Matrix[wn][wn+1] / Matrix[wn][wn]; for (int i = wn-1; i >0; i--) { 
         for (int j = i+1; j <= wn; j++) { 
         Matrix[i][wn+1] -= Matrix[i][j] * b[j]; } b[i] = Matrix[i][wn+1] / Matrix[i][i]; } for (int i = 1; i <= wn; i++) { 
         printf("%f ",b[i]); } return 0; } 

测试1：

y=4.714284-2.785712X+0.5X^2

测试2：

y=1.099828+2.200162X+3.299967X^2 验算： 当x=1.2时 1.099828+2.6401944+4.75195248=8.49197488≈8.492