数值分析算法之反幂法

这是数值分析中的反幂法代码，有需要的可以下载试用

#include 
#include 
#include 
#define N 3
#define D 0.005

double Abs（double x）
{
	return x >= 0 ? x : -x;
}

double Fan_two（double x[]）
{
	double sum = 0;
	int i;
	for （i = 0; i < N; i++）
	{
		sum += pow（x[i]2）;
	}
	return sqrt（sum）;
}



void * Solve（double （*A）[N] double *bdouble *xdouble （*T）[N]double *tb）
{
	double y[N];
	double s[N];
	double sum = 0;
	double temp;
	double max;
	int i j kt;
	int M[N];

	for （i = 0; i < N; i++）
	{
		for （j = 0; j < N; j++）
			T[i][j] = A[i][j];
	}

	for （i = 0; i < N; i++）
		tb[i] = b[i];

	//计算数组U L
	for （k = 0; k < N; k++）
	{
		//计算中间量//////////////////////////////
		i = k;
		while （i < N）
		{
			for （t = 0; t <= k - 1; t++）
			{
				sum += A[i][t] * A[t][k];
			}
			s[i] = A[i][k] - sum;
			i++;
			sum = 0;
		}
		
		////////////////////////////////////////////

		//选行号///////////////////////////////////////
		max = Abs（s[k]）;
		M[k] = k;
		for （i = k+1; i < N; i++）
		{	
			if （max < Abs（s[i]））
			{
				M[k] = i;
				max = Abs（s[i]）;
			}
		}

		///////////////////////////////////////////////

		//交换//////////////////////////////////////////////
		if （k != M[k]）
		{
			for （t = 0; t <= k - 1; t++）
			{
				temp = A[k][t];
				A[k][t] = A[M[k]][t];
				A[M[k]][t] = temp;
			}

			for （t = k; 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;

		}
		////////////////////

		//计算///////////////////
		j = k + 1;
		i = k + 1;
		A[k][k] = s[k];
		while （1）
		{
			if （k == N - 1）
				break;
			if （N  == j）
				break;
			for （t = 0; t <= k - 1; t++）
			{
				sum += A[k][t] * A[t][j];
			}
			A[k][j] = A[k][j] - sum;
			j++;
			sum = 0;
			
			if （i == N）
				break;
			A[i][k] = s[i] / A[k][k];
			i++;

		}
	}

	//求Qb/////////////////////////////
	for （k = 0; k < N - 1; k++）
	{
		t = M[k];
		temp = b[k];
		b[k] = b[t];
		b[t] = temp;
	}
	//////////////////////////////////////

	//计算解向量
	sum = 0;
	y[0] = b[0];
	for （i = 1; i < N; i++）
	{
		for （t = 0; t <= i - 1; t++）
		{
			sum += A[i][t] * y[t];
		}
		y[i] = b[i] - sum;
		sum = 0;
	}

	x[N-1] = y[N-1] / A[N-1][N-1];
	for （i = N - 2; i >= 0; i--）
	{
		for （t = i+1; t <= N-1; t++）
		{
			sum += A[i][t] * x[t];
		}
		x[i] = （y[i] - sum） / A[i][i];
		sum = 0;
	}

	for （i = 0; i < N; i++）
	{
		for （j = 0; j < N; j++）
			A[i][j] = T[i][j];
	}
	for （i = 0; i < N; i++）
		b[i] = tb[i];
	return NULL;
}

void main（）
{
	double u[N] = {1.0 1.0 1.0};
	double A[N][N] = {6.0 -12.0 6.0 -21.0 -3.0 24.0 -12.0 -12.0 51};
	double T[N][N];
	double ty[N];
	double y[N];
	double ita;
	double bta = 1;
	double sum = 0;
	double temp = 0;
	int ijk = 0;
	while （Abs（1/bta - 1/temp） / Abs（1/bta） > D）
	{
		printf（“k = %dbta = %lf\n“kbta）;
		for （i = 0; i < N; i++）
		{
			pri

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      33792  2011-10-23 13:48  反幂法\Debug\vc60.idb

     文件      45056  2011-10-23 13:48  反幂法\Debug\vc60.pdb

     文件     217145  2011-10-23 13:48  反幂法\Debug\反幂法.exe

     文件     226700  2011-10-23 13:48  反幂法\Debug\反幂法.ilk

     文件       9472  2011-10-23 13:48  反幂法\Debug\反幂法.obj

     文件     223872  2011-10-21 21:32  反幂法\Debug\反幂法.pch

     文件     459776  2011-10-23 13:48  反幂法\Debug\反幂法.pdb

     文件       3643  2011-10-23 13:48  反幂法\反幂法.cpp

     文件       3403  2011-10-21 16:30  反幂法\反幂法.dsp

     文件        520  2011-10-21 16:51  反幂法\反幂法.dsw

     文件      41984  2011-10-23 13:48  反幂法\反幂法.ncb

     文件      48640  2011-10-23 13:48  反幂法\反幂法.opt

     文件       1233  2011-10-23 13:48  反幂法\反幂法.plg

     目录          0  2011-10-23 13:48  反幂法\Debug

     目录          0  2011-10-23 13:48  反幂法

----------- ---------  ---------- -----  ----

              1315236                    15