/*********************************************/
/* */
/* 曲线拟合 */
/* */
/*********************************************/
#include<math.h>
#include<stdio.h>
#include<malloc.h>
#define type "%lf" /*-定义数据类型-*/
typedef double Dtype;
Dtype max(Dtype **a,int k,int n,int *ik,int *jk)
{
/*-用*ik和*jk返回最大值所在的行号和列号-*/
int i,j;
Dtype maxValue;
for(i=k;i<=n;i++)
for(j=k;j<=n;j++)
if(i==k&&j==k)
{
maxValue=a[i][j];
*ik=i;*jk=j;
}
else if(fabs(a[i][j])>fabs(maxValue))
{
maxValue=a[i][j];
*ik=i;*jk=j;
}
return maxValue;
}/*-max-*/
int xmap(int origin[],int n,int row)
{
/*-返回与xrow的映射-*/
int i;
for(i=0;i<=n;i++)
if(origin[i]==row) return i;
}/*-xmap-*/
void main()
{
int i,j,k,m,n,ik,jk,tcol;
int *origin;/*-原始未知数的顺序-*/
Dtype *x,*y,**A,**C,*b,m1,sum,temp,*a;
printf("Please input m,n(m>0,n>0):\n");
scanf("%d%d",&m,&n);
/*动态的分配空间*/
origin=(int*)malloc(sizeof(int)*(n+1));
x=(Dtype*)malloc(sizeof(Dtype)*m);
y=(Dtype*)malloc(sizeof(Dtype)*m);
b=(Dtype*)malloc(sizeof(Dtype)*(n+1));
a=(Dtype*)malloc(sizeof(Dtype)*(n+1));
A=(Dtype**)malloc(sizeof(Dtype*)*(n+1));
C=(Dtype**)malloc(sizeof(Dtype*)*m);
for(i=0;i<n+1;i++)
A[i]=(Dtype*)malloc(sizeof(Dtype)*(n+1));
for(j=0;j<m;j++)
C[j]=(Dtype*)malloc(sizeof(Dtype)*(n+1));
/*输入x的值*/
printf("Please input array x[1..%d]:\n",m);
for(i=0;i<m;i++)
scanf(type,&x[i]);
/*输入y的值*/
printf("Please input array y[1..%d];\n",m);
for(i=0;i<m;i++)
scanf(type,&y[i]);
/*生成矩阵C*/
for(i=0;i<m;i++)
{ C[i][0]=1;
for(j=1;j<n+1;j++)
C[i][j]=pow(x[i],(double)j);
}
/*生成矩阵A=CT*C*/
for(i=0;i<n+1;i++)
for(j=0;j<n+1;j++)
{
A[i][j]=0;
for(k=0;k<m;k++)
A[i][j]+=C[k][i]*C[k][j];
}
/*生成向量b=CT*Y*/
for(i=0;i<n+1;i++)
{
b[i]=0;
for(k=0;k<m;k++)
b[i]+=C[k][i]*y[k];
}
/*利用完全主元素削去法求a[0..n]*/
/*-对原始顺序进行记录-*/
for(i=0;i<=n;i++)
origin[i]=i;
for(k=0;k<=n-1;k++)
{
max(A,k,n,&ik,&jk);
if(ik!=k) /*-交换行ik<-->k-*/
{
for(j=0;j<=n;j++)
{
temp=A[ik][j];
A[ik][j]=A[k][j];
A[k][j]=temp;
}
/*-交换b[ik]<-->b[k]-*/
temp=b[ik];b[ik]=b[k];b[k]=temp;
}
if(jk!=k) /*-交换列jk<--->k-*/
{
for(i=0;i<=n;i++)
{
temp=A[i][jk];
A[i][jk]=A[i][k];
A[i][k]=temp;
}
tcol=origin[jk]; /*-记录未知数的新位置-*/
origin[jk]=origin[k];
origin[k]=tcol;
}
for(i=k+1;i<=n;i++) /*-消元计算-*/
{
m1=A[i][k]/A[k][k];
for(j=k+1;j<=n;j++)
A[i][j]-=m1*A[k][j];
b[i]-=m1*b[k];
}
}
/*-进行回代-*/
a[n]=b[n]/A[n][n];
for(i=n-1;i>=0;i--)
{
sum=0;
for(j=i+1;j<=n;j++)
sum+=A[i][j]*a[j];
a[i]=(b[i]-sum)/A[i][i];
}
/*-打印输出x1,x2,x3......xn-*/
for(i=0;i<=n;i++)
{
j=xmap(origin,n,i);/*-求映射xi<-->xj-*/
printf("x%d=%.8lf\n",i,a[j]);
}
}
C程序.rar （11个子文件） 
计算方法实验C程序 实验 线性方程组的直接法 LU 
LU.cpp 2KB FullMainElemDele FullMainElemDele.cpp 2KB Jacobi Jacobi.cpp 2KB Chase Chase.cpp 1KB Square Square.cpp 2KB 插值方法 Newton Newton.C 1KB Lagrange Lagrange.C 1012B 曲线拟合 MinTwoMult MinTwoMult.C 3KB 常微分方程的数值解 R_K R_K.cpp 825B 数值积分与数值微分 Newton_Cotes Simposn.C 997B Romberg Romberg.cpp 1KB资源评论
