线性方程组的迭代法(Jacobi 迭代法和Gauss-Seidel 迭代法) C++代码

Jacobi 迭代法

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

// 定义方程组的系数矩阵和常数向量
vector<vector<double>> A = {{20, 2, 3},
                            {1, 8, 1},
                            {2, -3, 15}};
vector<double> b = {24, 12, 30};

// 定义迭代次数和精度阈值
int maxIterations = 100;
double epsilon = 5e-5;
int iterations = 0;
// 雅可比迭代函数
vector<double> jacobiIteration(const vector<vector<double>>& A, const vector<double>& b, const vector<double>& x0) {
    int n = A.size();
    vector<double> x(x0);  // 初始解的近似值
    
    for (int k = 0; k < maxIterations; k++) {
        iterations++;
        vector<double> x_new(n, 0);
        
        for (int i = 0; i < n; i++) {
            double sum = 0;
            for (int j = 0; j < n; j++) {
                if (j != i) {
                    sum += A[i][j] * x[j];
                }
            }
            x_new[i] = (b[i] - sum) / A[i][i];
        }
        
        // 判断是否满足终止条件
        double diff = 0;
        for (int i = 0; i < n; i++) {
            diff += abs(x_new[i] - x[i]);
        }
        if (diff < epsilon) {
            break;
        }
        
        // 更新解的近似值
        x = x_new;
    }
    return x;
}

int main() {
    int n = A.size();
    vector<double> x0(n, 0);  // 初始解的近似值
    
    vector<double> x = jacobiIteration(A, b, x0);
    
    cout << "Solution: \n";
    for (int i = 0; i < n; i++) {
        cout << "x" << i+1 << " = " << x[i] << "\n";
    }
    cout << endl;
    cout << "iterations = " << iterations << '\n';
    return 0;
}

Gauss-Seidel 迭代法

#include <iostream>
#include <cmath>

#define N 3 // 线性方程组的未知数个数
#define MAX_ITERATIONS 100 // 最大迭代次数
#define EPSILON 0.00005 // 迭代停止的精度

void gaussSeidel(double coef[N][N], double b[N], double x[N]) {
    double x_new[N];

    // 初始化迭代结果
    for (int i = 0; i < N; i++) {
        x[i] = 0;
    }

    int iterations = 0;
    double error = EPSILON + 1;

    while (error > EPSILON && iterations < MAX_ITERATIONS) {
        for (int i = 0; i < N; i++) {
            double sum1 = 0;
            for (int j = 0; j < i; j++) {
                sum1 += coef[i][j] * x_new[j];
            }

            double sum2 = 0;
            for (int j = i + 1; j < N; j++) {
                sum2 += coef[i][j] * x[j];
            }

            x_new[i] = (b[i] - sum1 - sum2) / coef[i][i];
        }

        error = 0;
        for (int i = 0; i < N; i++) {
            error += std::abs(x_new[i] - x[i]);
            x[i] = x_new[i];
        }

        iterations++;
    }

    if (iterations < MAX_ITERATIONS) {
        std::cout << "Converged in " << iterations << " iterations." << std::endl;
    } else {
        std::cout << "Did not converge within the maximum number of iterations." << std::endl;
    }
}

int main() {
    double coef[N][N] = {{20, 2, 3}, {1, 8, 1}, {2, -3, 15}};
    double b[N] = {24, 12, 30};
    double x[N];

    gaussSeidel(coef, b, x);

    std::cout << "Solution:" << std::endl;
    for (int i = 0; i < N; i++) {
        std::cout << "x" << i << " = " << x[i] << std::endl;
    }

    return 0;
}

Jacobi 迭代法与Gauss-Seidel 迭代法的比较

#include <iostream>
#include <vector>
#include <cmath>

#define N 3 // 线性方程组的未知数个数
#define MAX_ITERATIONS 100 // 最大迭代次数
#define EPSILON 0.00005 // 迭代停止的精度

// 高斯-赛德尔迭代法
int gaussSeidel(const std::vector<std::vector<double>>& A, const std::vector<double>& b, std::vector<double>& x) {
    std::vector<double> x_new(N);

    // 初始化迭代结果
    for (int i = 0; i < N; i++) {
        x[i] = 0;
    }

    int iterations = 0;
    double error = EPSILON + 1;

    while (error > EPSILON && iterations < MAX_ITERATIONS) {
        for (int i = 0; i < N; i++) {
            double sum1 = 0;
            for (int j = 0; j < i; j++) {
                sum1 += A[i][j] * x_new[j];
            }

            double sum2 = 0;
            for (int j = i + 1; j < N; j++) {
                sum2 += A[i][j] * x[j];
            }

            x_new[i] = (b[i] - sum1 - sum2) / A[i][i];
        }

        error = 0;
        for (int i = 0; i < N; i++) {
            error += std::abs(x_new[i] - x[i]);
            x[i] = x_new[i];
        }

        iterations++;
    }

    return iterations;
}

// 雅可比迭代法
int jacobi(const std::vector<std::vector<double>>& A, const std::vector<double>& b, std::vector<double>& x) {
    std::vector<double> x_new(N);

    // 初始化迭代结果
    for (int i = 0; i < N; i++) {
        x[i] = 0;
    }

    int iterations = 0;
    double error = EPSILON + 1;

    while (error > EPSILON && iterations < MAX_ITERATIONS) {
        for (int i = 0; i < N; i++) {
            double sum = 0;
            for (int j = 0; j < N; j++) {
                if (j != i) {
                    sum += A[i][j] * x[j];
                }
            }

            x_new[i] = (b[i] - sum) / A[i][i];
        }

        error = 0;
        for (int i = 0; i < N; i++) {
            error += std::abs(x_new[i] - x[i]);
            x[i] = x_new[i];
        }

        iterations++;
    }

    return iterations;
}

int main() {
    std::vector<std::vector<double>> A = {{20, 2, 3}, {1, 8, 1}, {2, -3, 15}};
    std::vector<double> b = {24, 12, 30};
    std::vector<double> x(N);

    int gaussSeidelIterations = gaussSeidel(A, b, x);
    int jacobiIterations = jacobi(A, b, x);

    std::cout << "Gauss-Seidel iterations: " << gaussSeidelIterations << std::endl;
    std::cout << "Jacobi iterations: " << jacobiIterations << std::endl;

    return 0;
}