Poisson方程的五點差分格式例題求解-Matlab實作

1.前言

本文使用Matalab軟體，將應用偏微分方程數值解中橢圓形方程的五點差分格式求解一道Poisson方程的第一邊值例題，通過五點差分格式及邊值條件得到相應的差分方程組 K U = F KU=F KU=F后，采用Gauss-seidel迭代法對其求解即可得到數值解，并將數值解與真解作比較.
其中五點差分格式將直接給出，其構造程序從略.

2.例題

構造Poisson方程第一邊值如下：

采用五點差分格式求解并與決議解 u ( x , y ) = x 2 y 2 u(x,y)=x^2y^2 u(x,y)=x2y2進行比較.
( x 和 y x和y x和y方向均 n n n等分，取 h 1 = h 2 = h = 1 / n h_1=h_2=h=1/n h1?=h2?=h=1/n， x 和 y x和y x和y方向上的節點序號分別用 i , j i,j i,j表示)

3.符號說明

表1：

4.思路

主要思路是將二維問題當成一維問題求解，即二維的 ( n + 1 ) 2 (n+1)^2 (n+1)2個節點按順序分行放到 ( n + 1 ) 2 (n+1)^2 (n+1)2維向量 U U U中.

主要步驟如下：

1.確定步長后對區域進行網格均勻剖分（ x , y x,y x,y方向的步長 h h h相等）；

2.所求方程組為 K U = F KU=F KU=F，根據五點差分格式分別給 K K K和 F F F賦值（矩陣中包含邊界點）；

(說明： U U U是一個 ( n + 1 ) 2 (n+1)^2 (n+1)2維向量，即將求解區域中所有節點放到一個向量 U U U中，當作一維問題進行求解；其次，將邊界點放入方程組中可以使 K K K具有更好的性質，更好地防止求解中不收斂的情況，處理邊界點時將其當成未知的即可)

3.由邊界點值的算術平均值作為 U U U的初始值,以減少迭代次數；

4.給定迭代次數上限 N N N和誤差限，由 G a u s s ? s e i d e l Gauss-seidel Gauss?seidel迭代求解方程組 K U = F KU=F KU=F得到數值解向量 u u u；

5.為方便觀察，將解向量 u u u的元素放入 ( n + 1 ) ? ( n + 1 ) (n+1)*(n+1) (n+1)?(n+1)階矩陣 U 1 U_1 U1?中，并與由決議解 u ( x , y ) = x 2 y 2 u(x,y)=x^2y^2 u(x,y)=x2y2產生的解矩陣 U 2 U_2 U2?比較；

5.求解步驟

劃分網格數 n n n選取為7.

1.右端項為 f i j = ? 2 [ ( i h ) 2 + ( j h ) 2 ] f_{ij}=-2[(ih)^2+(jh)^2] fij?=?2[(ih)2+(jh)2]，則得到方程的五點差分格式為

相應的 G a u s s ? s e i d e l Gauss-seidel Gauss?seidel迭代公式為

2.構造求解矩陣時，對特殊型別的節點應根據邊界條件將格式適當化簡：
（1）左邊界點 u 0 j ( j = 0 , 1 , . . . , n ) : 4 u 0 j = 0 u_{0j}(j=0,1,...,n):4u_{0j}=0 u0j?(j=0,1,...,n):4u0j?=0
（2）下邊界點 u i 0 ( i = 0 , 1 , . . . , n ) : 4 u i 0 = 0 u_{i0}(i=0,1,...,n):4u_{i0}=0 ui0?(i=0,1,...,n):4ui0?=0
（3）右邊界點 u n j ( j = 0 , 1 , . . . , n ) : 4 u n j = 4 ( j h ) 2 u_{nj}(j=0,1,...,n):4u_{nj}=4(jh)^2 unj?(j=0,1,...,n):4unj?=4(jh)2
（4）上邊界點 u i n ( i = 0 , 1 , . . . , n ) : 4 u i n = 4 ( i h ) 2 u_{in}(i=0,1,...,n):4u_{in}=4(ih)^2 uin?(i=0,1,...,n):4uin?=4(ih)2
（5）左邊界點的相鄰內點 u 1 j ( j = 0 , 1 , . . . , n ) : u_{1j}(j=0,1,...,n): u1j?(j=0,1,...,n):
? ( u 2 j ? 2 u 1 j ) ? ( u 1 , j + 1 ? 2 u 1 j + u 1 , j ? 1 ) = h 2 f 1 j -(u_{2j}-2u_{1j})-(u_{1,j+1}-2u_{1j}+u_{1,j-1})=h^2f_{1j} ?(u2j??2u1j?)?(u1,j+1??2u1j?+u1,j?1?)=h2f1j?
（6）下邊界點的相鄰內點 u i 1 ( i = 0 , 1 , . . . , n ) : u_{i1}(i=0,1,...,n): ui1?(i=0,1,...,n):
? ( u i + 1 , 1 ? 2 u i 1 + u i ? 1 , 1 ) ? ( u i 2 ? 2 u i 1 ) = h 2 f i 1 -(u_{i+1,1}-2u_{i1}+u_{i-1,1})-(u_{i2}-2u_{i1})=h^2f_{i1} ?(ui+1,1??2ui1?+ui?1,1?)?(ui2??2ui1?)=h2fi1?
（7）右邊界點的相鄰內點 u n ? 1 , j ( j = 0 , 1 , . . . , n ) : u_{n-1,j}(j=0,1,...,n): un?1,j?(j=0,1,...,n):
? ( ? 2 u n ? 1 , j + u n ? 2 , j ) ? ( u n ? 1 , j + 1 ? 2 u n ? 1 , j + u n ? 1 , j ? 1 ) = h 2 f n ? 1 , j + ( j h ) 2 -(-2u_{n-1,j}+u_{n-2,j})-(u_{n-1,j+1}-2u_{n-1,j}+u_{n-1,j-1})=h^2f_{n-1,j}+(jh)^2 ?(?2un?1,j?+un?2,j?)?(un?1,j+1??2un?1,j?+un?1,j?1?)=h2fn?1,j?+(jh)2
（8）上邊界點的相鄰內點 u i , n ? 1 ( i = 0 , 1 , . . . , n ) : u_{i,n-1}(i=0,1,...,n): ui,n?1?(i=0,1,...,n):
? ( u i + 1 , n ? 1 ? 2 u i , n ? 1 + u i ? 1 , n ? 1 ) ? ( ? 2 u i , n ? 1 + u i , n ? 2 ) = h 2 f i , n ? 1 + ( i h ) 2 -(u_{i+1,n-1}-2u_{i,n-1}+u_{i-1,n-1})-(-2u_{i,n-1}+u_{i,n-2})=h^2f_{i,n-1}+(ih)^2 ?(ui+1,n?1??2ui,n?1?+ui?1,n?1?)?(?2ui,n?1?+ui,n?2?)=h2fi,n?1?+(ih)2
對于同時滿足上述多個情況的節點，只需根據實際對上述相應格式進行一定的疊加即可.

3.根據五點差分格式和邊界條件（步驟2）對 K , U 和 F K,U和F K,U和F賦值（其中選取邊界值的算術平均值 Q = 0.3242 Q=0.3242 Q=0.3242為初始近似值賦于未知向量 U U U.）

4.取最大迭代次數N=500,迭代誤差限 e p ep ep分別取 1 e ? 2 1e-2 1e?2和 1 e ? 4 1e-4 1e?4，采用 G a u s s ? s e i d e l Gauss-seidel Gauss?seidel迭代求解方程組 K U = F KU=F KU=F,將求解結果放到 ( n + 1 ) 2 (n+1)^2 (n+1)2維向量 u u u中.

5.將 u u u中的元素按一定次序存盤到 ( n + 1 ) (n+1) (n+1)階矩陣 U 1 U_1 U1?中，同時通過理論解 u ( x , y ) = x 2 y 2 u(x,y)=x^2y^2 u(x,y)=x2y2產生相應節點的函式值并存盤到 ( n + 1 ) (n+1) (n+1)階矩陣 U 2 U_2 U2?中,將 U 1 , U 2 U_1,U_2 U1?,U2?的元素進行比較和分析.

6.求解結果

表2：理論解得到的節點函式值（ U 2 中 元 素 U_2中元素 U2?中元素）

表3： e p = 1 e ? 2 ep=1e-2 ep=1e?2時得到的節點函式值（ U 1 中 元 素 U_1中元素 U1?中元素，此時迭代次數 k = 8 k=8 k=8）

表4： e p = 1 e ? 4 ep=1e-4 ep=1e?4時得到的節點函式值（ U 1 中 元 素 U_1中元素 U1?中元素，此時迭代次數 k = 29 k=29 k=29）

通過mesh函式分別將表2-表4的元素可視化，依次得到圖1，圖2，圖3如下.

圖1：

圖2：

圖3：

7.總結

通過結果可看出，當迭代誤差限不斷減小到時，迭代次數從8次上升到29次，方程組數值解也更接近于理論解.本次實驗選取大小適中，能夠更普遍地反映求解區域中各節點的數值解取值，同時也方便直觀地進行比較.可以推斷當誤差限趨于0時，迭代次數將趨于無窮，此時五點差分格式計算得到的數值解將無限趨近于理論解.但由于選取的誤差限較小，節點個數選取仍然較小，我們不易直觀地從圖1至圖3中直觀地觀察到兩次數值結果同理論解之間的差異.

8.Matlab代碼

1.主函式

%% 所求方程為KU=F...
...將二維問題當成一維問題求解，即二維的（n+1）^2個節點按順序分行放到（n+1）^2維向量U中

clc
n=7;       %網格數
h=1/n;     %步長
F=zeros((n+1)^2,1); %KU=F
K=zeros(size(F));   %K為（n+1）^2階方陣
U=ones((n+1)^2,1);  %K為（n+1）^2維向量
U0=zeros((n+1)^2,1);%U0用作存邊值條件，為初始向量做準備
for i=1:n+1    %對U0中邊界點的位置賦值
    U0(i)=0;            %下邊界點
    U0(n^2+n+i)=(i*h)^2;%右邊界點
    U0(1+(i-1)*(n+1))=0;%左邊界點
    U0(i*(n+1))=(i*h)^2;%上邊界點
end
%% 下面對K賦值
for i=1:(n+1)^2
    K(i,i)=4;
end
for i=1:n-1 %五點中的右點
    for j=2+i*(n+1):(n-1)+i*(n+1)
        K(j,j+1)=-1;
    end
end
for i=1:n-1 %五點中的左點
    for j=3+i*(n+1):n+i*(n+1)
        K(j,j-1)=-1;
    end
end
for i=1:n-2 %五點中的上點
    for j=2+i*(n+1):n+i*(n+1)
        K(j,j+n+1)=-1;
    end
end
for i=2:n-1 %五點中的下點
    for j=2+i*(n+1):n+i*(n+1)
        K(j,j-n-1)=-1;
    end
end
%% 下面對F賦值
%下面賦值非特殊的點
for i=1:n-1
    for j=2+i*(n+1):n+i*(n+1)
        F(j)=-2*((floor(j/(n+1)))^2+(mod(j,n+1)-1)^2)*h^4;
    end
end
%下面賦值與邊界點相鄰的內點(左邊和下邊為齊次，不用管)
for i=1:n-1 %右邊的點
    F(n+i*(n+1))=F(n+i*(n+1))+(i*h)^2;  
end
for i=1:n-1%上邊的點
    F((n+1)*(n-1)+i+1)=F((n+1)*(n-1)+i+1)+(i*h)^2;
end
%下面賦值F中邊界點
for i=i:n+1
    F(i)=0;%下邊界點
end
for i=1:n-1
    F(1+i*(n+1))=0;%左邊界點
    F((i+1)*(n+1))=4*(i*h)^2;%右邊界點
end
for i=(n+1)*n+1:(n+1)^2 %上邊界點
    if i==(n+1)^2
        F(i)=4*(n^2)*(h^2);
    else
        F(i)=4*((mod(i,n+1)-1)^2)*(h^2);
    end
end

%% 下面對U賦初值
d=sum(sum(U0))/4/n;  %使用邊界點值的算術平均值作為U的初值以加快迭代次數
for i=1:(n+1)^2
    U(i)=d;
end
%% 使用高斯賽德迭代求解KU=F
ep=1e-4;  %誤差限，可根據需要更改，本次我們使用1e-2和1e-4
N=500;    %最大迭代次數
u=Gauss(K,F,U,ep,N);  %將求解結果存到向量u中
%% 下面將解u放到矩陣U1中，并與真解U2作比較
U1=zeros(n+1);%U1用于存放向量u中元素二維化后的資料
U2=zeros(n+1);%U2用于存放真解
for i=1:n+1%u的一維資料放到二維矩陣U1中
    for j=1+(i-1)*(n+1):i*(n+1)
        U1(i,j-(i-1)*(n+1))=u(j);
    end
end
for i=1:n+1%真解產生的節點函式值放到U2中
    for j=1:n+1
        U2(i,j)=(((i-1)*h)^2)*(((j-1)*h)^2);
    end
end
%% 輸出解
U1
U2
%% 畫圖
X=linspace(0,1,n+1);
Y=linspace(0,1,n+1);
mesh(X,Y,U2);%當要畫數值解的圖時，U2改為U1即可
hold on;
title('理論解影像');%當要畫數值解的圖時，'理論解影像'改為'數值解影像'即可

2.Gauss-seidel迭代

function x=Gauss(A,b,x0,ep,N)

%用于Gauss-seidel迭代法解線性方程組Ax=b
%A，b,x0分別為系數矩陣，右端向量和初始向量（初始向量默認為零向量）
%ep為精度（1e-3），N為最大迭代次數（默認500次），x回傳數值解向量

n=length(b);
if nargin<5
    N=5000;
end
if nargin<4
    ep=1e-6;
end
if nargin<3
    x0=zeros(n,1);
end
x=zeros(n,1);
k=0;
while k<N
    for i=1:n
        if i==1
            x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
        elseif i==n
                x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);
        else
                x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);
        end
    end
    if norm(x-x0,inf)<ep
        break;
    end
    x0=x;
    k=k+1;
end
if k==N
    Warning('已到達迭代次數上限');
end
disp(['k=',num2str(k)])
end