mintwomult.c

数值计算方法中几个重要的算法用VC++实现

/*********************************************/
/*                                           */
/*                 曲线拟合                  */
/*                                           */
/*********************************************/
#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]);
	}
}