non-linearequation.cpp

牛顿法求解非线性方程组

#include <stdio.h>
#include <math.h>

void product(double z[5],double x,double y,double b[5]);
void differential(double diffa[5][5],double z[5],double x,double y);
double infiniteparadigm(double deltz[5]);
void ludecomposition(double a[5][5],int m[5]);
void linearsolution(double a[5][5],int m[5],double b[5],double x[5]);

void main()
{
	double z[5]={2.,2.,2.,2.,2.};
	double deltz[5];
	double diffa[5][5];
	double x=1.;
	double y=1.;
	double b[5];
	int    m[5];
	int    k,i;
	for(k=1;k<=1000;k++)
	{
		product(z,x,y,b);
		differential(diffa,z,x,y);
		ludecomposition(diffa,m);
		linearsolution(diffa,m,b,deltz);
		for(i=0;i<=4;i++)
		{
			z[i]+=deltz[i];
		}
		if(infiniteparadigm(deltz)<=1.e-12) break;
	}
	product(z,x,y,b);
}

void product(double z[5],double x,double y,double b[5])
{
	b[0]=0.-(exp(0.-z[0])+exp(-2.*z[1])+z[2]-2.*z[3]+x*z[4]-5.3);
	b[1]=0.-(exp(-2.*z[0])+exp(0.-z[1])-2.*z[2]+y*z[3]-z[4]+25.6);
	b[2]=0.-(x*z[0]+3.*z[1]+exp(0.-z[2])-3.*z[4]+37.8);
	b[3]=0.-(2.*z[0]+y*z[1]+z[2]-exp(0.-z[3])+2.*exp(-2.*z[4])-31.3);
	b[4]=0.-(z[0]-2.*z[1]-3.*x*z[2]+exp(-2.*z[3])+3.*exp(0.-z[4])+42.1);
}
void differential(double diffa[5][5],double z[5],double x,double y)
{
	diffa[0][0]=0.-exp(0.-z[0]);
	diffa[0][1]=-2.*exp(-2.*z[1]);
	diffa[0][2]=1.;
	diffa[0][3]=-2.;
	diffa[0][4]=x;

	diffa[1][0]=-2.*exp(-2.*z[0]);
	diffa[1][1]=0.-exp(0.-z[1]);
	diffa[1][2]=-2.;
	diffa[1][3]=y;
	diffa[1][4]=-1;

	diffa[2][0]=x;
	diffa[2][1]=3;
	diffa[2][2]=0.-exp(0.-z[2]);
	diffa[2][3]=0.;
	diffa[2][4]=-3;

	diffa[3][0]=2.;
	diffa[3][1]=y;
	diffa[3][2]=1;
	diffa[3][3]=exp(0.-z[3]);
	diffa[3][4]=-4.*exp(-2.*z[4]);

	diffa[4][0]=1.;
	diffa[4][1]=-2.;
	diffa[4][2]=-3*x;
	diffa[4][3]=-2.*exp(-2*z[3]);
	diffa[4][4]=-3*exp(0.-z[4]);

}
double infiniteparadigm(double deltz[5])
{
	double result;
	int i;
	result=fabs(deltz[0]);
	for(i=1;i<=4;i++)
	{
		if(result<fabs(deltz[i]))
		{
			result=fabs(deltz[i]);
		}
	}
	return result;
}
void ludecomposition(double a[5][5],int m[5])
{
	double	temp[5];
	double	s[5];
	double	modulemax;
	double	exchange;
	int		k,i,j,t;

	for(k=0;k<=4;k++)
	{
		for(i=k;i<=4;i++)
		{
			s[i]=a[i][k];
			for(t=0;t<=k-1;t++)
			{
				s[i]-=a[i][t]*a[t][k];
			}
		}
		modulemax=fabs(s[k]);
		m[k]=k;
		for(i=k+1;i<=4;i++)
		{
			if(modulemax<s[i])
			{
				m[k]=i;
				modulemax=s[i];
			}
		}
		if(m[k]!=k)
		{
			exchange=s[k];
			s[k]=s[m[k]];
			s[m[k]]=exchange;
			for(t=0;t<=4;t++)
			{
				temp[t]=a[k][t];
			}
			for(t=0;t<=4;t++)
			{
				a[k][t]=a[m[k]][t];
			}
			for(t=0;t<=4;t++)
			{
				a[m[k]][t]=temp[t];
			}
		}
		a[k][k]=s[k];
		for(j=k+1;j<=4;j++)
		{
			for(t=0;t<=k-1;t++)
			{
				a[k][j]-=a[k][t]*a[t][j];
			}
		}
		for(i=k+1;i<=4;i++)
		{
			a[i][k]=s[i]/a[k][k];
		}
	}
}

void linearsolution(double a[5][5],int m[5],double b[5],double x[5])
{
	int		k,i,t;
	double	exchange;
	double	y[5];
	for(k=0;k<=3;k++)
	{
		if(m[k]>k)
		{
			exchange=b[k];
			b[k]=b[m[k]];
			b[m[k]]=exchange;
		}
	}
	y[0]=b[0];
	for(i=1;i<=4;i++)
	{
		y[i]=b[i];
		for(t=0;t<=i-1;t++)
		{
			y[i]-=a[i][t]*y[t];
		}
	}
	x[4]=y[4]/a[4][4];
	for(i=3;i>=0;i--)
	{
		x[i]=y[i];
		for(t=i+1;t<=4;t++)
		{
			x[i]-=a[i][t]*x[t];
		}
		x[i]/=a[i][i];
	}
}