多重网格方法

多重网格方法

基本思想

一般的迭代法是在一种固定的网格上进行迭代,当网格比较细时,计算量十分大。多重网格说的是,在计算细网格上的精确解时,其初值是比它粗一些网格上的精确解构造的,因而迭代次数少。当然,求较粗的网格上的精确解,它的初值就是有更粗一些的网格上的精确解所构造的,如此往复,直到一个相当粗的网格位置。多重网格法的求解过程是一个递归的过程。

算法过程

  • 预平滑
    选定一个初值,先采用一般方法,比如说权雅克比方法迭代一两次,求出一个近似解 u l u_l ul
  • 粗网格校正
    1、计算剩余量 r l = f l − A l u l r_l = f_l - A_lu_l rl=flAlul
    2、限制剩余量 r l − 1 = R l , l − 1 r l r_{l-1}=R_{l,l-1}r_l rl1=Rl,l1rl,其中的 R R R为限制算子,它实质上是一个插值算子。
    3、在粗网格上,求下列问题精确解:
    A l − 1 ∗ e l − 1 = r l − 1 A_{l-1}* e_{l-1}= r_{l-1} Al1el1=rl1
    其中, A l − 1 A_{l-1} Al1是在粗网格上离散后得到的矩阵(如果是有限元方法,就是刚度矩阵)。如果在这个粗网格上精确解不好求,继续放粗,递归调用整个过程。
    4、将 e l − 1 e_{l-1} el1延拓到粗网格上,以求解 e l e_l el e l = P l − 1 , l e l − 1 e_l = P_{l-1,l}e_{l-1} el=Pl1,lel1,其中 P P P为延拓算子,其实它是一个插值算子。
    5、计算粗网格上的近似解,仍然用 u l ul ul表示:
    u l = u l + e l u_l = u_l + e_l ul=ul+el
  • 后平滑
    以算得的 u l u_l ul为初值,采用一般迭代法再迭代个一两次,求解近似解。

一个简单的一维多重网格算法流程如下所示,这里的一般迭代方法选用的Weighted Jacobi方法。

在这里插入图片描述

数值实验与结果

以一维举例,取 N l = 2 l + 1 − 1 N_l = 2^{l+1}-1 Nl=2l+11,表示划分出的点数(内点),即不包括边界。划分的网格数目为 2 l + 1 2^{l+1} 2l+1 ,则网格大小为 h l = 1 / 2 l + 1 h_l = 1/2^{l+1} hl=1/2l+1,取 A l = h l − 2 ∗ tridiag ( − 1 , 2 , 1 ) A_l = h_l^{-2}*\text{tridiag}(-1,2,1) Al=hl2tridiag(1,2,1),这时,限制算子 R l , l − 1 R_{l,l-1} Rl,l1 N l − 1 ∗ N l N_{l-1}*N_l Nl1Nl的矩阵,磨光算子 P l − 1 , l P_{l-1,l} Pl1,l N l ∗ N l − 1 N_l*N_{l-1} NlNl1的矩阵,他们长成形如下面这个样子:
在这里插入图片描述

据此,编写的matlab代码如下,使用的是环境是Matlab 2014a(盗版)。为了方便,我这里将主函数和函数放在一块了,实际使用时需要拆开。

clc
clear
l = 5;
fl = fl_gen(l);
N = 2^(l+1) - 1;
u0 = ones(N,1);
u_true = Al_gen(l)\fl;
m = 8;%算法迭代次数
ul = u0;
error_MG = norm((u_true - ul),2)/N;
for i=1:m
    ul = MG(l,fl,ul);
    error_MG(end+1) = norm((u_true - ul),2)/N;
end;

mm = 1000;
method = 'weighted_jacobi';%Richardson weighted_jacobi Gauss_Seidel SOR
w = 1/2;
ul = u0;
error_si_sol = norm((u_true - ul),2)/N;
for i = 1:m
    for j = 1:mm
        options = struct('A',Al_gen(l),'u_old',ul,'f',fl,'method',(method),'w',w);
        ul = simple_iter_solvers(options);
    end
    error_si_sol(end+1) = norm((u_true - ul),2)/N;
end

p1 = semilogy(error_MG);
hold on;
p2 = semilogy(error_si_sol);
legend('MG','weighted\_jacobi')
set(p1,'MarkerFaceColor',[0 0 1],'MarkerEdgeColor',[0 1 0],...
    'MarkerSize',10,...
    'Marker','o',...
    'LineStyle','--',...
    'Color',[1 0 0],...
    'DisplayName','MG');
set(p2,'MarkerFaceColor',[0 1 0],'MarkerEdgeColor',[1 0 0],...
    'MarkerSize',10,...
    'Marker','o',...
    'LineStyle','--',...
    'Color',[1 0 1]);

xlabel('迭代步数');
ylabel('误差');
title('MG和一般迭代法随迭代步增加收敛性质比较');


function ul = MG(l,fl,ul)
Al = Al_gen(l);
%rl = fl-Al*ul;
%if norm(rl,2)/norm(fl,2) < 10e-6 || l<1
if l == 1
  %  ul = zeros(2^(l+1)-1,1);
   % fl = fl_gen(1);
    ul = Al\fl;
    return;
end
%% 迭代方法的选择
method = 'weighted_jacobi';%可以改成其他方法Gauss_Seidel weighted_jacobi

%% 预平滑
%for i = 1:4
ul0 = ul;
options = struct('A',Al,'u_old',ul,'f',fl,'method',method,'w',1/2);
ul = simple_iter_solvers(options);
%end
%% 限制
[Pl,Rl] = PR_gen(l);
rl_1 = Rl*(fl-Al*ul0);
%rl_1 = Rl*(fl-Al*ul);
%% 粗网格校正
el_1 = MG(l-1,rl_1,zeros(2^l-1,1));
%% 延拓
ul = ul + Pl*el_1;
%% 后平滑
%for i = 1:4
options = struct('A',Al,'u_old',ul,'f',fl,'method',method,'w',1/2);
ul = simple_iter_solvers(options);
%end
end


function u_new = simple_iter_solvers(options)
%帮助信息:
%输入参数为options = struct('A',A,'u_old',u_old,'f',f,'method','method','w',w);
if nargin < 1, help(mfilename),
    printf('输入参数错误。')
end
%options = varargin;
A = options.A;
u_old = options.u_old;
f = options.f;
method = options.method;
D = diag(diag(A));
L = tril(A,-1);
U = triu(A,1);

if strcmp(method,'Richardson')==1
    w = options.w;
    u_new = u_old + w*(f-A*u_old);
end
if strcmp(method,'weighted_jacobi')==1
    w = options.w;
    %   u_new = u_old + w*D^-1*(f-A*u_old);
    u_new = u_old + w*(D\(f-A*u_old));
    %    u_new = u_new/4;
end
if strcmp(method,'Gauss_Seidel')==1
    u_new = u_old + (D+L)\(f-A*u_old);
end
if strcmp(method,'SOR')==1
    w = options.w;
    u_new = (D+w*L)\(w*f - (w*U+(w-1)*D)*u_old);
end

end



function [P,R] = PR_gen(l)
% 我们用PR_gen来生成限制算子R和磨光算子P,输入参数l表示较细网格划分了l次
N = 2^(l+1) - 1;
%h = 1/2^(l+1);
N_ = 2^l- 1;
%h_ = 1/2^(l);
p = 1/2*[1;2;1;zeros(N-3,1)];
P = zeros(N,N_);
for i = 0:N_-1
    P(:,i+1) = P(:,i+1) + circshift(p,2*i);
end
R = 1/2*P';
end


function fl = fl_gen(l)
NN = 2^(l+1);
h = 1/NN;
x = h:h:1-h;
fl = fun(x);
end
function fx = fun(x)
fx = sin(2*pi*x);
fx = fx';
%fx = ones(length(x),1);
%fx = zeros(length(x),1);
end



function A = Al_gen(l)
%% Al生成器
N = 2^(l+1) - 1;
h = 1/2^(l+1);
A = (1/h^2)*(diag(repmat(-1,N-1,1),-1)+diag(repmat(2,N,1))+diag(repmat(-1,N-1,1),1));
end

README:程序使用方法为打开main函数,设置好参数(层数l,一般迭代方法method、迭代步数m,初值u0等),直接运行。f函数在fl_gen子函数中设置,MG所用的一般迭代方法的选择在MG子函数中设置,可供选择的方法有Richardson、weighted_jacobi、Gauss_Seidel、SOR等。

我们想看到MG方法的收敛速度远快与一般迭代法,这里设置 l = 5 l=5 l=5,MG方法一直迭代到第1层,在第一层求解精确解。一般迭代法在程序中作为参数是可选的,我这里选的weighted jacobi方法。为了比较出效果和差距,这里的一般方法,我让其迭代1000次,算作一轮。MG算法和一般方法都跌打8轮,结果如下。

在这里插入图片描述

从图中可以看出,MG方法的误差,随着迭代步数的增加,是呈现指数级别降低的(这里plot对y坐标做了log处理,用了semilogy函数画图)。也就是说,一般迭代法,想要达到MG方法相同的精度,要做大几个量级数目的迭代。

打开matlab的探查,我们来看一下时间消耗。

在这里插入图片描述

从这可以看到,一般迭代法所用的时间消耗是MG方法的20多倍,但是从前面那个图上可以看到,在这种情况下,一般方法的收敛得依然没有MG方法快。当然,这只是一个大概的估计,因为在MG中也用了一般跌打法,且MG中也用到了 A l A_l Al生成函数,这是比较费时间的。由此,可以得到的一个结论是,在相同时间消耗的情况下,MG方法达到的精度是一般方法无法企及的。

其它一些东西

1、之前程序问题了,花了好长时间查错。后来发现一直用w*A\bw*(A^-1*b),这是不对的。因为w*A\b等价于(w*A)\b,求逆符号\和乘号*是同级的,是按从左往右算,一个一个括号的问题,耗费了我大量时间。

2、在求Restriction步骤时,我发现 R l , l − 1 ( f l − A l u l ) R_{l,l-1}(f_l-A_lu_l) Rl,l1(flAlul)中的 u l u_l ul使用预平滑之前的值,即初值,结果会更精确,我并不知道为什么。

3、matlab中semilogy函数,只不过是在刻度显示上,关于指数均匀增长,而y值是没有改变的。semilogy控制的只是坐标轴的剖分,按照指数等刻度均匀划分,并没有实际对y值取log。

©️2020 CSDN 皮肤主题: 代码科技 设计师:Amelia_0503 返回首页