matlab 多重循环在最外层加断点_多重网格法 (Multigrid Method) 简述和示例

1. 简介

多重网格法是迭代方法的一种，可用于求解线性方程组。在偏微分数值解上常用于求解稳态 (Steady-State) 和边界值问题 (Boundary Value Problem)。

不同于 Jacobi 迭代法和 Gauss-Seidel 迭代法，Multigrid Method 没有具体的公式形式，而是一整套求解的思路，求解过程中会利用到如 Jacobi 和 Gauss-Seidel 等基本的迭代方法。

多重网格法可以通过粗网格和细网格之间的转换有效地提高迭代算法的收敛速度。

2. Jacobi 和 Gauss-Seidel 迭代法的缺陷

迭代法求解线性方程组

得到

，误差

可被傅里叶分解为一组三角函数的积分，这些三角函数的频率 (wavenumber) 有大有小。对于任意固定大小的等距网格，Jacobi 和 Gauss-Seidel 迭代法消减频率较大的误差的速度较快，而消减频率较小的误差的速度往往十分缓慢。因此，在迭代过程中，频率较大的误差往往会率先消减掉，最后留下频率较小的误差。

• 为了更详细地说明这一点，观察下面这个边界值问题：

设置

个等距网格，设置初始值为

，利用 Gauss-Seidel 迭代法进行迭代，公式为：

不难计算，初始误差切好是微分方程的右边，两个频率不同的三角函数之和。

MATLAB 代码见附录。绘制误差随迭代次数变化的图像如下：

可以较为明显地看到，收敛行为分为两段，先收敛速度快，后收敛速度慢。这是因为同一等距网格在对付频率较高的误差时收敛速度较快，而在对付频率较低的误差时收敛速度较慢。

同样的，我们也可以绘制在不同迭代次数时误差在

轴上分布的图像：

由上图同样可以观察到，波动频率较大的误差随着迭代次数的增加率先衰减，而频率较小的误差则衰减得非常慢。

综上所述，面对频率小的误差时，Jacobi 和 Gauss-Seidel 迭代法的收敛速度极慢，消耗的计算资源较大。这一点虽然可以通过使用更粗的网格来克服（增大网格的尺度之后，相对于网格而言，误差的频率变大了），但如果遇到的误差中既有频率大的部分也有频率小的部分（如上例中的误差），那么加粗网格尽管提高了收敛速度，却使得精度大大下降了。

3. 多重网格法的基本思想

通过 2 中对 Jacobi 和 Gauss-Seidel 迭代法的观察，我们不难意识到交替使用粗细网格可以解决这些迭代方法的缺陷。先在细网格上迭代使得频率大的误差衰减，剩下误差中频率较小的部分；再将误差中频率较小的部分作为方程组的右侧，在粗网格上迭代，使得频率较小的误差也能快速地衰减；最后为了保证精度，将粗网格上迭代的结果和细网格上的结果加总之后，作为初始猜测值再在细网格上迭代。如此循环往复。因为先由细到粗，后由粗到细，网格转换的路径形似 V 字，所以该方法被称为 V-Cycle Multigrid。除此之外还有 F-cycle multigrid 和 W-cycle multigrid 。他们的基本思想都是相同。

4. 多重网格法的步骤

以最简单的 V-Cycle 为例

1. 设置合理的初始猜测值，在网格间距为
   的细网格上对
   进行
   次迭代，得到
   ，计算误差
2. 将误差
   转换为网格间距为
   的粗网格上的误差
   （restriction），在粗网格上对
   进行
   次迭代，得到
   ，计算误差
3. 用粗网格上的
   插值得到细网格上的
   ，将
   作为初始猜测值迭代求解
   ，迭代
   次之后得到新的
   ，计算新的误差
   。如果新的误差
   小于设定的阈值，则
   即为所求；否则转 2

以上内容为叙述方便只涉及 2 层网格，实际中往往使用更多重的网格

不难发现，1-2 是由细到粗，2-3 是由粗到细，要实现更多重的网格例如

重，每个 V-Cycle 需要进行

次 1-2，再进行

次 2-3。

5. 多重网格法示例

• 用多重网格法求解 2 中的边界值问题：

最细的一层设置

个等距网格，同样设置初始值为

，利用 Multigrid Method 进行迭代，每步的迭代次数

。参考了 Wikipedia 上的伪代码，MATLAB 代码见附录。

每运行一遍 V Cycle 会执行 Gauss-Seidel 迭代 10 次，运次附录中的程序显示为达到收敛条件（

）一共执行了 5 遍 V Cycle，所以多重网格法一共执行了 50 次 Gauss-Seidel 迭代。

对比 Gauss-Seidel 迭代法，绘制误差随迭代次数变化的图像如下：

由上图可直观地得知，Multigrid Acceleration 比 Gauss-Seidel Method 要高效得多，迭代次数更少，计算开销也更小。

Appendix (MATLAB Codes)

2. Gauss-Seidel Method

close all
clear all

N = 64; L = 1;
dx = L/N;

phi = zeros(1,N+1);
new = zeros(1,N+1);
tmp = (sin(pi*[1:N-1]*dx)+sin(16*pi*[1:N-1]*dx))/2;
r0 = zeros(1,N+1); r10 = zeros(1,N+1); r100 = zeros(1,N+1);
r0(2:N) = tmp;
resi(1) = max(abs(tmp));

for t = 1:100
    for j = 2:N
        new(j) = (phi(j+1)+new(j-1)-dx^2*tmp(j-1))/2;
    end
    new(1) = 0; new(N+1) = 0;
    r = tmp-(new(1:N-1)-2*new(2:N)+new(3:N+1))/dx^2;
    resi(t+1) = max(abs(r));
    phi = new;
    if t == 10
        r10(2:N) = r;
    elseif t == 100
        r100(2:N) = r;
    end
end

figure
plot([0:length(resi)-1],resi,'+-');
xlabel('Number of Iterations')
ylabel('max(|r_j|)')
title('Convergence Curve')

x = linspace(0,1,N+1);
figure
plot(x,r0,'-',x,r10,'+-',x,r100,'x-')
legend('0 iterations','10 iterations','100 iterations')
xlabel('x_j')
ylabel('r_j')
title('r_j against x_j')

5. Multigrid Method

close all
clear all

N = 64; L = 1;
h = L/N;
phi = zeros(1,N+1);
f = (sin(pi*[0:N]*h)+sin(16*pi*[0:N]*h))/2;

for cnt = 1:1000
    phi = V_Cycle(phi,f,h);
    r = residual(phi,f,h);
    if max(abs(r)) < 0.001
        break
    end
end

function phi = V_Cycle(phi,f,h)
 % Recursive V-Cycle Multigrid for solving the Poisson equation (nabla^2 phi = f) on a uniform grid of spacing h

 % Pre-Smoothing
 phi = smoothing(phi,f,h);
 
 % Compute Residual Errors
 r = residual(phi,f,h);
 
 % Restriction
 rhs = restriction(r);

 eps = zeros(size(rhs));

 % stop recursion at smallest grid size, otherwise continue recursion
 if length(eps)-1 == 2
         eps = smoothing(eps,rhs,2*h);
 else        
         eps = V_Cycle(eps,rhs,2*h);        
 end
 
 % Prolongation and Correction
 phi = phi + prolongation(eps);
 
 % Post-Smoothing
 phi = smoothing(phi,f,h);    
end

function res = smoothing(phi,f,h)
    N = length(phi)-1;
    res = zeros(1,N+1);
    for j = 2:N
        res(j) = (phi(j+1)+res(j-1)-h^2*f(j))/2;
    end
end

function res = residual(phi,f,h)
    N = length(phi)-1;
    res = zeros(1,N+1);
    res(2:N) = f(2:N)-(phi(1:N-1)-2*phi(2:N)+phi(3:N+1))/h^2;
end

function res = restriction(r)
    N = (length(r)-1)/2;
    res = zeros(1,N+1);
    for j = 2:N
        res(j) = (r(2*j-2)+2*r(2*j-1)+r(2*j))/4;
    end
end

function res = prolongation(eps)
    N = (length(eps)-1)*2;
    res = zeros(1,N+1);
    for j = 2:2:N
        res(j) = (eps(j/2)+eps(j/2+1))/2;
    end
    for j = 1:2:N+1
        res(j) = eps((j+1)/2);
    end
end