多重网格方法
基本思想
一般的迭代法是在一种固定的网格上进行迭代,当网格比较细时,计算量十分大。多重网格说的是,在计算细网格上的精确解时,其初值是比它粗一些网格上的精确解构造的,因而迭代次数少。当然,求较粗的网格上的精确解,它的初值就是有更粗一些的网格上的精确解所构造的,如此往复,直到一个相当粗的网格位置。多重网格法的求解过程是一个递归的过程。
算法过程
- 预平滑
选定一个初值,先采用一般方法,比如说权雅克比方法迭代一两次,求出一个近似解 u l u_l ul。 - 粗网格校正
1、计算剩余量 r l = f l − A l u l r_l = f_l - A_lu_l rl=fl−Alul。
2、限制剩余量 r l − 1 = R l , l − 1 r l r_{l-1}=R_{l,l-1}r_l rl−1=Rl,l−1rl,其中的 R R R为限制算子,它实质上是一个插值算子。
3、在粗网格上,求下列问题精确解:
A l − 1 ∗ e l − 1 = r l − 1 A_{l-1}* e_{l-1}= r_{l-1} Al−1∗el−1=rl−1
其中, A l − 1 A_{l-1} Al−1是在粗网格上离散后得到的矩阵(如果是有限元方法,就是刚度矩阵)。如果在这个粗网格上精确解不好求,继续放粗,递归调用整个过程。
4、将 e l − 1 e_{l-1} el−1延拓到粗网格上,以求解 e l e_l el, e l = P l − 1 , l e l − 1 e_l = P_{l-1,l}e_{l-1} el=Pl−1,lel−1,其中 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+1−1,表示划分出的点数(内点),即不包括边界。划分的网格数目为
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=hl−2∗tridiag(−1,2,1),这时,限制算子
R
l
,
l
−
1
R_{l,l-1}
Rl,l−1为
N
l
−
1
∗
N
l
N_{l-1}*N_l
Nl−1∗Nl的矩阵,磨光算子
P
l
−
1
,
l
P_{l-1,l}
Pl−1,l为
N
l
∗
N
l
−
1
N_l*N_{l-1}
Nl∗Nl−1的矩阵,他们长成形如下面这个样子:
据此,编写的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\b 求 w*(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,l−1(fl−Alul)中的 u l u_l ul使用预平滑之前的值,即初值,结果会更精确,我并不知道为什么。
3、matlab中semilogy函数,只不过是在刻度显示上,关于指数均匀增长,而y值是没有改变的。semilogy控制的只是坐标轴的剖分,按照指数等刻度均匀划分,并没有实际对y值取log。