Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://num-anal.srcc.msu.ru/lib_na/cat/as_htm_c/ass2r_c.htm
Дата изменения: Fri Dec 4 10:59:17 2015 Дата индексирования: Sun Apr 10 03:08:48 2016 Кодировка: Windows-1251 |
Текст подпрограммы и версий ass2r_c.zip , ass2d_c.zip |
Тексты тестовых примеров tass2r_c.zip , tass2d_c.zip |
Решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя (матрица системы представлена в формате RR (LU) U) .
Сокращенное название формата RR (LU) U происходит от английского словосочетания "Row - wise Representation, Lower - Upper, Unordered" (строчное представление, нижний треугольник - верхний треугольник, неупорядоченное).
Данный формат используется для прямоугольных (квадратных) матриц, у которых все или большинство диагональных элементов не равны нулю. Матрица A в этом формате представляется в виде суммы L + D + U, где L - нижняя треугольная, D - диагональная, U - верхняя треугольная матрицы. Диагональные элементы хранятся в отдельном одномерном массиве AD, а элементы матриц L и U содержатся в одномерном массиве AN и связаны между собой списками (портретами) IA и JA.
Например, пусть задана матрица
| 4 0 1 0 2 | | 1 2 0 0 0 | A = | 0 0 2 1 0 | | 1 1 0 1 1 | | 0 0 0 0 16 |
Тогда в формате RR (LU) U данная матрица с точностью до упорядоченности ее ненулевых внедиагональных элементов имеет вид:
IA = ( 1; 3; 4; 5; 8; 8 ) JA = ( 5, 3; 1; 4; 5, 1, 2 ) AN = ( 2, 1, 1, 1, 1, 1, 1 ) AD = ( 4, 2, 2, 1, 16 )
Пусть квадратная матрица A порядка n линейной системы Ax = b задана в формате RR (LU) U и пусть известно, что все диагональные элементы матрицы A не равны нулю.
Метод Гаусса - Зейделя заключается в следующем. Перепишем i - е уравнение системы в виде:
i -1 n ∑ ai j x j + ai i x i + ∑ ai j x j = b i j =1 j= i +1
Поскольку по предположению все ai i ≠ 0, то это уравнение можно записать в виде:
i -1 n x = ( b i - ∑ ai j x j - ∑ ai j x j ) / ai i j =1 j= i +1
Это уравнение решается итерационно, если задать начальное приближение к решению xj (1) ( j = 1, 2, ..., n). В случае, когда начальное приближение xj (1) заранее неизвестно, то можно положить
xj(1) = b j / aj j , j = 1, 2, ..., n .
Следует отметить, что скорость сходимости итерационного процесса Гаусса - Зейделя может быть значительно увеличена, если начальное приближение к решению выбрано удачно.
Далее на m - м итерационном шаге процесса Гаусса - Зейделя имеем (m = 2, 3, 4, ...)
i -1 xi(m) = ( bi - ∑ ai j xj(m) - j =1 n - ∑ ai j xj(m -1) ) / ai i j= i +1
Заметим, что здесь вновь вычисленные компоненты решения на m - м шаге сразу же используются для вычисления следующих компонент.
В некоторых случаях сходимость метода Гаусса - Зейделя можно ускорить, если использовать прием верхней релаксации:
xi(m) = xi(m -1) + Q( xi(m) - xi(m -1) )
Релаксационный множитель Q обычно выбирается в пределах от 1 до 2. Если матрица положительно определенная, то Q выбирается в пределах от 0 < Q < 1. Иногда в качестве Q используется диагональная матрица, позволяющая подбирать релаксационные множители индивидуально для каждого уравнения системы.
Метод Гаусса - Зейделя хорошо сходится, если матрица A является нижней или "почти нижней" треугольной. Для его сходимости не требуется диагонального преобладания в матрице A. Необходимым и достаточным условием этого метода является следующее: все корни уравнения
| a1 1λ a1 2 a1 3 ... a1 n | | a2 1λ a2 2λ a2 3 ... a2 n | det | . . . . . . . . . . . . . . . = 0 | an 1λ an 2λ an 3λ ... an nλ |
должны быть по модулю меньше 1. Если матрица A является вещественной симметричной положительно определенной матрицей, то метод Гаусса - Зейделя всегда сходится. Скорость сходимости будет быстрее, если матрица A близка к диагональной.
Описанный метод Гаусса - Зейделя реализован в виде подпрограммы ass2r_c. В этой подпрограмме не требуется задавать на входе начальное приближение к решению: оно вычисляется в ass2r_c по формуле xi (1) = bi /ai i. Подпрограмма предполагает задание максимально допустимого количества итераций. Контроль точности ведется по абсолютной погрешности EPS: текущая итерация считается приемлемым приближением к решению, если | xi (m) - xi (m - 1) | < EPS для всех i.
Если за заданное максимальное количество итераций требуемая точность не достигнута, то пользователю подпрограмма выдает соответствующее сообщение. В случае, когда решение системы вычислено с заданной точностью, можно получить информацию о реально выполненном количестве итераций.
Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.
int ass2r_c (integer *ia, integer *ja, real *an, real *ad, integer *n, real *b, real *x, real *q, real *eps, integer *itmax, integer *iflag)
Параметры
ia, ja, - an | заданные портрет и ненулевые внедиагональные элементы матрицы A в формате RR (LU) U; |
ad - | вещественный одномерный массив длины n, содержащий диагональные элементы матрицы A; требуется, чтобы ни один диагональный элемент не был равен нулю; |
n - | заданный порядок матрицы A (тип: целый); |
b - | вещественный одномерный массив длины n, содержащий вектор правой части системы; |
x - | вещественный одномерный массив длины n, на выходе из подпрограммы содержащий вычисленные компоненты вектора решения системы; |
q - | заданный релаксационный множитель (тип: вещественный); |
eps - | заданная абсолютная погрешность, с которой требуется вычислить решение; |
itmax - | заданное максимальное количество итераций по методу Гаусса - Зейделя; |
iflag - | целая переменная, служащая для сообщения о режиме окончания работы подпрограммы: |
iflag=0 - | решение системы вычислено с заданной точностью; |
iflag=1 - | решение не получено с заданной точностью за заданное максимальное количество итераций. |
Версии
ass2d_c - | решение невырожденной разреженной линейной системы итерационным методом Гаусса - Зейделя в режиме удвоенной точности; при этом параметры an, ad, b, x, q и eps должны иметь тип double. |
Вызываемые подпрограммы нет
Замечания по использованию
Подпрограммы ass2r_c и ass2d_c имеют внешнюю структуру с именем ass2rr_ , содержащую элемент целого типа iter. Значение переменной iter полагается равным количеству в действительности выполненных итераций, если iflag = 0. Если iflag = 1, то iter = itmax. |
struct { int iter; } ass2rr_; #define ass2rr_1 ass2rr_ int main(void) { /* Initialized data */ static int ia[6] = { 1,2,3,5,6,8 }; static int ja[7] = { 5,1,2,1,2,3,1 }; static float an[7] = { 1.f,1.f,1.f,1.f,1.f,1.f,2.f }; static float ad[5] = { 4.f,2.f,2.f,8.f,16.f }; static float b[5] = { 1.f,1.f,1.f,1.f,1.f }; /* System generated locals */ int i__1; /* Local variables */ extern int ass2r_c(int *, int *, float *, float *, int *, float *, float *, float *, float *, int *, int *); static int i__, n; static float q; static int iflag; static float x[5]; static int itmax; static float eps; n = 5; q = 1.5f; i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { /* l1: */ x[i__ - 1] = b[i__ - 1] / ad[i__ - 1]; } eps = .001f; itmax = 500; ass2r_c(ia, ja, an, ad, &n, b, x, &q, &eps, &itmax, &iflag); printf("\n %5i \n", iflag); printf("\n %5i \n", ass2rr_1.iter); printf("\n %15.6e %15.6e %15.6e %15.6e %15.6e \n", x[0], x[1], x[2], x[3], x[4]); return 0; } /* main */ Результаты: iflag = 0 , iter = 7 , x = (0.245396, 0.377041, 0.188364, 0.0778308, 0.0203379)