#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <math.h>
#define p(k) (p+(k)*m)
#define pkx(k) (pkx+(k)*N)
double *x,*f,*p,*a,*b,*d,*pkx,*sx,*xishu;
// x 存放节点的值 f存放节点对应的函数值 p存放各多项式在各节点的值
// a b 分别存放αβ的值 d 存放a*的值 pkx存放各多项式Pn(x)的系数
// sx 存放拟合函数在各节点的函数值 xishu 存放拟合多项式的系数
int m,n,N;
void GetDate() //分配存储空间和输入拟合数据
{
int i;
i=sizeof(double);
// 分配存储空间
x=(double *)malloc(i*m);
f=(double *)malloc(i*m);
p=(double *)malloc(i*m*n);
a=(double *)malloc(i*n);
b=(double *)malloc(i*n);
d=(double *)malloc(i*n);
pkx=(double *)malloc(i*n*n);
sx=(double *)malloc(i*m);
xishu=(double *)malloc(i*n);
if(!(x&&f&&p&&a&&b&&d&&pkx&&sx&&xishu))
{
printf("分配内存失败! 请按任意建结束程序\n");
getch();
exit(0);
}
fflush(stdin); // 清除缓冲区的数据
printf("请输入节点的值 :\n");
for(i=0;i<m;i++)
scanf("%lf",&x[i]);
fflush(stdin); // 清除缓冲区的数据
printf("请输入前面节点对应的函数值 :\n");
for(i=0;i<m;i++)
scanf("%lf",&f[i]);
} // GetDate end
double NeiJi(double *F,double *G,double *W) //内积函数
{
double s=0;
int i;
for (i=0;i<m;i++)
s+=F[i]*G[i]*W[i];
return s;
} // NeiJi end
void GetPk(int k) // 求正交多项式Pk(x)在节点 x(k) 的值
{
int i;
for (i=0;i<m;i++)
p(k)[i]=(x[i]-a[k])*p(k-1)[i]-b[k-1]*p(k-2)[i];
}
void GetXiShu(int k) // 求拟合多项式的系数
{
int i,j;
double t;
for (i=0;i<k;i++)
{
t=0;
for (j=0;j<k;j++)
{
t+=d[j]*pkx(j)[i];
}
xishu[i]=t;
}
} // GetXiShu end
double Wucha(int k) // 计算平方误差
{
double r1;
int i,j;
for (i=0;i<m;i++)
{
sx[i]=0;
for (j=0;j<k;j++)
sx[i]=sx[i]+xishu[j]*pow(x[i],j); // 计算拟合函数在各节点的值
}
r1=fabs(NeiJi(f,f,p(0))-NeiJi(f,sx,p(0))); // 计算平方误差
return r1;
} // Wucha end
void main()
{
int k,i,j;
double t,r=0,R;
printf("\n本程序是用正交函数作最小二乘拟合\n");
printf("请输入节点个数:");
scanf("%d",&m);
fflush(stdin); // 清除缓冲区的数据
printf("请选择拟合控制条件\n");
printf("按 n 根据正交多项式最高次数控制\n按 其它键 根据拟合的平方误差控制: ");
if (getchar()=='n')
{
printf("\n请输入 n 的值(0-%d): ",m-1);
fflush(stdin); // 清除缓冲区的数据
scanf("%d",&n);
n=min(n+1,m);
N=n;
}
else
{
n=m;
N=n;
printf("\n请输入平方误差: ");
fflush(stdin); // 清除缓冲区的数据
scanf("%lf",&r);
}
GetDate(); // 输入数据
for (i=0;i<n;i++)
{
for (j=0;j<n;j++)
pkx(i)[j]=0; // pkx 初始化
}
for (k=0;k<m;k++)
{
p(0)[k]=1; // P0(x)=1;
}
pkx(0)[0]=1;
d[0]=NeiJi(f,p(0),p(0))/NeiJi(p(0),p(0),p(0));
GetXiShu(1);
R=Wucha(1);
if(R<=r) n=1;
if (n>1&&R>r)
{
a[1]=NeiJi(x,p(0),p(0))/NeiJi(p(0),p(0),p(0)); //α1=(xP1,P1)/(P1,P1)
pkx(1)[0]=-a[1];pkx(1)[1]=1;
for (k=0;k<m;k++)
{
p(1)[k]=x[k]-a[1]; // P1(x)=(x-α1)P0(x) P0(x)=1
}
d[1]=NeiJi(f,p(1),p(0))/NeiJi(p(1),p(1),p(0));
GetXiShu(2);
R=Wucha(2);
if(R<=r) n=2;
}
for (k=1;k<n-1&&R>r;k++)
{
t=NeiJi(p(k),p(k),p(0));
a[k+1]=NeiJi(p(k),p(k),x)/t; // αk+1=(xPk,Pk)/(Pk,Pk)
b[k]=t/NeiJi(p(k-1),p(k-1),p(0)); //βk=(Pk,Pk)/(Pk-1,Pk-1)
GetPk(k+1); // 计算Pk(x)在各节点的函数值
d[k+1]=NeiJi(f,p(k+1),p(0))/NeiJi(p(k+1),p(k+1),p(0));
// ak=(f,Pk)/(Pk,Pk)
pkx(k+1)[0]=-a[k+1]*pkx(k)[0]-b[k]*pkx(k-1)[0];
for (i=1;i<n;i++)
{ // Pk+1=(x-αk+1)Pk(x)-βkPk-1(x)
pkx(k+1)[i]=pkx(k)[i-1]-a[k+1]*pkx(k)[i]-b[k]*pkx(k-1)[i];
}
GetXiShu(k+1); // 求拟合多项式的系数
R=Wucha(k+1); //计算平方误差
if(R<=r)
{
n=k+1;
break;
}
}
if (R>r)
{
GetXiShu(n); // 求拟合多项式的系数
R=Wucha(n); //计算平方误差
}
//////////////下面是输出结果/////////////////////
printf("\n正交多项式Pk(x)为: \n");
for (k=0;k<n;k++)
{
printf("P%d(x)=%lG",k,pkx(k)[0]);
for (i=1;i<k+1;i++)
if(pkx(k)[i]>=0)
printf("+%lGx^%d",pkx(k)[i],i);
else
printf("%lGx^%d",pkx(k)[i],i);
printf("\n");
}
printf("\n平方误差为: %lf\n",R);
printf("\n拟合多项式为: \n");
printf("F(x)=%lG",xishu[0]);
for (k=1;k<n;k++)
{
if(xishu[k]>=0)
printf("+%lGx^%d",xishu[k],k);
else
printf("%lGx^%d",xishu[k],k);
}
printf("\n完成! 按任意键结束程序\n");
getch();
free(x);free(f);free(a);free(p);free(b);
free(d);free(pkx);free(sx);free(xishu);
}
