non-linearequation.cpp

求解高 阶 线性方程组

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

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		a[5][5]={{0.8,10.2,3.5,-1.8},{8.3,2.1,-1.2,0.5},{1.2,0.2,-4.,-0.5},{-0.2,0.3,0.4,-2.},{1.,1.,1.,1.,1.}};
//	double		b[5]={4.79,-3.02,-6.72,8.89,-2.5};
	double		a[5][5]={{2.,1.,-5,1},{1,-3,0,-6,0},{0.,2.,-1.,2.},{1,4,-7.,6.,0.},{1,1,1,1,1}};
	double		b[5]={8,9,-5,0,-1};
	int			m[5];
	double		x[5];
	ludecomposition(a,m);
	linearsolution(a,m,b,x);
}

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];
	}
}