多重网格法的Matlab实现

多重网格法的Matlab实现

  在我的上一篇博客两网格法的Matlab实现中，首先介绍了两网格法的详细算法流程以及两网格法的核心组件：提升算子以及限制算子，接着通过Matlab编程实现了两网格法对一个简单的偏微分方程问题的近似计算，最后得出了两网格法可以提升计算效率，但是迭代步数往往较高的结论，为了解决迭代步数较高的问题，我们可以尝试使用多重网格法。

多重网格法简介

  多重网格法（Multigrid Method）是一种高效的数值算法，用于求解大规模线性和非线性偏微分方程。该方法通过利用不同精细程度的网格层次来加速求解过程，特别适合于那些直接方法或标准迭代方法难以高效处理的问题。多重网格法的核心思想是在多个尺度上同时处理误差，从而实现快速收敛。

多重网格法工作原理

  多重网格法基于这样一个观察：传统的迭代方法（如雅可比或高斯-赛德尔方法）倾向于快速消除误差的高频部分，而低频部分则消除得较慢。多重网格法通过以下步骤来解决这个问题：

1.平滑操作（Smoothing）：在细网格上使用传统迭代方法几步，以快速消除高频误差。
2.粗化（Coarsening）：将问题从细网格转移到更粗的网格上。这通常通过一个称为限制算子的过程完成，它将细网格上的残差映射到粗网格上。
3.求解或递归：在粗网格上求解或递归使用多重网格方法来处理问题，这有助于处理那些在细网格上处理不了的低频误差部分。
4.细化（Refinement）或插值（Prolongation）：将粗网格上的解或误差校正项传回细网格。使用提升算子（也称为插值算子）来实现，这个步骤有助于改善细网格上的解。
5.后平滑：再次在细网格上应用几步平滑操作，以进一步提炼解。

  多重网格方法可以设计为两层（两网格法）或多层的形式。在多层形式中，问题被递归地传递到越来越粗的网格上，每个层次都可以使用更粗层次的网格进行校正。最粗的网格通常足够小，可以直接求解，其解可以逐级传回到最细的网格层次。

多重网格法应用，优点及挑战

  多重网格法非常适合于大规模稀疏矩阵问题，如那些来源于偏微分方程的离散化。它的优点包括：

1.高效性：与传统的线性求解器相比，多重网格方法可以更快地达到解的精度，计算复杂度通常与问题大小成线性关系。
2.灵活性：适用于多种类型的方程，包括椭圆型、抛物型和双曲型偏微分方程。
3.鲁棒性：对不同类型的问题和边界条件都能保持良好的性能。
  多重网格方法的主要挑战在于正确设计各层次网格的结构、选择合适的平滑器、限制和提升算子，这些设计决策对算法的总体效率和效果有决定性影响。如果想要对多重网格法有更加细致的了解，可以参照陈龙教授的UCI，里面对两网格法，多重网格法以及其它的有限元方法都有详细讲解。

多重网格法的Matlab实现

  在了解了多重网格法的工作原理以及优缺点后，接下来我们对两网格法的Matlab实现中的问题使用多重网格法进行解决，具体的Matlab代码如下：

% 定义Gauss_Seidel迭代函数
function x=Gauss_Seidel(A,u0,b,num)
    D=diag(diag(A));
    L=D-tril(A);
    U=D-triu(A);
    M=(D-L)\U;
    g=(D-L)\b;
    for i=1:num
        u=M*u0+g;
        u0=u;
    end
    x=u;
end   

clc,clear 
% 清空之前的变量，以便对此程序产生不必要的干扰

format long;

H = [1/16,1/32,1/64,1/128];
% 将h的可能值用集合表示，方便接下来的使用

v1 = 1;
% 限制过程中Gauss_Seidel迭代的次数

for i =1:4
% 最外层循环，由于H的可能值为四个，因此进行四次循环
    h = H(i);
    N = 1/h;
    A = diag(2/h*ones(N-1,1),0)-diag(1/h*ones(N-2,1),1)-diag(1/h*ones(N-2,1),-1);
    % 由有限元方法的矩阵表示可以得出Ah的具体形式，此处记为A

    f = zeros(N-1,1);
    for j =1:N-1
        f(j) = pi*pi*sin(pi*j*h); % 右端项
    end
    B = diag(h/6*ones(N-2,1),1)+diag(h/6*ones(N-2,1),-1)+diag(2*h/3*ones(N-1,1),0);
    f1 = B*f;
    % 此处为求解有限元问题的右端项fh，可以对其积分进行近似求解

    u0 = ones(N-1,1);
    uh = Gauss_Seidel(A,u0,f,v1);
    % 定义迭代初值，并对其进行v1次Gauss_Seidel迭代得到uh

    rh = f1 - A*uh;
    % 求出迭代误差

    l = i+2;
    % 定义内层循环的次数，由于最粗网格的h为1/2，因此分别进行限制和提升（i+2)次

    C = cell(10,1);
    C{1} = f1;
    % C表示一个元组，用于存储右端项矩阵

    D = cell(10,1);
    D{1} = A;
    % D表示另一个元组，用于存储系数矩阵

    E = cell(10,1);
    E{1} = uh;
    % E表示另一个元组，用于存储迭代后的Uh

    count = 0;
    %count表示计数器，用于记循环次数

    while norm (rh,2)>1e-8
        for k =1:l
            I2h_h=zeros(N/2^k-1,N/2^(k-1)-1);
            Ih_2h=zeros(N/2^(k-1)-1,N/2^k-1);
            for n=1:1:N/2^k-1
                I2h_h(n,2*n-1)=1/4;
                I2h_h(n,2*n)=2/4;
                I2h_h(n,2*n+1)=1/4;
                Ih_2h(2*n-1,n)=1/2;
                Ih_2h(2*n,n)=1;
                Ih_2h(2*n+1,n)=1/2;
            end
            % 求出提升算子与限制算子

            if k<l
               f2 = I2h_h*rh;
               C{k+1} = f2;
               u1 = zeros(N/2^k-1,1);
               A1 = I2h_h*A*Ih_2h;
               A = A1;
               D{k+1} = A1;
               u1h = Gauss_Seidel(A1,u1,f2,v1);
               E{k+1} = u1h;
               r1h = f2-A1*u1h;
               rh = r1h;
            end
            % 当h<1/2时，按照算法步骤对其进行限制，通过Jacobi迭代求出uh的近似值

            if k == l
               f2 = I2h_h*rh;
               C{k+1} = f2;
               A1 = I2h_h*A*Ih_2h;
               A = A1;
               D{k+1} = A1;
               u1h = A1\f2;
               E{k+1} = u1h;
               r1h = f2-A1*u1h;
               rh = r1h;
            end
            % 当h=1/2时，对uh进行精确求解
            
       end

       v2 = 1;
       % 定义V2，用于表示接下来每次Jacobi迭代的次数

       A = diag(2/h*ones(N-1,1),0)-diag(1/h*ones(N-2,1),1)-diag(1/h*ones(N-2,1),-1);
       % 由于上面的循环导致A已经发生改变，而后面的求解又要使用初始的A，因此对A重新赋值

       for m =1:l
           I2h_h=zeros(N/2^(l-m+1)-1,N/2^(l-m)-1);
           Ih_2h=zeros(N/2^(l-m)-1,N/2^(l-m+1)-1);
           for t=1:1:N/2^(l-m+1)-1
               I2h_h(t,2*t-1)=1/4;
               I2h_h(t,2*t)=2/4;
               I2h_h(t,2*t+1)=1/4;
               Ih_2h(2*t-1,t)=1/2;
               Ih_2h(2*t,t)=1;
               Ih_2h(2*t+1,t)=1/2;
           end
           u2h = E{l-m+1}+Ih_2h*E{l-m+2};
           u3h = Gauss_Seidel(D{l-m+1},u2h,C{l-m+1},v2);
           E{l-m+1} = u3h;
           % 提升算法，将步长从1/2提升至h

       end
       rh = f1-A*u3h;
       % 再次求出误差，并对其进行判断是否小于1e-8，若不小于，则继续进行限制与提升，直到满足为止

       count = count+1;
       % 通过计数器记录其循环次数

    end
    fprintf('步长为%d/%d 迭代次数为 %d 此时误差为 %e \n',1,N,count,norm(rh));
    % 输出步长对应的迭代次数，以及最后一次循环之后的误差

    u=[0;u3h;0];
    x=0:h:1;
    plot(x,u,"LineWidth",1.5);
    % 做出误差和横坐标x的图像

    hold on;
    % 将图像在一个框里面输出

end

legend("N=16","N=32","N=64","N=128")
xlabel("x")
ylabel("u")
% 对图像横坐标，纵坐标，以及N进行标注

运行程序，可以得到如下的结果：

  由以上的结果可以得出以下两个重要结论：首先是随着步长的减小，迭代的次数逐渐增加，其次根据之前编过的两网格法结果可知：多重网格方法的迭代次数要远少于两网格法的迭代次数，也就是说多重网格法可以更快的使误差小于1e-8。同时我们可以得到多重网格法近似的结果图像为：

  由上图可以看出只需要迭代十几步便可以得到与原函数$sin(\pi x)$特别相近的图像，由此可见多重网格法的强大！但是多重网格法的效果很大程度上依赖于对网格层次、平滑器、限制算子和提升算子的正确选择和配置。对于不同的问题，找到最佳的组合可能需要大量的试验和错误，这可能导致初始设置和调试过程变得复杂和耗时。同时多重网格法理论比较复杂，搞懂理论便要花费很长的时间，因此在实际应用中有很多限制。