Программирование. Принципы и практика использования C++ Исправленное издание
вернуться

Страуструп Бьерн
Шрифт:

void classical_elimination(Matrix& A,Vector& b)

{

const Index n = A.dim1;

// проходим от первого столбца до последнего,

// обнуляя элементы, стоящие ниже диагонали:

for (Index j = 0; j<n–1; ++j) {

const double pivot = A(j, j);

if (pivot == 0) throw Elim_failure(j);

// обнуляем элементы, стоящие ниже диагонали в строке i

for (Index i = j+1; i<n; ++i) {

const double mult = A(i, j) / pivot;

A[i].slice(j) = scale_and_add(A[j].slice(j),

–mult, A[i].slice(j));

b(i) –= mult * b(j); // изменяем вектор b

}

}

}

Опорным называется элемент, лежащий на диагонали в строке, которую мы в данный момент обрабатываем. Он должен быть ненулевым, потому что нам придется на него делить; если он равен нулю, то генерируется исключение.

Vector back_substitution(const Matrix& A, const Vector& b)

{

const Index n = A.dim1;

Vector x(n);

for (Index i = n – 1; i >= 0; ––i) {

double s = b(i)–dot_product(A[i].slice(i+1),x.slice(i+1));

if (double m = A(i, i))

x(i) = s / m;

else

throw Back_subst_failure(i);

}

return x;

}

24.6.2. Выбор ведущего элемента

Для того чтобы избежать проблем с нулевыми диагональными элементами и повысить устойчивость алгоритма, можно переставить строки так, чтобы нули и малые величины на диагонали не стояли. Говоря “повысить устойчивость”, мы имеем в виду понижение чувствительности к ошибкам округления. Однако по мере выполнения алгоритма элементы матрицы будут изменяться, поэтому перестановку строк приходится делать постоянно (иначе говоря, мы не можем лишь один раз переупорядочить матрицу, а затем применить классический алгоритм).

void elim_with_partial_pivot(Matrix& A, Vector& b)

{

const Index n = A.dim1;

for (Index j = 0; j < n; ++j) {

Index pivot_row = j;

// ищем подходящий опорный элемент:

for (Index k = j + 1; k < n; ++k)

if (abs(A(k, j)) > abs(A(pivot_row, j))) pivot_row = k;

// переставляем строки, если найдется лучший опорный

// элемент

if (pivot_row != j) {

A.swap_rows(j, pivot_row);

std::swap(b(j), b(pivot_row));

}

// исключение:

for (Index i = j + 1; i < n; ++i) {

const double pivot = A(j, j);

if (pivot==0) error("Решения нет: pivot==0");

onst double mult = A(i, j)/pivot;

A[i].slice(j) = scale_and_add(A[j].slice(j),

–mult, A[i].slice(j));

b(i) –= mult * b(j);

}

}

}

Для того чтобы не писать циклы явно и привести код в более традиционный вид, мы используем функции

swap_rows
и
scale_and_multiply
.

24.6.3. Тестирование

Очевидно, что мы должны протестировать нашу программу. К счастью, это сделать несложно.

void solve_random_system(Index n)

{

Matrix A = random_matrix(n); // см. раздел 24.7

Vector b = random_vector(n);

cout << "A = " << A << endl;

cout << "b = " << b << endl;

try {