Шпоры по МПиПА / Численные методы / метод Якоби решения СЛАУ / Исходник / SLAU

#include "stdio.h"
#include "iostream.h"
#include "math.h"

// приведение матрицы коэффициентов к виду с ненулевой диагональю и соответствующее изменение вектора правых частей
// возвращает true - если приведение - успешно
bool getDiagonal( double **coefficients, double *rightPart, int currColumn, int numberOfEquation ) {
	bool result = false;
	int i, row = currColumn;
	double tempItem;

	// для матрицы 1x1 ответом является ненулевость ее единственного элемента
	if ( currColumn == numberOfEquation - 1 ) {
		result = coefficients[currColumn][currColumn] != 0;
	}
	else {
		// пока не найдена перестановка приводящая матрицу к виду с ненулевой диагональю и пока мы можем еще просматривать строки ищем перестановку
		while ( !result && row < numberOfEquation ) {
			// если текущий элемент на диагонали нулевой ищем в столбце дальше ненулевой
			if ( coefficients[row][currColumn] != 0 ) {
				// если элемент ненулевой и не лежит на диагонали меняем местами сотвествующие строки в матрице и элементы в векторе прваых частей
				if ( row != currColumn ) {
					for ( i = 0; i < numberOfEquation; i++ ) {
						tempItem = coefficients[currColumn][i];
						coefficients[currColumn][i] = coefficients[row][i];
						coefficients[row][i] = tempItem;
					}				
					tempItem = rightPart[currColumn];
					rightPart[currColumn] = rightPart[row];
					rightPart[row] = tempItem;					
				}
				// рекурсивный вызов фактически для подматрицы с размерностью на 1 меньше
				result = getDiagonal( coefficients, rightPart, currColumn + 1, numberOfEquation );
				if ( result ) {
					break;
				}
			}
			row ++;
		}
	}

	return result;
}

// было ли найдено решение, если да - итог в параметре solution
int iteration( double **coefficients, double *rightPart, int numberOfEquation, double *solution, double precision ) {
	bool result;
	int i, j, k, step = 1;
	double temp;
	double* tempSolution;

	tempSolution = new double[numberOfEquation];
	
	// приведение матрицы коэффициентов к виду с ненулевой диагональю и соответствующее изменение вектора правых частей
	result = getDiagonal( coefficients, rightPart, 0, numberOfEquation );	

	// если приведение успешно - работаем дальше
	if ( result ) {
		double fault = precision + 1;

		// преобразуем матрицу и правую часть для дальнейшего решения
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				if ( i != j ) {
					coefficients[i][j] = - coefficients[i][j] / coefficients[i][i];
				}
			}
			rightPart[i] = rightPart[i] / coefficients[i][i];
			coefficients[i][i] = 0;
		}

		// первое приближение решения - преобразованный вектор правых частей
		for ( i = 0; i < numberOfEquation; i ++ ) {
			solution[i] = rightPart[i];
		}

		// пока не найдено решение с допустимй погрешнстью или пока не исчерпан лимит шагов... если все расходится например
		while ( fault > precision && step <= 1000 ) {

			// поиск новой итерации с ее "самоиспользованием" при расчетах			
			for ( j = 0; j < numberOfEquation; j ++ ) {
				tempSolution[j] = 0;
			}
			for ( i = 0; i < numberOfEquation; i ++ ) {
				for ( j = 0; j < numberOfEquation; j ++ ) {
					tempSolution[i] = tempSolution[i] + coefficients[i][j] * solution[j]; 
				}
				tempSolution[i] = tempSolution[i] + rightPart[i];
			}

			// расчет погрешности
			fault = 0.0;
			for ( j = 0; j < numberOfEquation; j ++ ) {
				fault = fault + (solution[j] - tempSolution[j])*(solution[j] - tempSolution[j]);
			}
			fault = sqrt( fault );

			// сохранение полученной новой итерации
			for ( j = 0; j < numberOfEquation; j ++ ) {
				solution[j] = tempSolution[j];
			}
			step++;
		}
	}
	else {
		result = 1001;
	}
	

	return step;
}


void main() {
	int i, j;
	int size;
	double **coefficients, *rightPart, *solution, precision;

    cout << "Simple iteration method.\nEnter system dimension: ";
    cin >> size;

	coefficients = new double*[size];
	for ( i = 0; i < size; i++ ) {
		coefficients[i] = new double[size];
	}
	rightPart = new double[size];
	solution = new double[size];

	for ( i = 0; i < size; i ++ ){
		cout << "Enter " << i + 1 << " row: ";
		for ( j = 0; j < size; j ++ ){
			cin >> coefficients[i][j];
		}
	}

	cout << "Enter right part: ";
	for ( j = 0; j < size; j ++ ){
		cin >> rightPart[j];
	}

	cout << "Enter precision: ";
		cin >> precision;

	int steps = iteration( coefficients, rightPart, size, solution, precision );
	if ( steps > 1000 ) {
		cout << "Solution for this matrix of coefficients not exist or not found";
	}
	else {
		cout << "Solution is:\n";
		for ( j = 0; j < size; j ++ ){
			cout << solution[j] << "\n";
		}
		cout << "Number of step: " << steps;
	}

    cout << "\nPress \"Enter\" to continue..." << endl; 
    getchar();	
}