gauss.cpp

数值计算的例子和参考书籍

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

const double eps = 1e-10;

const int MAXN = 1000;

bool gauss(double a[][MAXN], int n, double b[])
{
	int i, j, k, ik;
	double det = 1.0;

	for(k = 1; k <= n - 1; k++)
	{
		double max = 0;
		ik = 0;
		for(i = k; i <= n; i++)
		{
			if(fabs(a[i][k]) > max)
			{
				max = fabs(a[i][k]);
				ik = i;
			}
		}
		if(max < eps)
		{
			det = 0;
			return false;
		}
		if(ik != k)
		{
			det = -det;
			for(i = 1; i <= n; i++) swap(a[k][i], a[ik][i]);
			swap(b[k], b[ik]);
		}
		for(i = k + 1; i <= n; i++)
		{
			a[i][k] /= a[k][k];
			for(j = k + 1; j <= n; j++) a[i][j] -= a[i][k] * a[k][j];
			b[i] -= a[i][k] * b[k];
		}
		det = a[k][k] * det;
	}
	if(fabs(a[n][n]) < eps)
	{
		det = 0;
		return false;
	}
	//回带求解过程
	b[n] = b[n] / a[n][n];
	for(i = n - 1; i >= 1; i--)
	{
		double sum = 0.0;
		for(j = i + 1; j <= n; j++) sum += a[i][j] * b[j];
		b[i] = (b[i] - sum) / a[i][i];
	}
	det = a[n][n] * det;
	
	for(i = 1; i <= n; i++) printf("%lf\n", b[i]);

	return true;
}

double a[MAXN][MAXN], b[MAXN];
int n;
int main(void)
{
	int i, j;
	while(scanf("%d", &n) != EOF)
	{
		for(i = 1; i <= n; i++) for(j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
		for(i = 1; i <= n; i++) scanf("%lf", &b[i]);
		gauss(a, n, b);
	}
	return 0;
}