Метод Якобі

Метод Якобі — класичний ітераційний метод розв'язку системи лінійних рівнянь

Опис методу

Для квадратної системи з n лінійних рівнянь:

$A\mathbf {x} =\mathbf {b}$

де:

$A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.$

Матрицю A можна розкласти на два доданки: діагональну матрицю D, та все інше R:

$A=D+R,\qquad D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}},\qquad R={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\a_{21}&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}}.$

Систему лінійних рівнянь можна переписати в вигляді:

$D\mathbf {x} =\mathbf {b} -R\mathbf {x}$

Ітераційний метод Якобі виражається формулою:

$\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)}).$

чи

$x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.$

Збіжність

Метод є збіжним, коли матриця A має домінантну головну діагональ:

$\left|a_{ii}\right|>\sum _{i\neq j}{\left|a_{ij}\right|}.$

другою умовою збіжності є, те щоб спектральний радіус матриці не перевищував одиницю:

$\rho (D^{-1}R)<1.$

Алгоритм

Реалізація на С++

#include <vector>
#include <iostream>
#include <cmath>
using namespace std;

void solve(const vector<float> a, //квадратна матриця
           const vector<float> b, //вектор вільних елементів
           vector<float>& x, //сюди буде записано розв'язок
           const float allowed_error) //допустима похибка
{
    const unsigned n = x.size();
    vector<float> tmp_x(n);
    
    float error;
    
    do
    {
        error = 0;
        
        tmp_x = b;
        for (unsigned i = 0; i < n; ++i)
        {
            for (unsigned j = 0; j < n; ++j)
            {
                if (i != j)
                {
                    tmp_x[i] -= a[i * n + j] * x[j];
                }
            }
            
            //оновити x[i] та порахувати похибку
            const float x_updated = tmp_x[i] / a[i * (n + 1)];
            const float e = fabs(x[i] - x_updated);
            x[i] = x_updated;
            if (e > error) { error = e; }
        }
    }
    while (error > allowed_error);
}

//приклад використання. Користувач вводить вхідні дані.
int main()
{
    unsigned n;
    
    cout << "\nВведіть розмір сиcтеми:\nn = ";
    cin >> n;
    
    vector<float> x(n);
    vector<float> a(n * n);
    vector<float> b(n);
    
    cout << "\nВведіть вектор вільних елементів: \n";
    for (auto& b_elem : b)
    {
        cin >> b_elem;
    }
    
    cout << "\nВведіть коефіцієнти системи: \n";
    for (auto& a_elem : a)
    {
        cin >> a_elem;
    }
    
    float allowed_error;
    cout << "\nВведіть допустиму похибку: \n";
    cin >> allowed_error;
    
    solve(a, b, x, allowed_error);
    
    cout << "\nРозв'язок системи:: \n";
    for (unsigned i = 0; i < n; ++i)
    {
        cout << "\nx[" << i << "]=  " << x[i];
    }
}

Реалізація на С

#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>

void main(){
	int n, i, j, count = 0;
	float **a, *b, *x, *tmp_x, exp, e;

	printf("Введіть розміри сиcтеми:\n");
	printf("n = ");
	scanf("%i", &n);

	a = (float**)malloc(n*sizeof(float*));
	for(i = 0; i < n; i++){
		a[i] = (float*)malloc(n*sizeof(float));		
	}


	b = (float*)malloc(n*sizeof(float));
	x = (float*)malloc(n*sizeof(float));
        tmp_x = (float*)malloc(n*sizeof(float));

	printf("Введіть коефіцієнти системи:\n");
	for(i = 0; i < n; i++){
		for(j = 0; j < n; j++){
                	scanf("%f", &a[i][j]);
		}
	}

	printf("Введіть вектор вільних елементів:\n");
	for(i = 0; i < n; i++){
		scanf("%f", &b[i]);
                x[i] = 0;
	}

	printf("Введіть точність обчислення: ");
	scanf("%f", &e);
	                        
	do{
		count++;

		for(i = 0; i < n; i++){
			tmp_x[i] = 0.0;
			for(j = 0; j < n; j++){
				if(i != j){
					tmp_x[i] = tmp_x[i] + (a[i][j] * x[j]);
				}
			}
			tmp_x[i] = (b[i] - tmp_x[i]) / a[i][i];
		}

		exp = 0;

		for(i = 0; i < n; i++){
			if(fabs(x[i] - tmp_x[i]) > exp){
				exp = fabs(x[i] - tmp_x[i]);
			}
			x[i] = tmp_x[i];
		}
	}while(exp > e);
	
	free(tmp_x);
	for(i = 0; i < n; i++){
		free(a[i]);
	}
	free(a);
	free(b);

        printf("Розв'язок системи:\n");
	for(i = 0; i < n; i++){
		printf("x[%d] = %.6f\n", i+1, x[i]);
	}

	free(x);
}

Див. також

• Ітераційні методи розв'язування СЛАР
• Метод Гауса — Зейделя

Джерела

• Гантмахер Ф. Р. Теорія матриць. — 2024. — 400+ с.(укр.)