点击刷新
- 算法特征:
①. pre-smoothing提取低频残差; ②. 向下插值计算残差补偿; ③.向上插值填充残差补偿; ④. post-smoothing降低整体残差 - 算法推导:
Part Ⅰ: 算法原理
考虑一般线性系统:
\begin{equation}
Ax = b
\label{eq_1}
\end{equation}
给定某初始值$x^{0}$, 残差为:
\begin{equation}
r^{0} = b - Ax^{0}
\label{eq_2}
\end{equation}
以该残差为基础, 计算残差补偿:
\begin{equation}
\delta x^{0} = A^{-1}r^{0}
\label{eq_3}
\end{equation}
填充残差补偿, 式(\ref{eq_1})解表示如下:
\begin{equation}
x = x^{0} + \delta x^{0}
\label{eq_4}
\end{equation}
Part Ⅱ: 误差(即残差补偿)特征
现以Jacobi迭代为例, 分解系数矩阵$A$, 式(\ref{eq_1})转换为
\begin{equation}
(L + D + U)x = b
\label{eq_5}
\end{equation}
其中, $L$、$D$、$U$分别为系数矩阵$A$的下三角矩阵、对角矩阵及上三角矩阵. Jacobi迭代公式表示如下:
\begin{equation}
x^{k+1} = Jx^k + D^{-1}b
\label{eq_6}
\end{equation}
其中, $J = -D^{-1}(L + U)$, 代表Jacobi迭代矩阵. 定义第$k$步迭代误差:
\begin{equation}
e^k = x^k - x
\label{eq_7}
\end{equation}
根据式(\ref{eq_6})可知, 迭代误差满足如下条件:
\begin{equation}
e^{k+1} = Je^k
\label{eq_8}
\end{equation}
结合上式, 对矩阵$J$进行特征分解, 则有
\begin{equation}
e^k = V\Lambda^kV^{-1}e^0
\label{eq_9}
\end{equation}
其中, 矩阵$V$和$\Lambda$分别由矩阵$J$的本征矢和本征值组成. 显然, Jacobi迭代的收敛性由相应迭代矩阵的本征值决定. 若矩阵$J$所有本征值绝对值均小于1, 则Jacobi迭代收敛; 否则, 不收敛.
现以有限差分法离散化1维Poisson's equation($\Delta u = f$)为例, 采用Jacobi迭代求解该线性系统, 并分析误差成分:
\begin{equation}
\begin{split}
A &= \begin{bmatrix}-2 & 1 & & \\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ & & 1 & -2 \end{bmatrix} \frac{1}{\Delta x^2} \\
D &= -2I\frac{1}{\Delta x^2} \\
L + U &= \begin{bmatrix} 0 & 1 & & \\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ & & 1 & 0 \end{bmatrix} \frac{1}{\Delta x^2} \\
J &= \begin{bmatrix} 0 & 1/2 & & \\ 1/2 & \ddots & \ddots & \\ & \ddots & \ddots & 1/2 \\ & & 1/2 & 0 \end{bmatrix}
\end{split}
\label{eq_10}
\end{equation}
矩阵$J$之本征矢与本征值分别为:
\begin{equation}
\begin{split}
V_{ij} &= \sin\frac{ij\pi}{n+1} \\
\lambda_j &= \cos\frac{j\pi}{n+1}
\end{split}
\label{eq_11}
\end{equation}
其中, $i$、$j$取$[1, n]$之间的整数. 显然, 矩阵$J$之本征值$\lambda \in (-1, 1)$, Jacobi迭代在该线性系统上收敛.
观察本征矢取值, 当$j$逐渐增大时, 频率逐渐升高. 由此可将误差大致分为低频组分与高频组分. 低频组分, 误差分布平滑, 利于近似表达; 高频组分, 误差分布不平滑, 不利于近似表达.
观察本征值取值, 当$j$逐渐增大时, 其变化趋势为: $1 \rightarrow 0 \rightarrow -1$. 低频部分, 误差收敛速度慢; 中频部分, 误差收敛速度快; 高频部分, 误差收敛速度慢.
对于低频部分, 误差分布平滑, 可采用向下插值技术将离散点上的残差由精细网格过渡至粗糙网格(即降低$n + 1$的取值), 以提升误差组分的频率, 达到快速收敛之目的, 再由向上插值将离散点上的误差由粗糙网格过渡至精细网格, 获取一个可以接受的残差补偿; 对于中频部分, 误差收敛速度快, 无需特殊处理; 对于高频部分, 由于误差收敛速度慢且分布不平滑, 需要引入under-relaxation技术.
Part Ⅲ: under-relaxation技术
结合式(\ref{eq_6}), 采用如下凸组合改写Jacobi迭代:
\begin{equation}
x^{k+1} = (1 - \alpha)x^k + \alpha [-D^{-1}(L + U)x^k + D^{-1}b], \quad \text{where }\alpha \in (0, 1)
\label{eq_12}
\end{equation}
根据式(\ref{eq_7})之定义, 误差方程转换如下:
\begin{equation}
e^{k+1} = [(1 - \alpha )I - \alpha D^{-1}(L + U)]e^k
\label{eq_13}
\end{equation}
此时, 迭代矩阵及其本征值分别为:
\begin{equation}
\begin{cases}
J_\alpha = (1 - \alpha)I - \alpha D^{-1}(L + U) \\
\lambda_\alpha = 1 - \alpha + \alpha\lambda
\end{cases}
\label{eq_14}
\end{equation}
根据式(\ref{eq_11}), 本征值具体取值为:
\begin{equation}
\lambda_j^\alpha = 1 - \alpha + \alpha \cdot \cos\frac{j\pi}{n+1} \in (1 - 2\alpha, 1)
\label{eq_15}
\end{equation}
若$\alpha = 0.5$, 则本征值取值区间为$(0, 1)$. 此时已成功消除分布不平滑且收敛速度慢的误差组分, 使Jacobi迭代在多重网格上求解Poisson's equation时快速收敛. - 代码实现:
Part Ⅰ:
现以Poisson's equation为例进行算法实施, 原始图像、Laplace变换图像及边界条件图像分别如下方左、中、右三图所示:
Part Ⅱ:
Multigrid实现如下:
View Code1 # 多重网格法求解Poisson's equation 2 3 import copy 4 import numpy 5 from PIL import Image 6 from matplotlib import pyplot as plt 7 from matplotlib import animation 8 9 10 11 def get_F(picName): 12 ''' 13 根据图片数据获取Poisson's equation等式右侧值 14 ''' 15 img = Image.open(picName) 16 arr = numpy.array(img).astype(float) 17 img.close() 18 19 dx = 1 / (arr.shape[1] - 1) 20 21 F = numpy.zeros(arr.shape) 22 for i in range(1, F.shape[0]-1): 23 for j in range(1, F.shape[1]-1): 24 F[i, j] = ((arr[i-1, j] + arr[i+1, j] + arr[i, j-1] + arr[i, j+1]) - 4 * arr[i, j]) / dx ** 2 25 return F 26 27 28 29 def get_B(picName): 30 ''' 31 根据图片数据获取Poisson's equation边界条件, 非边界值填充0 32 ''' 33 img = Image.open(picName) 34 B = numpy.array(img).astype(float) 35 img.close() 36 B[1:-1, 1:-1] = 0 37 return B 38 39 40 41 def jacobi(U, F): 42 ''' 43 Jacobi迭代法求解Poisson's equation: ΔU = F 44 U: 迭代解初始值, 含边界条件 45 F: Poisson's equation右侧值, 边界值无效 46 横轴长度固定取1 47 ''' 48 dx = 1 / (U.shape[1] - 1) 49 UNew = copy.deepcopy(U) 50 for i in range(1, UNew.shape[0]-1): 51 for j in range(1, UNew.shape[1]-1): 52 UNew[i, j] = (U[i-1, j] + U[i+1, j] + U[i, j-1] + U[i, j+1] - dx ** 2 * F[i, j]) / 4 53 return UNew 54 55 56 57 class Multigrid(object): 58 ''' 59 多重网格法之实现 60 ''' 61 def __init__(self, B, F, alpha=0.5): 62 self.__B = B # 边界条件, 矩阵非边界值无效 63 self.__F = F # Poisson's equation等式右侧 64 self.__alpha = alpha # under-relaxation凸组合因子 65 66 self.__UList = list() # 记录每帧迭代解 67 68 69 def get_solu(self): 70 return self.__UList 71 72 73 def optimize(self, maxIter=3000, epsilon=1.e-1): 74 ''' 75 maxIter: 最大迭代次数 76 epsilon: 收敛判据, 迭代解趋于稳定则收敛 77 ''' 78 self.__UList.clear() 79 80 U0 = copy.deepcopy(self.__B) 81 F0 = copy.deepcopy(self.__F) 82 alpha = self.__alpha 83 84 self.__UList.append(numpy.rint(U0).astype("int")) 85 for idx in range(maxIter): 86 R0 = self.__calc_residual(U0, F0) 87 R0Norm = numpy.linalg.norm(R0) 88 print("iterCnt: {:3d}, ResidualNorm: {}".format(idx, R0Norm)) 89 if R0Norm <= epsilon: 90 break 91 92 U0 = self.__do_multigrid(U0, F0, alpha) 93 self.__UList.append(numpy.rint(U0).astype("int")) 94 95 return U0 96 97 98 def __do_multigrid(self, U0, F0, alpha): 99 # pre-smoothing 100 Uf = self.__smoothing(U0, F0, alpha, 10) 101 Rf = self.__calc_residual(Uf, F0) # 残差计算 102 # downward interpolation - find grid to coarse grid 103 Rc = self.__downward_interpolate(Rf) # 向下插值 - 残差计算 104 deltaU0 = numpy.zeros(Rc.shape) 105 if (Rc.shape[0] >= 5 and Rc.shape[1] >= 5): 106 deltaUc = self.__do_multigrid(deltaU0, Rc, alpha) 107 else: 108 deltaUc = self.__smoothing(deltaU0, Rf, alpha, 10) 109 # upward interpolation - coarse grid to fine grid 110 deltaUf = self.__upward_interpolate(deltaUc, Uf) # 向上插值 - 残差补偿计算 111 Uf += deltaUf 112 # post-smoothing 113 U0 = self.__smoothing(Uf, F0, alpha, 10) 114 return U0 115 116 117 def __upward_interpolate(self, deltaUc, Uf): 118 deltaUf = numpy.zeros(Uf.shape) 119 deltaUf[2:-2:2, 2:-2:2] = deltaUc[1:-1, 1:-1] 120 deltaUf[1:-1:2, 2:-2:2] = (deltaUf[0:-2:2, 2:-2:2] + deltaUf[2::2, 2:-2:2]) / 2 121 deltaUf[:, 1:-1:2] = (deltaUf[:, 0:-2:2] + deltaUf[:, 2::2]) / 2 122 return deltaUf 123 124 125 def __downward_interpolate(self, Rf): 126 ''' 127 向下插值 - R中边界值无效 128 ''' 129 RcOri = Rf[2:-2:2, 2:-2:2] 130 RcTra = numpy.zeros((RcOri.shape[0]+2, RcOri.shape[1]+2)) 131 RcTra[1:-1, 1:-1] = RcOri 132 return RcTra 133 134 135 def __calc_residual(self, U, F): 136 ''' 137 残差计算 138 ''' 139 dx = 1 / (U.shape[1] - 1) 140 R = numpy.zeros(U.shape) 141 142 for i in range(1, U.shape[0]-1): 143 for j in range(1, U.shape[1]-1): 144 R[i, j] = F[i, j] - (U[i-1, j] + U[i+1, j] + U[i, j-1] + U[i, j+1] - 4 * U[i, j]) / dx ** 2 145 return R 146 147 148 def __smoothing(self, U, F, alpha, iterNum=5): 149 for idx in range(iterNum): 150 U = (1 - alpha) * U + alpha * jacobi(U, F) # under-relaxation 151 return U 152 153 154 155 class MultigridPlot(object): 156 157 fig = None 158 ax1 = None 159 line = None 160 txt = None 161 UList = None 162 163 @classmethod 164 def plot_ani(cls, mgObj): 165 mgObj.optimize() 166 167 cls.UList = mgObj.get_solu() 168 cls.fig = plt.figure(figsize=(5, 4)) 169 cls.ax1 = plt.subplot() 170 cls.line = cls.ax1.imshow(cls.UList[0], cmap="gray") 171 cls.txt = cls.ax1.text(x=100, y=100, s="") 172 173 ani = animation.FuncAnimation(cls.fig, cls.update, frames=numpy.arange(len(cls.UList)), blit=True, interval=1000, repeat=True) 174 175 cls.ax1.xaxis.set_visible(False) 176 cls.ax1.yaxis.set_visible(False) 177 cls.fig.tight_layout() 178 ani.save("plot_ani.gif", writer="PillowWriter", fps=3, dpi=1000) 179 plt.close() 180 181 182 @classmethod 183 def update(cls, frame): 184 cls.line.set_data(cls.UList[frame]) 185 cls.txt.set_text("iterCnt: {}".format(frame)) 186 return cls.line, cls.txt 187 188 189 @staticmethod 190 def plot_fig(mgObj): 191 U = mgObj.optimize() 192 img = Image.fromarray(U).convert("L") 193 img.save("plot_fig.png") 194 img.close() 195 196 197 198 if __name__ == "__main__": 199 picName = "./telescope-gray.jpg" 200 F = get_F(picName) 201 B = get_B(picName) 202 203 img = Image.fromarray(F).convert("L") 204 img.save("F.png") 205 img.close() 206 207 img2 = Image.fromarray(B).convert("L") 208 img2.save("B.png") 209 img2.close() 210 211 mgObj = Multigrid(B, F) 212 MultigridPlot.plot_fig(mgObj) 213 # MultigridPlot.plot_ani(mgObj)
-
结果展示:
显然, 通过Multigrid方法, 根据中间Laplace变换图像及右侧边界条件图像即可实现对左侧原始图像的快速求解, 迭代速度较常规迭代方法(Jacobi迭代、Gauss-Seidel迭代等)提升明显. -
使用建议:
①. 适用于结构化网格, 计算过程在精细网格与粗糙网格间多次变换;
②. 适用于求解椭圆型方程, 尤其是线性椭圆型方程. -
参考文档:
Numerical Methods for PDE 2016 by Professor Qiqi Wang & Jacob White from MIT
© 著作权归作者所有
举报
相关热门文章
发表评论
0/200
