| Flame of Soul |
04.02.2009 13:54 |
Метод простых итераций довольно медленно сходится. Для его ускорения существует метод Зейделя
Цитата:
Входные данные - матрица коэффициентов СЛАУ, вектор правых частей, размерность системы, заданная точность.
Выходные данные - вектор решения.
Вспомогательная функция getDiagonal( double **coefficients, double *rightPart, int currColumn, int numberOfEquation ) преобразует матрицу СЛАУ к виду с ненулевыми диагональными элементами с параллельным преобразованием вектора правых частей.
Не всякая система может быть приведена к системе с матрицей с ненулевой диагональю. Если такое приведение не найдено ( в данном случае оно ищется перестановкой строк ) – система не решается.
Быстрее метода простых итераций.
Даже если система приведена к виду с ненулевой диагональю решение может быть не найдено. Это зависит от нормы матрицы и первого приближения решения.
|
Исходный код:
PHP код:
#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 < i; j ++ ) {
tempSolution[i] = tempSolution[i] + coefficients[i][j] * tempSolution[j];
}
for ( j = i; j < numberOfEquation; j ++ ) {
tempSolution[i] = tempSolution[i] + coefficients[i][j] * solution[j];
}
tempSolution[i] = tempSolution[i] + rightPart[i];
}
// расчет погрешности
fault = 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();
}
Примеры использования (распечатка листинга):
Цитата:
Simple iteration method.
Enter system dimension: 3
Enter 1 row: 10 1 1
Enter 2 row: 2 10 1
Enter 3 row: 2 2 10
Enter right part: 12 13 14
Enter precision: 0.1
Solution is:
1.00068
0.997944
1.00028
Number of step: 3
Press "Enter" to continue...
|
PS: Больше не чем помочь немогу
|
