最小二乘法三次多项式曲线拟合 算法 C++ 实现，该怎么处理

Answer 1

/******************************************
//参考 ã€Šå¸¸ç”¨ç®—法程序集 ï¼ˆC语言描述 ç¬¬ä¸‰ç‰ˆ)》
 
//最小二乘法
//x[n] y[n] å·²çŸ¥è¾“å…¥
//n输入点个数
//a[m] è¿”回m-1次拟合多项式的m个系数
//m  æ‹Ÿåˆå¤šé¡¹å¼çš„项数,即拟合多项式的最高次为m-1
//dt[3] dt[0]返回拟合多项式与各数据点误差的平方和，
        dt[1]返回拟合多项式与各数据点误差的绝对值之和
        dt[2]返回拟合多项式与各数据点误差的绝对值的最大值
 
//
//拟合多项式的输出
//Y(x) = a0 + a1(x-X) + a2(x-X)^2 + â€¦â€¦ am(x-X)^m
// å…¶ä¸­X为已知点x的平均值
******************************************/
 
#include "math.h"
void pir1(x,y,n,a,m,dt)
int n,m;
double x[],y[],a[],dt[];
{ 
  int i,j,k;
  double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
  for (i=0; i<=m-1; i++) a[i]=0.0;
  if (m>n) m=n;
  if (m>20) m=20;
  z=0.0;
  for (i=0; i<=n-1; i++) z=z+x[i]/(1.0*n);
  b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
  for (i=0; i<=n-1; i++)
  { p=p+(x[i]-z); c=c+y[i];}
  c=c/d1; p=p/d1;
  a[0]=c*b[0];
  if (m>1)
  { t[1]=1.0; t[0]=-p;
  d2=0.0; c=0.0; g=0.0;
  for (i=0; i<=n-1; i++)
  { q=x[i]-z-p; d2=d2+q*q;
  c=c+y[i]*q;
  g=g+(x[i]-z)*q*q;
  }
  c=c/d2; p=g/d2; q=d2/d1;
  d1=d2;
  a[1]=c*t[1]; a[0]=c*t[0]+a[0];
  }
  for (j=2; j<=m-1; j++)
  { s[j]=t[j-1];
  s[j-1]=-p*t[j-1]+t[j-2];
  if (j>=3)
    for (k=j-2; k>=1; k--)
      s[k]=-p*t[k]+t[k-1]-q*b[k];
    s[0]=-p*t[0]-q*b[0];
    d2=0.0; c=0.0; g=0.0;
    for (i=0; i<=n-1; i++)
    { q=s[j];
    for (k=j-1; k>=0; k--)
      q=q*(x[i]-z)+s[k];
    d2=d2+q*q; c=c+y[i]*q;
    g=g+(x[i]-z)*q*q;
    }
    c=c/d2; p=g/d2; q=d2/d1;
    d1=d2;
    a[j]=c*s[j]; t[j]=s[j];
    for (k=j-1; k>=0; k--)
    { a[k]=c*s[k]+a[k];
    b[k]=t[k]; t[k]=s[k];
    }
  }
  dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
  for (i=0; i<=n-1; i++)
  { q=a[m-1];
  for (k=m-2; k>=0; k--)
    q=a[k]+q*(x[i]-z);
  p=q-y[i];
  if (fabs(p)>dt[2]) dt[2]=fabs(p);
  dt[0]=dt[0]+p*p;
  dt[1]=dt[1]+fabs(p);
  }
  return;
}