jgssor.cpp

该程序用于数值分析中计算J迭代

#include <stdio.h>

double Ab[50][50], x0[50], exact[50];
double Jres[50][50], GSres[50][50], SORres[50][50];
double Jpre[50], GSpre[50], SORpre[50];
int num, N; // num denote the number of coefficient
			// N denote the iterate times
double w;

void input( )
{
	scanf( "%d", &num );
	int i, j;

	for( i = 0; i < num; i++ )
		for( j = 0; j <= num; j++ )
			scanf( "%lf", &Ab[i][j] );
	for( i = 0; i < num; i++ )
		scanf( "%lf", &x0[i] );
	for( i = 0; i < num; i++ )
		scanf( "%lf", &exact[i] );

	scanf( "%d", &N );
	scanf( "%lf", &w ); 
}
 
double getAbs( double d )
{
	return d < 0 ? -d : d;
}

void JIterate( )
{
	int iter, i;

	for( iter = 1; iter <= N; iter++ ) {
		double max;

		for( i = 0; i < num; i++ ) {
			int bi, ai; // bi denote the coefficient before iterator i
						// ai the same
			double total = 0.0;
			for( bi = 0; bi < i; bi++ )
				total += Ab[i][bi] * Jres[iter-1][bi];
			for( ai = i + 1; ai < num; ai++ )
				total += Ab[i][ai] * Jres[iter-1][ai];

			Jres[iter][i] = Ab[i][num] - total;

			// the following is used to find the precision this iterate
			if( i == 0 ) max = getAbs( Jres[iter][0] - exact[0] );
			else {
				double tmp = getAbs( Jres[iter][i] - exact[i] );
				max = max < tmp ? tmp : max;
			}
		}

		Jpre[iter] = max;
	}
}

void GSIterate( )
{
	int iter, i;

	for( iter = 1; iter <= N; iter++ ) {
		double max;

		for( i = 0; i < num; i++ ) {
			int bi, ai; // bi denote the coefficient before iterator i
						// ai the same
			double total = 0.0;
			for( bi = 0; bi < i; bi++ )
				total += Ab[i][bi] * GSres[iter][bi];
			for( ai = i + 1; ai < num; ai++ )
				total += Ab[i][ai] * GSres[iter-1][ai];

			GSres[iter][i] = Ab[i][num] - total;

			// the following is used to find the precision this iterate
			if( i == 0 ) max = getAbs( GSres[iter][0] - exact[0] );
			else {
				double tmp = getAbs( GSres[iter][i] - exact[i] );
				max = max < tmp ? tmp : max;
			}
		}
 
		GSpre[iter] = max;
	}
}

void SORIterate( )
{
	int i, j, iter;

	for( i = 0; i < num; i++ )
		for( j = 0; j <= num; j++ )
			Ab[i][j] *= w;

	for( iter = 1; iter <= N; iter++ ) {
		double max;

		for( i = 0; i < num; i++ ) {
			int bi, ai; // bi denote the coefficient before iterator i
						// ai the same
			double total = 0.0;
			for( bi = 0; bi < i; bi++ )
				total += Ab[i][bi] * SORres[iter][bi];
			for( ai = i; ai < num; ai++ )
				total += Ab[i][ai] * SORres[iter-1][ai];

			SORres[iter][i] = Ab[i][num] - total + SORres[iter-1][i];

			// the following is used to find the precision this iterate
			if( i == 0 ) max = getAbs( SORres[iter][0] - exact[0] );
			else {
				double tmp = getAbs( SORres[iter][i] - exact[i] );
				max = max < tmp ? tmp : max;
			}
		}
 
		SORpre[iter] = max;
	}
}

void init( )
{
	int i, j;

	for( i = 0; i < num; i++ ) {
		Jres[0][i] = GSres[0][i] = SORres[0][i] = 0;
	}
	for( i = 0; i < num; i++ ) {
		double tmp = Ab[i][i];

		for( j = 0; j <= num; j++ )
			Ab[i][j] /= tmp;
	}
}

void output( double rarr[50][50], double parr[50] )
{
	for( int i = 1; i <= N; i++ ) {
		for( int j = 0; j < num; j++ )
			printf( "%lf ", rarr[i][j] );
		printf( "%lf\n", parr[i] );
	}
}

int main( )
{
	freopen( "test.in", "r", stdin );
	freopen( "test.out", "w", stdout );

	input( );
	init( );

	JIterate( );
	output( Jres, Jpre );

	putchar( '\n' );
	GSIterate( );
	output( GSres, GSpre );

	putchar( '\n' );
	SORIterate( );
	output( SORres, SORpre ); 

	return 0;
}