Метод Якобі — класичний ітераційний метод розв'язку системи лінійних рівнянь
Опис методу
Для квадратної системи з n лінійних рівнянь:
де:
Матрицю A можна розкласти на два доданки: діагональну матрицю D, та все інше R:
Систему лінійних рівнянь можна переписати в вигляді:
Ітераційний метод Якобі виражається формулою:
чи
Збіжність
Достатньою, але не необхідною, є умова, коли матриця A має домінантну головну діагональ:
Необхідною і достатньою умовою збіжності є, те щоб спектральний радіус матриці не перевищував одиницю:
Алгоритм
Реалізація на С++
#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Введіть розмір системи:\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("Введіть розміри системи:\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);
}
Див. також
Джерела