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

最新推荐文章于 2023-07-18 09:56:59 发布

老汤不脑瘫

最新推荐文章于 2023-07-18 09:56:59 发布

阅读量1.6k

点赞数 1

文章标签： c语言最小二乘法矩阵

本文链接：https://blog.csdn.net/qq_52292165/article/details/125156597

版权

多项式拟合高斯消元 C语言矩阵运算数值解

关键词由CSDN通过智能技术生成

在这里插入图片描述

注意：本代码仅针对本公式：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("请输入已知点个数：\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;
}