数值分析——Jacobi & Gauss-Seidel迭代

怀帝阍而不见

已于 2022-01-27 15:56:43 修改

阅读量1.5k

点赞数

分类专栏：计算数学文章标签： c++

于 2022-01-27 15:55:02 首次发布

本文链接：https://blog.csdn.net/lk32767/article/details/122687732

版权

计算数学专栏收录该内容

4 篇文章 3 订阅

订阅专栏

本系列整理自博主21年秋季学期本科课程数值分析I 的编程作业，内容相对基础，参考书: David Kincaid, Ward Cheney - Numerical Analysis Mathematics of Scientific Computing (2002, Americal Mathematical Society)

算法引入

关于求解线性方程组Ax=b，由高斯消元法得到矩阵的LU分解，称为直接法（direct method）

非直接的方法一般为迭代法（iterative method），适用于矩阵A为稀疏矩阵的情形（sparse system）

引入splitting matrix Q，迭代形式如下：

$Qx^{(k)}=(Q-A)x^{(k-1)}+b$

或写为

$x^{(k)}=Gx^{(k-1)}+b'$

Jacobi Method

Q=diag(A)

分量形式：

$\small x_1^{(k)}=-\frac{1}{a_{11}}(a_{12}x_2^{(k-1)}+...+a_{1n}x^{(k-1)}_n-b_1)$

$\small x_2^{(k)}=-\frac{1}{a_{22}}(a_{21}x_1^{(k-1)}+...+a_{2n}x^{(k-1)}_n-b_2)$

$\small \vdots$

$\small x_n^{(k)}=-\frac{1}{a_{nn}}(a_{n1}x_2^{(k-1)}+...+a_{n,n-1}x^{(k-1)}_{n-1}-b_n)$

代码

1.辅助部分

#include<iostream>
#include<math.h>
#include<iomanip>
#define N 30//dimension
#define K 200//maximum step
using namespace std;
void print(int k, double x[N]) {//output the iteration step and x^(k)
	cout << setw(2) << "[" << setw(2) << k << "]";
	for (int i = 0; i < N; i++)
		cout << setw(10) << x[i];
	cout << endl; cout << endl;
}
int stopcondition(double x1[N], double x2[N]) {
	double temple = 0, M = 0;
	for (int i = 0; i < N; i++) {
		temple = abs(x1[i] - x2[i]);
		if (temple > M)
			M = temple;
	}
	if (M <= 0.00001)
		return 1;
	return 0;
}
void Multiple(double A[N][N], double xk[N]) {
	double b0[N] = { 0 };
	cout << endl;
	cout << "  Ax=";
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			b0[i] += A[i][j] * xk[j];
		}
		cout << setw(10) << b0[i];
	}
	cout << endl; 
}

2. Jacobi

void Jacobi(double A[N][N], double b[N], double x0[N]) {//x0 is the given initial value
	cout << "Jacobi iteration" << endl;
	cout << setw(5) << "[ k]";
	for (int i = 0; i < N; i++)
		cout << setw(4) << "x" << i + 1;
	cout << endl;
	double x1[N] = { 0 }, x2[N] = { 0 }, s = 0;
	//in kth step, x1 stands for x^(k-1),x2 stands for x^(k)
	for (int i = 0; i < N; i++)
		x1[i] = x0[i];
	print(0, x1);
	for (int k = 1; k < K; k++) {
		for (int i = 0; i < N; i++) {
			for (int j = 0; j < N; j++) {
				if (j != i)
					x2[i] -= x1[j] * A[i][j];
			}
			x2[i] += b[i];
			x2[i] = x2[i] / A[i][i];
		}
		print(k, x2);
		if (stopcondition(x1, x2))
			break;
		for (int i = 0; i < N; i++) {
			x1[i] = x2[i];
			x2[i] = 0;
		}

	}
	Multiple(A, x2);
}

Gauss-Seidel Method

Q=A的下三角部分（包括对角）

分量形式

$\small x_1^{(k)}=-\frac{1}{a_{11}}(a_{12}x_2^{(k-1)}+a_{13}x_3^{(k-1)}+...+a_{1n}x^{(k-1)}_n-b_1)$

$\small x_2^{(k)}=-\frac{1}{a_{22}}(a_{21}{\color{Red} x_1^{(k)}}+a_{23}x_3^{(k-1)}+...+a_{2n}x^{(k-1)}_n-b_2)$

$\small x_3^{(k)}=-\frac{1}{a_{33}}(a_{31}{\color{Red} x_1^{(k)}}+a_{32}{\color{Red} x_2^{(k)}}+...+a_{3n}x^{(k-1)}_n-b_3)$

$\small \vdots$

$\small x_n^{(k)}=-\frac{1}{a_{nn}}(a_{n1}{\color{Red} x_1^{(k)}}+a_{n2}{\color{Red} x_2^{(k)}}+...+a_{n,n-1}{\color{Red} x^{(k)}_n}-b_n)$

代码

void GaussSeidel(double A[N][N], double b[N], double x0[N]) {//x0 is the given initial value
	cout << "Gauss-Seidel iteration" << endl;
	cout << setw(5) << "[ k]";
	for (int i = 0; i < N; i++)
		cout << setw(4) << "x" << i + 1;
	cout << endl;
	double x1[N] = { 0 }, x2[N] = { 0 }, s = 0;
	for (int i = 0; i < N; i++)
		x2[i] = x0[i];
	print(0, x2);
	for (int k = 1; k < K; k++) {
		for (int i = 0; i < N; i++) {
			x1[i] = x2[i];//in kth step, x1 stands for x^(k-1)
			s = 0;
			for (int j = 0; j < N; j++) {
				if (j != i)
                //the first i components of x2 equal to that of x^(k)，else equal to x^(k-1)
					s -= x2[j] * A[i][j];
			}
			s += b[i];
			x2[i] = s / A[i][i];
		}
		print(k, x2);
		if (stopcondition(x1, x2))
			break;
	}
	Multiple(A, x2);
}

拓展

（1）Richard Method

Q=I，即

$x^{(k)}=(I-A)x^{(k-1)}+b=x^{(k-1)}+r^{(k-1)}$

（2）SOR Method

Q=α D-C，α∈R，D为任意正定Hermite阵，C为任意满足C+C*=D-A的矩阵。

一般取D=diag(A)，C为A的下三角部分（不包括对角），记1/ω=α，称ω为松弛因子

分量形式

$\small x_1^{(k)}=x_1^{(k-1)}-\frac{\omega}{a_{11}}(a_{11}x_1^{(k-1)}+a_{12}x_2^{(k-1)}+a_{13}x_3^{(k-1)}+...+a_{1n}x^{(k-1)}_n-b_1)$

$\small x_2^{(k)}=x_2^{(k-1)}-\frac{\omega}{a_{22}}(a_{21}{\color{Red} x_1^{(k)}}+a_{22}x_2^{(k-1)}+a_{23}x_3^{(k-1)}+...+a_{2n}x^{(k-1)}_n-b_2)$

$\small x_3^{(k)}=x_3^{(k-1)}-\frac{\omega}{a_{33}}(a_{31}{\color{Red} x_1^{(k)}}+a_{32}{\color{Red} x_2^{(k)}}+a_{33}x_3^{(k-1)}+...+a_{3n}x^{(k-1)}_n-b_3)$

$\small \vdots$

$\small x_n^{(k)}=x_n^{(k-1)}-\frac{\omega}{a_{nn}}(a_{n1}{\color{Red} x_1^{(k)}}+a_{n2}{\color{Red} x_2^{(k)}}+...+a_{n,n-1}{\color{Red} x_{n-1}^{(k)}}+a_{nn}x_n^{(k-1)}+-b_n)$

怀帝阍而不见

关注

0
点赞
踩
11

收藏

觉得还不错? 一键收藏
0
评论
数值分析——Jacobi & Gauss-Seidel迭代

本系列整理自博主21年秋季学期本科课程数值分析I 的编程作业，内容相对基础，参考书: David Kincaid, Ward Cheney - Numerical Analysis Mathematics of Scientific Computing (2002, Americal Mathematical Society)
复制链接

扫一扫