离散Poisson系统的常用处理方法以及比较(Python实现)

最新推荐文章于 2024-04-11 13:55:14 发布

gdx6662020

最新推荐文章于 2024-04-11 13:55:14 发布

阅读量815

点赞数

分类专栏：偏微分方程数值解文章标签：算法

本文链接：https://blog.csdn.net/liuxinglongqi/article/details/114591896

版权

偏微分方程数值解专栏收录该内容

14 篇文章 7 订阅

订阅专栏

介绍

在二维Poisson方程五点差分格式及简单求解方法Python实现一文中给出了二维齐次 Dirichlet 型边界条件的 Poisson 方程的标准五点差分格式离散, 得到了离散系统:
$A X = b$
其中 $A$ 为离散后的微分矩阵:
对于一维情形:
$\frac{1}{h^{2}}A_{1} = \frac{1}{h^{2}} \begin{pmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{pmatrix}$
对于二维情形:
$A_{2} = \frac{1}{h_{x}^{2}}I_{N} \otimes A_{1} + \frac{1}{h_{y}^{2}}A_{1} \otimes I_{N}$

在上文中采取最简单粗暴的方式直接生成 $A$ , 并调用库函数直接对线性系统进行求解, 实际上这样的方法不适用于大规模问题,并且效率较低,本文介绍 Poisson 系统的一些常用的解法,这些解法适用于规模较大的一些问题.

简单迭代法

考虑线性系统:
$A x = b$
将上系统改写为如下等价形式
$x = (I - C A) x + C b$
对上述形式进行不动点迭代得如下迭代格式:
$x^{(k+1)} = (I - CA)x^{(k)} + Cb$
其中 $C$ 待定, 由不动点迭代原理可知, 如果能够找到某个范数 $\|\cdot\|$ , 使得: $\| I - CA\| < 1$ , 则上述迭代格式收敛,(并且 $\| I - CA\| 越小时,收敛速度越快.$ ) 因此可知引入 $C$ 的目的即通过找到合适的 $C$ , 使得上述迭代格式收敛, $C$ 的选择应当尽可能的简洁.( 显然,取 $C = A^{-1}$ 时, $\| I - CA\| = 0$ , 从而迭代效果最好, 但显然这是不现实的,因为求矩阵逆将比求解一个线性方程组复杂的多.)
通过上面的思路得到的迭代方法统称为简单迭代法 , 最经典的两种简答迭代法即: Jacobi 迭代, Gauss-Seidel 迭代 (下面简记为 GS 迭代) .

Jacobi 迭代

Jacobi 迭代原理

对于矩阵 $A$ , 如果满足对角线元素的量级相对于该行其它元素较大时(一般而言,严格主对角占优即可), 可以使用 $J a c o b i$ 迭代, $J a c o b i$ 迭代取上面中的 $C$ 为:
$C = D^{-1} :$
其中 $D : = d i a g (A)$ 表示由 $A$ 的对角元素构成的对角矩阵.
则 $J a c o b i$ 迭代方法迭代格式为:
$x^{(k+1)} = (I - D^{-1}A)x^{(k)} + D^{-1}b$
一般的,对迭代格式需要做一些修改, 记:
$A = L + D + U$
其中 $L$ 表示 $A$ 的下三角部分, $U$ 表示 $A$ 的上三角部分, $D$ 表示 $A$ 的对角线部分, 则容易知道, $J a c o b i$ 迭代格式可以改写为:
$x^{(k+1)} = -D^{-1}\left( L + U \right) x^{(k)} + D^{-1}b$
可以知道, 当 $A$ 具有较好的稀疏性时, 上述迭代能够以 $O (N)$ 的计算量完成, 反之, 当 $A$ 为满矩阵时, 则上述计算过程仍然有 $O(N^{2})$ 的复杂度, 如果是这样, $J a c o b i$ 迭代和 $L U$ 分解比起来就没有很明显的优势了, 但本文关注 $P o i s s o n$ 系统的求解 (甚至其它的 PDE有很大一部分情况下,离散出来的矩阵具有比较好的稀疏性).

由收敛性条件, 可知如果 : $\| D^{-1}\left( L + U\right) \| < 1$ 时, $J a c o b i$ 迭代收敛, 因此前面要求对角线元素量级较大.

Jacobi 迭代 Python 实现

import numpy as np	# 调用numpy包

def jacobi_iter(A, b, x0 = None, tol = 1e-14, itmax = 1000):
 
	L = np.tril(A, -1)	# 获取矩阵 A 下三角部分
	U = np.triu(A, 1)  # 获取矩阵 A 上三角部分
	D_inv = 1 / np.diag(A) # 获取矩阵 A 对角元的逆 (D_inv 是一维数组)
	B = - np.einsum('i,ij->ij',D_inv, L + U) # 创建迭代矩阵 B
	g = D_inv * b # 
	
	iter_count = 0	 	# 设置初始迭代次数为 0 
	iter_err = 1		# 设置初始迭代误差为 1
	iter_max = itmax  # 设置最大迭代次数
	iter_tol = tol  # 设置停机条件
	iter_val =  np.random.random(len(b)) if x0 is None else x0 # 如果没有赋初值,随机生成初值.
	
	while iter_err >= iter_tol and iter_count < iter_max: # 进入循环迭代
		iter_rec = iter_val  # 记录上一次迭代变量
		iter_val = B @ iter_val + g # 获取下一个迭代值
		iter_count += 1 # 增加迭代次数
		iter_err = np.max(np.abs(iter_val - iter_rec)) # 计算迭代误差
		print('iteration numbers = %d, iteration error = %e'%(iter_count, iter_err))
	
	# 退出循环, 处理结果并进行相应输出
	if iter_count == iter_max: 
		print("exceeding the maximum iteration numbers") # 打印提示,超出最大迭代次数
		return None
		
	print("iteration numbers: %d, error = %e"%(iter_count, iter_err)) 
	return iter_val  # 如果在上一步没有 return 表示成功找到解,打印最终迭代信息,并返回方程组的解.

用以下程序测试迭代法:

def iter_test():
    A = np.array([[10, 5, 1], [2, 10, 3], [3, 3, 10]], dtype=np.float64)
    b = np.array([3, 2, 6])
    x = jacobi_iter(A, b)
    print(x)
    x = np.linalg.inv(A) @ b
    print(x)

输出结果如下:

iteration numbers: 56, error = 8.992806e-15
[ 0.25150421 -0.00842359  0.52707581] 	# Jacobi 迭代的结果
[ 0.25150421 -0.00842359  0.52707581]   # 用逆矩阵乘以 b 得到的结果

GS 迭代

GS 迭代算法原理

$G S$ 迭代是对 $J a c o b i$ 迭代的改进, 考虑 $J a c o b i 迭代的标量形式 :$

$x_{i}^{(k+1)} = \dfrac{1}{a_{ii}} \left[ b_{i} - \sum_{j < i} a_{ii} x_{j}^{(k)} - \sum_{j > i}a_{ij}x_{j}^{(k)} \right]$
注意, 当按照 $1,2,\cdots,n$ 计算时, 当计算 $\sum_{j < i}a_{ii}x_{j}^{(k)}$ 时, 已经计算出了 $x_{1}^{(k+1)}, \cdots, x_{i-1}^{(k+1)},$ 因此这里可以用已经计算出的结果进行改进, 改进得到的就是 $G S$ 迭代法:

$x_{i}^{(k+1)} = \dfrac{1}{a_{ii}} \left[ b_{i} - \sum_{j < i} a_{ij} x_{j}^{\color{red}{(k+1)}} - \sum_{j > i}a_{ij}x_{j}^{(k)} \right], \quad i = 1,2, \cdots, n$

理论上,由表达式可知, $G S$ 迭代和 $J a c o b i$ 迭代的复杂度相差不大,并且直观上 $G S$ 迭代具有更好的收敛性.

下面来看一下 $G S$ 迭代的矩阵形式, 写出 $G S$ 迭代的矩阵形式如下:

$Dx^{(k+1)} = b - Lx^{(k+1)} - Ux^{(k)}$

整理得:

$D + L)x^{(k+1)} = -Ux^{(k)} + b$

注意 $D + L$ 为三角矩阵

注意编程时候不要用矩阵形式,不要用矩阵形式,不要用矩阵形式…不要对上三角矩阵求逆,直接用标量形式, 或者直接用回代方法求解三角方程组.

GS迭代 Python 实现:

import numpy as np
def gs_iter(A, b, x0 = None, tol = 1e-14, itmax = 1000):

	L = np.tril(A, -1)	# 获取矩阵 A 下三角部分
	U = np.triu(A, 1)  # 获取矩阵 A 上三角部分
	D = np.diag(np.diag(A))
	P = D + L 
	n = len(b)  # 计算 b 长度
	
	iter_count = 0	 	# 设置初始迭代次数为 0 
	iter_err = 1		# 设置初始迭代误差为 1
	iter_max = itmax  # 设置最大迭代次数
	iter_tol = tol  # 设置停机条件
	iter_val =  np.random.random(n) if x0 is None else x0 # 如果没有赋初值,随机生成初值.
	
	while iter_err >= iter_tol and iter_count < iter_max: # 进入循环迭代
		iter_rec = iter_val  # 记录上一次迭代变量
		 # 获取下一个迭代值
		g = -U @ iter_val + b
		x = np.zeros(n)		# 注意 Python 的拷贝机制, 这里不能直接对 iter_val 中的元素逐个修改.
		 # 利用循环进行回代来求解下三角矩阵方程组 Pxk+1 = g 
		for i in range(n): 
			x[i] = (g[i] - (P[i, :i] * x[:i]).sum()) / P[i, i]
		iter_val = x
		iter_count += 1 # 增加迭代次数
		iter_err = np.max(np.abs(iter_val - iter_rec)) # 计算迭代误差
		print('iteration numbers = %d, iteration error = %e'%(iter_count, iter_err))
	
	# 退出循环, 处理结果并进行相应输出
	if iter_count == iter_max: 
		print("exceeding the maximum iteration numbers") # 打印提示,超出最大迭代次数
		return None
		
	print("iteration numbers: %d, error = %e"%(iter_count, iter_err)) 
	return iter_val

同样,使用上面的测试脚本,得到结果如下:

iteration numbers: 22, error = 5.939693e-15
[ 0.25150421 -0.00842359  0.52707581]
[ 0.25150421 -0.00842359  0.52707581]

可以看到, GS 需要的迭代次数比 Jacobi 迭代少

显然, $G S$ 迭代收敛性要比 $J a c o b i$ 迭代好, 但是注意, $J a c o b i$ 迭代能够并行计算, 而 $G S$ 迭代需要求解前 $i - 1$ 个分量,逐个刷新 $x^{(k+1)}$ 的每个分量, 因此 $G S$ 迭代和 $J a c o b i$ 迭代各自都拥有着非常明显的优势, 因此这两种迭代都是非常重要的.