正交多项式拟合函数

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

实验要求：

给定函数f(x)在 m 个采样点处的值f()以及每个点的权重（1<= i <= m），求曲线拟合的正交多项式，满足最小二乘误差。

函数接口定义：

int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps )

在接口定义中：f是给定函数f(x)，m为采样点个数，x传入m个采样点，w传入每个点的权重；c传回求得的曲线拟合的正交多项式的系数，eps传入最小二乘误差上限TOL，并且传回拟合的实际最小二乘误差。

另外，在裁判程序中还定义了常数MAXN，当拟合多项式的阶数达到MAXN却仍未达到精度要求时，即返回 n = MAXN时的结果。

#include<math.h>
 
#define MAXM 200 /* 采样点个数上限*/
#define MAXN 5 /* 拟合多项式阶数上限*/
 
/* 这里仅给出2个测试用例，正式裁判程序可能包含更多测试用例 */
double qi_x[MAXN+1][MAXN+1];//存储基底函数
double ak[MAXN+1];//存储拟合函数对应基地的系数
 
 
double f1(double x)
{
  return sin(x);
}
 
double f2(double x)
{
  return exp(x);
}
 
 
//用于求x的n次方
double power(double x,int n)
{
    int i;
    double y=1;
    for(i=0;i<n;i++)
    {
        y=y*x;
    }
    return y;
}
 
//用于求解两个基底函数的内积
double get_Inner_product(double x[],double y1[],double y2[] ,double w[],int m)
{
    int i,j;
    double sum=0;
    double sum1,sum2;
    for(i=0;i<m;i++)
    {
        sum1=0;
        sum2=0;
        for(j=0;j<=MAXN;j++)
        {
            sum1=sum1+y1[j]*power(x[i],j);
            sum2=sum2+y2[j]*power(x[i],j);
        }
        sum=sum+w[i]*sum1*sum2;
    }
    return sum;
}
 
//用于完成x*a(x)的功能
 void carry_x(double x[],double *temp)
{
    int i;
 
    for(i=MAXN+1;i>0;i--)
    temp[i]=x[i-1];
 
    temp[0]=0;
}
 
//用于完成拟合函数和基底函数的内积
double get_function_inner_product(double (*f)(double t),double x[],double y1[],double w[],int m)
{
    int i,j;
    double sum=0;
    double sum1;
    for(i=0;i<m;i++)
    {
      sum1=0;
      for(j=0;j<=MAXN;j++)
      {
          sum1=sum1+y1[j]*power(x[i],j);
      }
      sum=sum+w[i]*f(x[i])*sum1;
    }
   // printf("%f\n",sum);
    return sum;
}
 
//求解正交基底
void Array_op(int n,double fg,double fg1)
{
    int i;
    carry_x(qi_x[n-1],qi_x[n]);
       if(n==1)
       {
           for(i=0;i<=n;i++)
           qi_x[n][i]=qi_x[n][i]-fg*qi_x[n-1][i];
       }
       else{
       for(i=0;i<=n;i++)
        {
            qi_x[n][i]=qi_x[n][i]-fg*qi_x[n-1][i]-fg1*qi_x[n-2][i];
        }
        //qi_x[n][0]=-fg*qi_x[n-1][0]-fg1*qi_x[n-2][0];
       }
}
 
//完成函数拟合并返回误差和拟合次数
int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps )
{
    double err;
    double sum=0;
    double m1,n;
    double temp[MAXN+1]={0};
    int i;
 
    qi_x[0][0]=1;
    ak[0]=get_function_inner_product(f,x,qi_x[0],w,m)/get_Inner_product(x,qi_x[0],qi_x[0],w,m);
 
    for(i=0;i<m;i++)
    {
        sum=sum+w[i]*f(x[i])*f(x[i]);
    }
 
    err=sum-ak[0]*get_function_inner_product(f,x,qi_x[0],w,m);
   // printf("%f\n",err);
 
    carry_x(qi_x[0],temp);
 
    m1=get_Inner_product(x,temp,qi_x[0],w,m)/get_Inner_product(x,qi_x[0],qi_x[0],w,m);
 
    Array_op(1,m1,0.0);
 
    ak[1]=get_function_inner_product(f,x,qi_x[1],w,m)/get_Inner_product(x,qi_x[1],qi_x[1],w,m);
    err=err-ak[1]*get_function_inner_product(f,x,qi_x[1],w,m);
    int k=1;
 
    while((k<MAXN)&&(fabs(err)>=*eps))
    {
        k=k+1;
        double temp1[MAXN+1]={0};
 
        carry_x(qi_x[k-1],temp1);
 
        m1=get_Inner_product(x,temp1,qi_x[k-1],w,m)/get_Inner_product(x,qi_x[k-1],qi_x[k-1],w,m);
 
        n=get_Inner_product(x,qi_x[k-1],qi_x[k-1],w,m)/get_Inner_product(x,qi_x[k-2],qi_x[k-2],w,m);
 
        Array_op(k,m1,n);
       // printf("%f\n",get_function_inner_product(f,x,qi_x[k],w,m));
 
        ak[k]=get_function_inner_product(f,x,qi_x[k],w,m)/get_Inner_product(x,qi_x[k],qi_x[k],w,m);
       // printf("%f\n",ak[k]);
        err=err-ak[k]*get_function_inner_product(f,x,qi_x[k],w,m);
    }
    *eps=err;
    return k;
}
 
//输出实验结果
void print_results( int n, double c[], double eps)
{ /* 用于输出结果 */
  int i;
  int m,k;
  for(m=0;m<MAXN+1;m++)
  {
      c[m]=0;
      for(k=0;k<MAXN+1;k++)
      c[m]=c[m]+ak[k]*qi_x[k][m];
  }
  printf("%d\n", n);
  for (i=0; i<=n; i++)
    printf("%8.4e ", c[i]);
  printf("\n");
  printf("error = %6.2e\n", eps);
  printf("\n");
}
 
int main()
{
  int m, i, n;
  double x[MAXM], w[MAXM], c[MAXN+1], eps;
 
  /* 测试1 */
  m = 90;
  for (i=0; i<m; i++) {
    x[i] = 3.1415926535897932 * (double)(i+1) / 180.0;
    w[i] = 1.0;
  }
  eps = 0.001;
  n = OPA(f1, m, x, w, c, &eps);
  print_results(n, c, eps);
 
  /* 测试2 */
  memset(qi_x,0,sizeof(qi_x));
  memset(ak,0,sizeof(ak));
  m = 200;
  for (i=0; i<m; i++) {
    x[i] = 0.01*(double)i;
    w[i] = 1.0;
  }
  eps = 0.001;
  n = OPA(f2, m, x, w, c, &eps);
  print_results(n, c, eps);
 
  return 0;
}