计算方法.cpp

本实验主要涉及解线性方程组的列主元消去法.通过编辑程序,能更深入地了解到列主元消去法的内涵.并且通过计算比较,也可以了解到数据的精度及舍入误差等一系列有关因素.

#include<stdio.h>
#include<math.h>
#define N 3/*方程阶数*/ 

void print(double a[][N],double *b); 
void gaosi(double a[][N],double *b); 
void linechange(double a[][N],double *b,int k);
void huidai(double a[][N],double *b);

double c[N][N] = {{0.2,0.1,-0.8},{0.3,-0.9,0.2},{-1,0.4,0.3}};

void main()
{
	double a[N][N] = {{0.2,0.1,-0.8},{0.3,-0.9,0.2},{-1,0.4,0.3}};
	double b[N] = {-50000,-30000,-60000};

	gaosi(a,b);
	print(a,b);
	huidai(a,b);

	printf( "\nAfter +5000\n\n" );
	for( int i = 0; i <= 2; i++ )
		b[i] -= 5000;
	huidai( a, b );
}

void print(double a[][N],double *b)/*打印运算矩阵a,b*/ 
{
	int i,j;
	printf("a is:\n");
	for(i = 0;i < N;i++)
	{
		for(j = 0;j < N;j++)
		{
			printf("%4.4f ",a[i][j]);
		}
	printf("\n");
	}
	printf("b is :\n");
	for(i = 0;i < N;i++)
	{
		printf("%4.4f\n",b[i]);
	}
}

void gaosi(double a[][N],double *b)/*高斯变换*/
{
	int i,j,k;
	double l;
	for(k = 0;k < N-1;k++)
	{
		linechange(a,b,k);
		for(i = k+1;i < N;i++)
		{
			l = a[i][k] / a[k][k];
			b[i]-=l*b[k];
			for(j=k;j<N;j++)
			{
				a[i][j]-=l*a[k][j];
			}
		}
	}
}

void linechange(double a[][N],double *b,int k)/*行变换*/ 
{
	int i,j,t;
	double tempa,tempb;
	t=k;
	tempa=a[t][k];
	for(i=k+1;i<N;i++)
	{
		if(fabs(tempa)<fabs(a[i][k]))
		{
			tempa=a[i][k];
			t=i;
		}
	}
	if(t!=k)
	{
		for(j=0;j<N;j++)
		{
			tempa=a[k][j];
			a[k][j]=a[t][j];
			a[t][j]=tempa;
		}
	tempb=b[k];
	b[k]=b[t];
	b[t]=tempb;
	}
}

void huidai(double a[][N],double *b)/*回代求解*/ 
{
	double x[N];
	int i,j;
	double sum;
	x[N-1]=b[N-1]/a[N-1][N-1];
	for(i=N-2;i>=0;i--)
	{
		sum=0;
		for(j=i+1;j<N;j++)
			sum+=a[i][j]*x[j];
		x[i]=(b[i]-sum)/a[i][i];
	}
	printf("resout:\n");
	for(i=0;i<N;i++)
	{
		printf("x[%d]=%4.4f\n",i+1,x[i]);
	}

	double used[3];
	double inc[3];

	used[0] = ( c[2][1] + c[2][2] ) * x[0];
	used[1] = ( c[1][0] + 1 + c[1][1] + c[1][2] ) * x[1];
	used[2] = ( c[0][0] + c[0][1] + 1 + c[0][2] ) * x[2];

	inc[0] = ( 1 - c[2][1] - c[2][2] ) * x[0];
	inc[1] = ( 0 - c[1][0] - c[1][1] - c[1][2] ) * x[1];
	inc[2] = ( 0 - c[0][0] - c[0][1] - c[0][2] ) * x[2];

	for( int count = 0; count <= 2; count++ )
		printf( "used[%d] = %4.4f\n",count+1,used[count] );
	for( int count2 = 0; count2 <= 2; count2++ )
		printf( "incased[%d] = %4.4f\n",count2+1,inc[count2]);
}