newton1.cpp

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

#include "Newton1.h"

//初始化f(t,u)
void IniTabFtu(POINT ftu[6][6])
{
	int i,j;
	for (i=0; i<6; i++)
	{
		for (j=0; j<6; j++)
		{
			ftu[i][j].a1 = 0.2 * i;
            ftu[i][j].a2 = 0.4 * j;
		}
	}
	ftu[0][0].a3 = -0.5;  ftu[0][1].a3 = -0.34; ftu[0][2].a3 = 0.14;
	ftu[0][3].a3 = 0.94;  ftu[0][4].a3 = 2.06;  ftu[0][5].a3 = 3.5;
	ftu[1][0].a3 = -0.42; ftu[1][1].a3 = -0.5;  ftu[1][2].a3 = -0.26;
	ftu[1][3].a3 = 0.3;	  ftu[1][4].a3 = 1.18;  ftu[1][5].a3 = 2.38;
	ftu[2][0].a3 = -0.18; ftu[2][1].a3 = -0.5;  ftu[2][2].a3 = -0.5;
	ftu[2][3].a3 = -0.18; ftu[2][4].a3 = 0.46;  ftu[2][5].a3 = 1.42;
	ftu[3][0].a3 = 0.22;  ftu[3][1].a3 = -0.34; ftu[3][2].a3 = -0.58;
	ftu[3][3].a3 = -0.5;  ftu[3][4].a3 = -0.1;  ftu[3][5].a3 = 0.62;
	ftu[4][0].a3 = 0.78;  ftu[4][1].a3 = -0.02; ftu[4][2].a3 = -0.5;
	ftu[4][3].a3 = -0.66; ftu[4][4].a3 = -0.5;  ftu[4][5].a3 = -0.02;
	ftu[5][0].a3 = 1.5;   ftu[5][1].a3 = 0.46;  ftu[5][2].a3 = -0.26;
	ftu[5][3].a3 = -0.66; ftu[5][4].a3 = -0.74; ftu[5][5].a3 = -0.5;
}

//计算关于t,u,v,w的方程组的每个函数-fi的值
//其中q[0]、q[1]、q[2]和q[3]分别代表t,u,v和w
void Function(double f[4], double q[4], double x, double y)
{
	f[0] = -1.0 * (0.5 * cos(q[0]) + q[1] + q[2] + q[3] - x - 2.67);
	f[1] = -1.0 * (q[0] + 0.5 * sin(q[1]) + q[2] + q[3] - y - 1.07);
	f[2] = -1.0 * (0.5 * q[0] + q[1] + cos(q[2]) + q[3] - x - 3.74);
	f[3] = -1.0 * (q[0] + 0.5 * q[1] + q[2] + sin(q[3]) - y - 0.79);
}

//求Jacobi矩阵
//其中q[0]、q[1]、q[2]和q[3]分别代表t,u,v和w
void Jacobi(double **J,double q[4])
{
	J[0][0] = -0.5 * sin(q[0]); J[0][1] = 1.0; J[0][2] = 1.0; J[0][3] = 1.0;
	J[1][0] = 1.0; J[1][1] = 0.5 * cos(q[1]); J[1][2] = 1.0; J[1][3] = 1.0;
	J[2][0] = 0.5; J[2][1] = 1.0; J[2][2] = -1.0 * sin(q[2]); J[2][3] = 1.0;
	J[3][0] = 1.0; J[3][1] = 0.5; J[3][2] = 1.0; J[3][3] = cos(q[3]);
}

//对矩阵A作选主元的LU分解操作
void LU(double **a, int *M, int n)
{
	//先作分解QA=LU
	//将L和U的值存入a中
	double *s = (double*)malloc(n * sizeof(double));
	int i,j,k,t;
	double temp;
	for (k=0; k<n; k++)
	{
		double s_i_k = 0.0;
		for (i=k; i<n; i++)
		{
			s[i] = a[i][k];
			for (t=0; t<=k-1; t++)
			{
				s[i] = s[i] - a[i][t] * a[t][k]; 
			}
			if (fabs(s[i]) > s_i_k) 
			{
				s_i_k = fabs(s[i]);
				M[k] = i;
			}
			
		}
		if (M[k] != k)
		{
			for (t=0; t<n; t++)
			{
				temp = a[k][t];
				a[k][t] = a[M[k]][t];
				a[M[k]][t] = temp;
			}
			temp = s[k];
			s[k] = s[M[k]];
			s[M[k]] = temp;
		}
		a[k][k] = s[k];
		for (j=k+1; j<n; j++)
		{
			for (t=0; t<=k-1; t++)
			{
				a[k][j] = a[k][j] - a[k][t] * a[t][j];
			}	
		}
		for (i=k+1; i<n; i++)
		{
			a[i][k] = s[i] / a[k][k];
		}		
	}


}
//用Doolittle分解法求解线性方程组
void Doolittle(double **a,double *b,int *M,double *x,int m)
{
	int i,k,t;
	double temp;
	//求Qb
	for (k=0; k<m-1; k++)
	{
		if (M[k] != k)
		{
			temp = b[k];
			b[k] = b[M[k]];
			b[M[k]] = temp;
		}
		
	}
	//下面求解Ly=Qb和Ux=y
	//将y的值存入b中
	for (i=1; i<m; i++)
	{
		for (t=0; t<=i-1; t++)
		{
			b[i] = b[i] - a[i][t] * b[t];
		}
	}
	x[m-1] = b[m-1] /a[m-1][m-1];
	for (i=m-2; i>=0; i--)
	{
		x[i] = b[i];
		for (t=i+1; t<=3; t++)
		{
			x[i] = x[i] - a[i][t] * x[t];
		}
		x[i] = x[i] / a[i][i];
	}
		
}

//求向量的无穷范数
double InfiFS(double x[4])
{
	int i;
	double FS = 0.0;
	for (i=0; i<4; i++)
	{
		if (fabs(x[i]) > fabs(FS))
		FS = fabs(x[i]);
	}
	return FS;
}

//用Newton法解非线性方程组
void Newton1(double q[4],double x,double y)
{
	int k,i;
	int Max = 1024;
	double eps = 1e-12;
	double f[4] = {0};
	double dx[4] = {0};
    double **J;
	J= (double**)malloc(4 * sizeof(double));
	int *M = (int*)malloc(4 * sizeof(int));
	for (i=0; i<4; i++)
	{
		J[i] = (double*)malloc(4 * sizeof(double));
	}
	for (k=0; k<Max; k++)
	{
		// cout<<"x:"<<x<<" y:"<<y<<endl;
		Function(f, q, x, y); 
//		for (i=0; i<4; i++)
//		{
//			cout<<f[i]<<endl;
//		}
		// cout<<f[0]<<endl;
        Jacobi(J,q);
//		for (i=0;i<4;i++)
//		{
//			for (j=0;j<4;j++)
//			{
//				cout<<J[i][j]<<" ";
//			}
//            cout<<endl;
//		}
		LU(J,M,4);
		Doolittle(J,f,M,dx,4);

//		for (i=0; i<4; i++)
//		{
//			cout<<dx[i]<<endl;
//		}
		if (InfiFS(dx) / InfiFS(q) < eps)
		{
		// 	cout<<"迭代成功！"<<endl;
		// 	cout<<q[0]<<" ";
			return;
		}
		
		for (i=0; i<4; i++)
		{
			q[i] = q[i] + dx[i];
		}
		
	}
	cout<<"迭代失败！"<<endl;
	return;
}