Главная
Новости
Строительство
Ремонт
Дизайн и интерьер




27.11.2021


27.11.2021


24.11.2021


24.11.2021


24.11.2021





Яндекс.Метрика





Метод Гаусса — Зейделя решения системы линейных уравнений

29.10.2021

Метод Гаусса — Зейделя (метод Зейделя, процесс Либмана, метод последовательных замещений) — является классическим итерационным методом решения системы линейных уравнений.

Назван в честь Зейделя и Гаусса.

Постановка задачи

Возьмём систему: A x → = b → {displaystyle A{vec {x}}={vec {b}}} , где A = ( a 11 … a 1 n ⋮ ⋱ ⋮ a n 1 … a n n ) , b → = ( b 1 ⋮ b n ) {displaystyle A=left({egin{array}{ccc}a_{11}&ldots &a_{1n}vdots &ddots &vdots a_{n1}&ldots &a_{nn}end{array}} ight),quad {vec {b}}=left({egin{array}{c}b_{1}vdots b_{n}end{array}} ight)}

Или { a 11 x 1 + … + a 1 n x n = b 1 a n 1 x 1 + … + a n n x n = b n {displaystyle left{{egin{array}{rcl}a_{11}x_{1}+ldots +a_{1n}x_{n}&=&b_{1}&&a_{n1}x_{1}+ldots +a_{nn}x_{n}&=&b_{n}end{array}} ight.}

И покажем, как её можно решить с использованием метода Гаусса-Зейделя.

Метод

Чтобы пояснить суть метода, перепишем задачу в виде

{ a 11 x 1 = − a 12 x 2 − a 13 x 3 − … − a 1 n x n + b 1 a 21 x 1 + a 22 x 2 = − a 23 x 3 − … − a 2 n x n + b 2 … a ( n − 1 ) 1 x 1 + a ( n − 1 ) 2 x 2 + … + a ( n − 1 ) ( n − 1 ) x n − 1 = − a ( n − 1 ) n x n + b n − 1 a n 1 x 1 + a n 2 x 2 + … + a n ( n − 1 ) x n − 1 + a n n x n = b n {displaystyle left{{egin{array}{lcr}a_{11}x_{1}&=&-a_{12}x_{2}-a_{13}x_{3}-ldots -a_{1n}x_{n}+b_{1}a_{21}x_{1}+a_{22}x_{2}&=&-a_{23}x_{3}-ldots -a_{2n}x_{n}+b_{2}ldots &&a_{(n-1)1}x_{1}+a_{(n-1)2}x_{2}+ldots +a_{(n-1)(n-1)}x_{n-1}&=&-a_{(n-1)n}x_{n}+b_{n-1}a_{n1}x_{1}+a_{n2}x_{2}+ldots +a_{n(n-1)}x_{n-1}+a_{nn}x_{n}&=&b_{n}end{array}} ight.}

Здесь в j {displaystyle j} -м уравнении мы перенесли в правую часть все члены, содержащие x i {displaystyle x_{i}} , для i > j {displaystyle i>j} . Эта запись может быть представлена как

( L + D ) x → = − U x → + b → , {displaystyle (mathrm {L} +mathrm {D} ){vec {x}}=-mathrm {U} ,{vec {x}}+{vec {b}},}

где в принятых обозначениях D {displaystyle mathrm {D} } означает матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы A {displaystyle A} , а все остальные нули; тогда как матрицы U {displaystyle mathrm {U} } и L {displaystyle mathrm {L} } содержат верхнюю и нижнюю треугольные части A {displaystyle A} , на главной диагонали которых нули.

Итерационный процесс в методе Гаусса — Зейделя строится по формуле

( L + D ) x → ( k + 1 ) = − U x → ( k ) + b → , k = 0 , 1 , 2 , … {displaystyle (mathrm {L} +mathrm {D} ){vec {x}}^{(k+1)}=-mathrm {U} {vec {x}}^{(k)}+{vec {b}},quad k=0,1,2,ldots }

после выбора соответствующего начального приближения x → ( 0 ) {displaystyle {vec {x}}^{(0)}} .

Метод Гаусса — Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения x → ( i ) {displaystyle {vec {x}}^{(i)}} используются здесь сразу же по мере получения, в то время как в методе Якоби они не используются до следующей итерации:

{ x 1 ( k + 1 ) = c 12 x 2 ( k ) + c 13 x 3 ( k ) + … + c 1 n x n ( k ) + d 1 x 2 ( k + 1 ) = c 21 x 1 ( k + 1 ) + c 23 x 3 ( k ) + … + c 2 n x n ( k ) + d 2 … x n ( k + 1 ) = c n 1 x 1 ( k + 1 ) + c n 2 x 2 ( k + 1 ) + … + c n ( n − 1 ) x n − 1 ( k + 1 ) + d n , {displaystyle left{{egin{array}{ccccccccccc}{x_{1}}^{(k+1)}&=&c_{12}{x_{2}^{(k)}}&+&c_{13}x_{3}^{(k)}&+&{ldots }&+&c_{1n}{x_{n}}^{(k)}&+&d_{1}{x_{2}}^{(k+1)}&=&c_{21}{x_{1}^{(k+1)}}&+&c_{23}x_{3}^{(k)}&+&{ldots }&+&c_{2n}{x_{n}}^{(k)}&+&d_{2}ldots &&&&&&&&&&{x_{n}}^{(k+1)}&=&c_{n1}{x_{1}^{(k+1)}}&+&c_{n2}{x_{2}^{(k+1)}}&+&{ldots }&+&c_{n(n-1)}{x_{n-1}}^{(k+1)}&+&d_{n}end{array}} ight.,}

где

c i j = { − a i j a i i , j ≠ i , 0 , j = i . d i = b i a i i , i = 1 , … , n . {displaystyle c_{ij}={egin{cases}-{frac {a_{ij}}{a_{ii}}},&j eq i,,&j=i.end{cases}}quad d_{i}={frac {b_{i}}{a_{ii}}},quad i=1,ldots ,n.}

Таким образом, i-я компонента ( k + 1 ) {displaystyle (k+1)} -го приближения вычисляется по формуле

x i ( k + 1 ) = ∑ j = 1 i − 1 c i j x j ( k + 1 ) + ∑ j = i n c i j x j ( k ) + d i , i = 1 , … , n . {displaystyle x_{i}^{(k+1)}=sum _{j=1}^{i-1}c_{ij}x_{j}^{(k+1)}+sum _{j=i}^{n}c_{ij}x_{j}^{(k)}+d_{i},quad i=1,ldots ,n.}

Например, при n = 3 {displaystyle n=3} :

x 1 ( k + 1 ) = ∑ j = 1 1 − 1 c 1 j x j ( k + 1 ) + ∑ j = 1 3 c 1 j x j ( k ) + d 1 {displaystyle x_{1}^{(k+1)}=sum _{j=1}^{1-1}c_{1j}x_{j}^{(k+1)}+sum _{j=1}^{3}c_{1j}x_{j}^{(k)}+d_{1}} , то есть x 1 ( k + 1 ) = c 11 x 1 ( k ) + c 12 x 2 ( k ) + c 13 x 3 ( k ) + d 1 , {displaystyle x_{1}^{(k+1)}=c_{11}x_{1}^{(k)}+c_{12}x_{2}^{(k)}+c_{13}x_{3}^{(k)}+d_{1},} x 2 ( k + 1 ) = ∑ j = 1 2 − 1 c 2 j x j ( k + 1 ) + ∑ j = 2 3 c 2 j x j ( k ) + d 2 {displaystyle x_{2}^{(k+1)}=sum _{j=1}^{2-1}c_{2j}x_{j}^{(k+1)}+sum _{j=2}^{3}c_{2j}x_{j}^{(k)}+d_{2}} , то есть x 2 ( k + 1 ) = c 21 x 1 ( k + 1 ) + c 22 x 2 ( k ) + c 23 x 3 ( k ) + d 2 , {displaystyle x_{2}^{(k+1)}=c_{21}x_{1}^{(k+1)}+c_{22}x_{2}^{(k)}+c_{23}x_{3}^{(k)}+d_{2},} x 3 ( k + 1 ) = ∑ j = 1 3 − 1 c 3 j x j ( k + 1 ) + ∑ j = 3 3 c 3 j x j ( k ) + d 3 {displaystyle x_{3}^{(k+1)}=sum _{j=1}^{3-1}c_{3j}x_{j}^{(k+1)}+sum _{j=3}^{3}c_{3j}x_{j}^{(k)}+d_{3}} , то есть x 3 ( k + 1 ) = c 31 x 1 ( k + 1 ) + c 32 x 2 ( k + 1 ) + c 33 x 3 ( k ) + d 3 . {displaystyle x_{3}^{(k+1)}=c_{31}x_{1}^{(k+1)}+c_{32}x_{2}^{(k+1)}+c_{33}x_{3}^{(k)}+d_{3}.}

Условие сходимости

Приведём достаточное условие сходимости метода.

Условие окончания

Условие окончания итерационного процесса Зейделя при достижении точности ε {displaystyle varepsilon } в упрощённой форме имеет вид:

∥ x ( k + 1 ) − x ( k ) ∥≤ ε {displaystyle parallel x^{(k+1)}-x^{(k)}parallel leq varepsilon }

Более точное условие окончания итерационного процесса имеет вид

∥ A x ( k ) − b ∥≤ ε {displaystyle parallel Ax^{(k)}-bparallel leq varepsilon }

и требует больше вычислений. Хорошо подходит для разреженных матриц.

Пример реализации на C++

#include <iostream> #include <cmath> using namespace std; // Условие окончания bool converge(double xk[10], double xkp[10], int n, double eps) { double norm = 0; for (int i = 0; i < n; i++) norm += (xk[i] - xkp[i]) * (xk[i] - xkp[i]); return (sqrt(norm) < eps); } double okr(double x, double eps) { int i = 0; double neweps = eps; while (neweps < 1) { i++; neweps *= 10; } int okr = pow(double(10), i); x = int(x * okr + 0.5) / double(okr); return x; } bool diagonal(double a[10][10], int n) { int i, j, k = 1; double sum; for (i = 0; i < n; i++) { sum = 0; for (j = 0; j < n; j++) sum += abs(a[i][j]); sum -= abs(a[i][i]); if (sum > a[i][i]) { k = 0; cout << a[i][i] << " < " << sum << endl; } else { cout << a[i][i] << " > " << sum << endl; } } return (k == 1); } int main() { setlocale(LC_ALL, ""); double eps, a[10][10], b[10], x[10], p[10]; int n, i, j, m = 0; int method; cout << "Введите размер квадратной матрицы: "; cin >> n; cout << "Введите точность вычислений: "; cin >> eps; cout << "Заполните матрицу А: " << endl << endl; for (i = 0; i < n; i++) for (j = 0; j < n; j++) { cout << "A[" << i << "][" << j << "] = "; cin >> a[i][j]; } cout << endl << endl; cout << "Ваша матрица А: " << endl << endl; for (i = 0; i < n; i++) { for (j = 0; j < n; j++) cout << a[i][j] << " "; cout << endl; } cout << endl; cout << "Заполните столбец свободных членов: " << endl << endl; for (i = 0; i < n; i++) { cout << "В[" << i + 1 << "] = "; cin >> b[i]; } cout << endl << endl; /* Ход метода, где: a[n][n] - Матрица коэффициентов x[n], p[n] - Текущее и предыдущее решения b[n] - Столбец правых частей Все перечисленные массивы вещественные и должны быть определены в основной программе, также в массив x[n] следует поместить начальное приближение столбца решений (например, все нули) */ for (int i = 0; i < n; i++) x[i] = 1; cout << "Диагональное преобладание: " << endl; if (diagonal(a, n)) { do { for (int i = 0; i < n; i++) p[i] = x[i]; for (int i = 0; i < n; i++) { double var = 0; for (int j = 0; j < n; j++) if(j!=i) var += (a[i][j] * x[j]); x[i] = (b[i] - var) / a[i][i]; } m++; } while (!converge(x, p, n, eps)); cout << "Решение системы:" << endl << endl; for (i = 0; i < n; i++) cout << "x" << i << " = " << okr(x[i], eps) << "" << endl; cout << "Итераций: " << m << endl; } else { cout << "Не выполняется преобладание диагоналей" << endl; } system("pause"); return 0; }


Пример реализации на Python

import numpy as np def seidel(A, b, eps): n = len(A) x = np.zeros(n) # zero vector converge = False while not converge: x_new = np.copy(x) for i in range(n): s1 = sum(A[i][j] * x_new[j] for j in range(i)) s2 = sum(A[i][j] * x[j] for j in range(i + 1, n)) x_new[i] = (b[i] - s1 - s2) / A[i][i] converge = np.sqrt(sum((x_new[i] - x[i]) ** 2 for i in range(n))) <= eps x = x_new return x