多重網格方法

阿新 • • 發佈：2018-10-31

多重網格方法

基本思想

一般的迭代法是在一種固定的網格上進行迭代，當網格比較細時，計算量十分大。多重網格說的是，在計算細網格上的精確解時，其初值是比它粗一些網格上的精確解構造的，因而迭代次數少。當然，求較粗的網格上的精確解，它的初值就是有更粗一些的網格上的精確解所構造的，如此往復，直到一個相當粗的網格位置。多重網格法的求解過程是一個遞迴的過程。

演算法過程

預平滑
選定一個初值，先採用一般方法，比如說權雅克比方法迭代一兩次，求出一個近似解 $u_l$ 。
粗網格校正
1、計算剩餘量 $r_l = f_l - A_lu_l$ 。
2、限制剩餘量 $r_{l-1}=R_{l,l-1}r_l$ ，其中的 $R$ 為限制運算元，它實質上是一個插值運算元。
3、在粗網格上，求下列問題精確解：
$A_{l-1}* e_{l-1}= r_{l-1}$
其中， $A_{l-1}$ 是在粗網格上離散後得到的矩陣（如果是有限元方法，就是剛度矩陣）。如果在這個粗網格上精確解不好求，繼續放粗，遞迴呼叫整個過程。
4、將 $e_{l-1}$ 延拓到粗網格上，以求解 $e_l$ ， $e_l = P_{l-1,l}e_{l-1}$ ，其中 $P$ 為延拓運算元，其實它是一個插值運算元。
5、計算粗網格上的近似解，仍然用 $ul$ 表示：
$u_l = u_l + e_l$
後平滑
以算得的 $u_l$ 為初值，採用一般迭代法再迭代個一兩次，求解近似解。

一個簡單的一維多重網格演算法流程如下所示，這裡的一般迭代方法選用的Weighted Jacobi方法。

在這裡插入圖片描述

數值實驗與結果

以一維舉例，取 $N_l = 2^{l+1}-1$ ，表示劃分出的點數（內點），即不包括邊界。劃分的網格數目為 $2^{l+1}$ ，則網格大小為 $h_l = 1/2^{l+1}$ ，取 $A_l = h_l^{-2}*\text{tridiag}(-1,2,1)$ ，這時，限制運算元 $R_{l,l-1}$ 為 $N_{l-1}*N_l$ 的矩陣，磨光運算元 $P_{l-1,l}$ 為 $N_l*N_{l-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$ ，MG方法一直迭代到第1層，在第一層求解精確解。一般迭代法在程式中作為引數是可選的，我這裡選的weighted jacobi方法。為了比較出效果和差距，這裡的一般方法，我讓其迭代1000次，算作一輪。MG演算法和一般方法都跌打8輪，結果如下。

在這裡插入圖片描述

從圖中可以看出，MG方法的誤差，隨著迭代步數的增加，是呈現指數級別降低的（這裡plot對y座標做了log處理，用了semilogy函式畫圖）。也就是說，一般迭代法，想要達到MG方法相同的精度，要做大幾個量級數目的迭代。

開啟matlab的探查，我們來看一下時間消耗。

在這裡插入圖片描述

從這可以看到，一般迭代法所用的時間消耗是MG方法的20多倍，但是從前面那個圖上可以看到，在這種情況下，一般方法的收斂得依然沒有MG方法快。當然，這只是一個大概的估計，因為在MG中也用了一般跌打法，且MG中也用到了 $A_l$ 生成函式，這是比較費時間的。由此，可以得到的一個結論是，在相同時間消耗的情況下，MG方法達到的精度是一般方法無法企及的。

其它一些東西

1、之前程式問題了，花了好長時間查錯。後來發現一直用w*A\b 求 w*(A^-1*b)，這是不對的。因為w*A\b等價於(w*A)\b，求逆符號\和乘號*是同級的，是按從左往右算，一個一個括號的問題，耗費了我大量時間。

2、在求Restriction步驟時，我發現 $R_{l,l-1}(f_l-A_lu_l)$ 中的 $u_l$ 使用預平滑之前的值，即初值，結果會更精確，我並不知道為什麼。

3、matlab中semilogy函式，只不過是在刻度顯示上，關於指數均勻增長，而y值是沒有改變的。semilogy控制的只是座標軸的剖分，按照指數等刻度均勻劃分，並沒有實際對y值取log。

多重網格方法