main.cpp

数值分析中实现对给定数表进行插值、曲面拟合、逼近

#include "ChaZhi.h"
#include "BiJin.h"

void main()
{
	int i,j,k,r,s;
	double q[4];
	POINT ftu[6][6],fxy[11][21],pxy[11][21];
	double **B,**G,**C;
	double de;
	B = (double **)malloc(11 * sizeof(double));
	G = (double **)malloc(21 * sizeof(double));
	IniTabFtu(ftu);
	IniTabFxy(fxy);
	IniTabFxy(pxy);
	//计算fxy
	for (i=0; i<11; i++)
	{
		for (j=0; j<21;j++)
		{
			for (k=0; k<4; k++)
			{
				q[k] = 0.8;
			}
			Newton1(q,0.08 * (double)i,0.5 + 0.05 * (double)j);
			fxy[i][j].a3 = PartInsertValue(ftu,q[0],q[1]);			
		}
	}
	//下面输出fxy
	cout<<"1、输出数表xi,yj,f(xi,yj)。（i=0,1,...,10;j=0,1,...,20）"<<endl;
	for (k=0; k<21/3; k++)
	{
		printf("  z y\t");
		for (j=k*3; j<(k+1)*3; j++)
		{
			printf("\t%2.2f\t\t",fxy[0][j].a2);
		}
		cout<<"  x"<<endl;
		for (i=0; i<11; i++)
		{
			printf("  %2.2f\t",fxy[i][0].a1);
			for (j=k*3; j<(k+1)*3; j++)
			{
				printf("%12.11e\t",fxy[i][j].a3);				
			}
		}
	}

	//动态分配矩阵C及每行的空间
	for (k=2; k<10; k++)
	{

		C = (double **)malloc(k * sizeof(double));
		for (i=0; i<k; i++)
		{
			C[i] = (double *)malloc(k * sizeof(double));
		}
		//初始化B,G
	 	IniMatrixBG(fxy,B,G,k);
		//计算C
		ComputeMatrixC(fxy,B,G,C,11,21,k);
		//计算p(x,y)
		for (i=0; i<11; i++)
		{
			for (j=0; j<21; j++)
			{
				pxy[i][j].a3 = 0.0;
				for (r=0; r<k; r++)				
					for (s=0; s<k; s++)					
						pxy[i][j].a3 = pxy[i][j].a3 + C[r][s] * pow(pxy[i][j].a1,r) *
									  pow(pxy[i][j].a2,s);
					//cout<<pxy[i][j].a3<<" ";
			}			
		}
		
		//计算拟和精度
		de = 0.0;
		for (i=0; i<11; i++)
		{
			for (j=0; j<21; j++)
			{
				de = de + (fxy[i][j].a3 - pxy[i][j].a3) * (fxy[i][j].a3 - pxy[i][j].a3); 				
			}
		}
		printf("k=%d时，σ=%12.11e\n",k,de);
		if(de <= 1e-7)
		{
			printf("k=%d时，拟合函数p(x,y)达到精度要求, 此时σ=%12.11e\n",k,de);
			printf("p(x,y)中的系数Crs:\n");
			for (i=0; i<k; i++)
			{
				for (j=0; j<k; j++)
				{
					printf("C%d%d=%12.11e",i,j,C[i][j]);
					if((k*i+j)%3==2)printf("\n");
					else printf("\t");
				}
			}
			printf("\n");
			break;
		}		
	}
	
}