计算方法实验六 最小二乘法

#include<math.h>
#include<stdio.h>
#include<string.h>
#include<algorithm>
using namespace std;

double num1[]={0,5,10,15,20,25,30,35,40,45,50,55};
double num2[]={0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64};

#define MAXN 1010
#define number 0.0001
#define mem(a) memset(a,0,sizeof(a))
const double eps=1e-5;

double temp[MAXN];
int len=3,index=12;

double arr[MAXN][MAXN];
double ans[MAXN];
int equ=len+1,var=len+1;///方程的个数，变量的个数

void deal()
{
    mem(temp),mem(ans);
    //处理第0个
    temp[0]=index;
    for(int i=0;i<index;i++)
        ans[0]+=num2[i]*number;

    double num[MAXN];
    for(int i=0;i<index;i++)
        num[i]=1;

    for(int i=1;i<=len*2;i++){
        for(int j=0;j<index;j++){
            num[j]*=num1[j];
            temp[i]+=num[j];
            ans[i]+=num[j]*num2[j]*number;
        }
    }

    //转换为矩阵
    for(int i=0;i<=len;i++)
        for(int j=0;j<=len;j++)
            arr[i][j]=temp[i+j];
}

int Gauss()
{
    int row=0,col=0;
    for(;row<equ&&col<var;row++,col++)
    {
        int max_r=row;
        for(int i=row+1;i<equ;i++)
            if(fabs(arr[i][col])>fabs(arr[max_r][col]))
                max_r=i;

        if(fabs(arr[max_r][col])<eps)
            return 0;

        if(row!=max_r){
            for(int j=col;j<var;j++)
                swap(arr[row][j],arr[max_r][j]);
            swap(ans[row],ans[max_r]);
        }

        ans[row]/=arr[row][col];
        for(int j=col+1;j<var;j++)
            arr[row][j]/=arr[row][col];
        arr[row][col]=1;

        for(int i=0;i<equ;i++){
            if(i!=row){
                ans[i]-=ans[row]*arr[i][row];
                for(int j=col+1;j<var;j++)
                    arr[i][j]-=arr[row][j]*arr[i][col];
                arr[i][col]=0;
            }
        }
    }
    return 1;
}

double calc(int t)
{
    double temp=0;
    for(int i=0;i<=len;i++)
        temp+=ans[i]*pow(t,i);

    return temp;
}

int main()
{
    deal();
    Gauss();
    int flag=0;
    printf("这个多项式就是:");
    for(int i=0;i<=len;i++){
        if(ans[i]){
            printf(" %c %e",(flag==0?'\n':'+'),ans[i]);
            if(i)
                printf("X^%d",i);
            flag=1;
        }
    }

    int t=0;
    printf("\n误差  ：\n");
    for(int i=0;i<index;i++){
        printf("当t==%d时候 误差：%lf \n",t,num2[i]-calc(t));
        t+=5;
    }
    return 0;
}