最小二乘法

最新推荐文章于 2024-07-28 08:40:13 发布

勿忘初心，方得始终。

最新推荐文章于 2024-07-28 08:40:13 发布

阅读量214

点赞数

分类专栏：数值分析 C 文章标签：最小二乘法算法

本文链接：https://blog.csdn.net/qq_29144629/article/details/122510893

版权

C 同时被 2 个专栏收录

9 篇文章 0 订阅

订阅专栏

数值分析

3 篇文章 0 订阅

订阅专栏

#include "stdio.h"
#include "math.h"
#include "stdlib.h"
typedef struct 
{
    double *Wi;
    double *G_X;
    double *G_Fx;
    int     G_Num;
    int     S_Order;
    /* data */
} Data_Typedef;


void Newton(void *Matrix,double *B,int Order,double *Result)
{
    // double *X = malloc(sizeof(double)*Order);
    for(int i = 0;i<Order-1;i++)
    {
        for(int j = i+1;j<Order;j++)
        {
            double k = ((double (*)[Order])Matrix)[j][i]/((double (*)[Order])Matrix)[i][i];
            for(int m = i;m<Order;m++)
            {
                ((double (*)[Order])Matrix)[j][m]-= k*((double (*)[Order])Matrix)[i][m];
            }
            B[j] -=k*B[i];
        }
    }
    for(int i = Order-1;i>=0;i--)
    {
        double Tmp = 0;
        for(int j = Order-1;j>i;j--)
        {
            Tmp += Result[j]*((double (*)[Order])Matrix)[i][j];
        }
        Result[i] = (B[i]-Tmp)/((double (*)[Order])Matrix)[i][i];
    }
}

void Least_Squares(Data_Typedef *Para)
{
    double *Matrix = malloc(sizeof(double)*(Para->S_Order+1)*(Para->S_Order+1));
    double *B = malloc(sizeof(double)*(Para->S_Order+1));
    double  *Poly = malloc(sizeof(double)*(Para->S_Order+1));
    for(int i = 0;i<Para->S_Order+1;i++)
    {
        for(int j = 0;j<Para->S_Order+1;j++)
        {
            double tmpA = 0;
            for(int m = 0;m<Para->G_Num;m++)
            {
                tmpA +=(pow(Para->G_X[m],(i+j))*Para->Wi[m]);
            }   
            ((double (*)[Para->S_Order+1])Matrix)[i][j] = tmpA;
            
        }

        double tmpB = 0;
        for(int m = 0;m<Para->G_Num;m++)
        {
            tmpB += pow(Para->G_X[m],i)*Para->G_Fx[m]*Para->Wi[m];
        }  
        B[i] = tmpB;
    }
    Newton(Matrix,B,Para->S_Order+1,Poly);
    printf("G(x) = ");
    for(int i = Para->S_Order+1;i>0;i--)
    {
        if(i != (Para->S_Order+1))
        {
            if(Poly[i-1]>=0)
            {
                printf("+");
            }
        }
        printf("%.10lf",Poly[i-1]);
        if(i-1>0)
        {
            printf("X^%d",i-1);
        }
    }
    printf("\n");
    for(int i = 0;i<Para->G_Num;i++)
    {
        // for(int j = 0;j < Para->G_Num;j++)
        // {}
        double Sum[4] = {0};
        double Al = 0;
        for(int m = 0;m < Para->S_Order+1;m++)
        {
            Sum[m] += Poly[m]*pow(Para->G_X[i],m);
            Al += Sum[m];
        }
        
        printf("G(%.0f) \t= %.9lf\t",Para->G_X[i],Al);
        printf("e(%.0f) \t= %.9lf\t",Para->G_X[i],Al-Para->G_Fx[i]);
        // if(i>0&&i%4==0)
        {
            printf("\n");
        }
    }
    free(Matrix);
    free(B);
    free(Poly);
}


int main(int arg)
{
    Data_Typedef My_D;
    double G_X[] =  {0,5,10,15,20,25,30,35,40,45,50,55};
    double G_Fx[] = {0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64};
    double Wi[]   = {1,1,1,1,1,1,1,1,1,1,1,1,};
    int     Order = 2;
    for(int i = 0;i<sizeof(G_Fx)/sizeof(double);i++)
    {
        G_Fx[i]/=10000;
    }

    My_D.G_Fx   = G_Fx;
    My_D.G_X    = G_X;
    My_D.Wi     = Wi;
    My_D.S_Order  = 3;
    My_D.G_Num  = sizeof(G_X)/sizeof(double);
    Least_Squares(&My_D);
    
}