Основы вычислительной математики. Денисова Э.В - 48 стр.

UptoLike

ВУЗ: 

Университет ИТМО | Санкт-Петербург

Составители: 

Рубрика: 

Программистский вариант метода Гаусса имеет три отличия от математического:
1.
индексы строк и столбцов матрицы начинаются с нуля, а не с единицы;
2.
недостаточно найти просто ненулевой элемент в столбце. В программировании все действия с
вещественными числами производятся приближенно, поэтому можно считать, что точного равенства
вещественных чисел вообще не бывает. Некоторые компиляторы даже выдают предупреждения на каждую
операцию проверки равенства вещественных чисел. Поэтому вместо проверки на равенство нулю числа a
ij
следует сравнивать его абсолютную величину a
ij
с очень маленьким числом ε (например, ε = 0.00000001).
Если a
ij
=< ε, то следует считать элемент a
ij
нулевым;
3.
при обнулении элементов j-го столбца, начиная со строки i + 1, мы к k-й строке, где k > i, прибавляем i-ю
строку, умноженную на коэффициент r = -a
kj
/a
ij
:
r = -a
kj
/a
ij
.
a
k
= a
k
+ r * a
i
Такая схема работает нормально только тогда, когда коэффициент r по абсолютной величине не превосходит
единицы. В противном случае, ошибки округления умножаются на большой коэффициент и, таким образом,
экспоненциально растут. Математики называют это явление неустойчивостью вычислительной схемы. Если
вычислительная схема неустойчива, то полученные с ее помощью результаты не имеют никакого отношения к
исходной
задаче. В нашем случае схема устойчива, когда коэффициент r = -a
kj
/a
ij
не превосходит по модулю
единицы. Для этого должно выполняться неравенство
a
ij
>= a
kj
при k > i
Отсюда следует, что при поиске разрешающего элемента в j-м столбце необходимо найти не первый
попавшийся ненулевой элемент, а максимальный по абсолютной величине. Если он по модулю не
превосходит ε, то считаем, что все элементы столбца нулевые; иначе меняем местами строки, ставя его на
вершину столбца, и затем обнуляем столбец элементарными преобразованиями
второго рода.
Ниже дан полный текст программы на Си, приводящей вещественную матрицу к ступенчатому виду.
Функция, реализующая метод Гаусса, одновременно подсчитывает и ранг матрицы. Программа вводит
размеры матрицы и ее элементы с клавиатуры и вызывает функцию приведения к ступенчатому виду. Затем
программа печатает ступенчатый вид матрицы и ее ранг. В
случае квадратной матрицы также вычисляется и
печатается определитель матрицы, равный произведению диагональных элементов ступенчатой матрицы.
При реализации метода Гаусса используется схема построения цикла с помощью инварианта,. В цикле
меняются две переменные -- индекс строки i, 0 =< i < m - 1, и индекс столбца j, 0 =< j < n - 1. Инвариантом
цикла является утверждение о том, что часть матрицы (математики говорят минор) в столбцах 0,1,...j - 1
приведена
к ступенчатому виду и что первый ненулевой элемент в строке i - 1 стоит в столбце с индексом
меньшим j. В теле цикла рассматривается только минор матрицы в строках i,...,m - 1 и столбцах j,...,n - 1.
Сначала ищется максимальный по модулю элемент в j-м столбце. Если он по абсолютной величине не
превосходит ε, то j увеличивается на единицу (считается, что
столбец нулевой). Иначе перестановкой строк
разрешающий элемент ставится на вершину j-го столбца минора, и затем столбец обнуляется элементарными
преобразованиями второго рода. После этого оба индекса i и j увеличиваются на единицу. Алгоритм
завершается, когда либо i = m, либо j = n. По окончании алгоритма значение переменной i равно числу
ненулевых строк ступенчатой матрицы, т.е. рангу исходной матрицы.
Для
вычисления абсолютной величины вещественного числа x типа double мы пользуемся стандарной
математической функцией fabs(x), описанной в стандартном заголовочном файле "math.h.
#include <stdio.h> // Описания функций ввода-вывода
#include <math.h> // Описания математических функций
#include <stdlib.h> // Описания функций malloc и free
// Прототип функции приведения матрицы
// к ступенчатому виду.
// Функция возвращает ранг матрицы
int gaussMethod(
int m, // Число строк матрицы
int n, // Число столбцов матрицы
double *a, // Адрес массива элементов матрицы
double eps // Точность вычислений
);
int main() {
int m, n, i, j, rank;
double *a;
double eps, det;
printf("Введите размеры матрицы m, n: ");

Страницы