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;
}