Possion方程的差分法求解

1、Possion 方程的差分方法課程名稱： 微分方程數(shù)值解 學(xué)生姓名： 張弘 一、問題描述二、問題分析I偏微分方程的數(shù)值解法主要有有限差分法和Galerkin有限元法。用差分法和有限元法將連續(xù)問題離散化的步驟是：1、對求解區(qū)域做網(wǎng)格剖分用有限個網(wǎng)格節(jié)點(diǎn)代替連續(xù)區(qū)域。2、微分算子離散化。3、把微分方程的定解問題化為線性代數(shù)方程組的求解問題。差分法和有限元方法的主要區(qū)別是離散化的第二步，差分法從定解問題的微分或積分形式出發(fā)，用數(shù)值微商或數(shù)值積分公式到處相應(yīng)的線性代數(shù)方程組，有限元法從定解問題的變分形式出發(fā)，用Ritz-Galerkin法導(dǎo)出相應(yīng)的線性代數(shù)方程組。差分法的基本問題有：（1） 對求解區(qū)域

2、做網(wǎng)格剖分一維情形為把區(qū)間分為等距或不等距的區(qū)間，二維情形是把區(qū)域分割成均勻或不均勻的矩形，其邊與坐標(biāo)軸平行，也可分割成小三角形或凸四邊形。（2） 構(gòu)造逼近微分方程定解問題的差分格式有兩種構(gòu)造差分格式的方法：直接差分化法和有限體積法。（3） 差分解的存在唯一性，收斂性和穩(wěn)定性研究（4） 差分方程的解法：逐次超松弛法，交替方向法，共軛梯度法。II（1）由題目可知，本題需要考慮矩形網(wǎng)的五點(diǎn)差分格式。題目給出的為第一邊值條件，將十字形圖形的中心放于坐標(biāo)原點(diǎn)處，如下圖所示：OS1G1S2由圖形可知，所給區(qū)域可以看出是由8個梯形G1通過旋轉(zhuǎn)、翻轉(zhuǎn)拼接而成。因此為了方便計算、減少計算量，只針對G1進(jìn)行網(wǎng)格

3、剖分，用5點(diǎn)差分格式進(jìn)行求解。但是由于G1是直角梯形，進(jìn)行網(wǎng)格剖分時會出現(xiàn)x軸方向網(wǎng)格點(diǎn)個數(shù)不同的現(xiàn)象，不利于有差分形成的系數(shù)矩陣的生成，所以將三角形S1和梯形G1合在一起形成一個矩形，其區(qū)域?yàn)?,3/2×0,1/2。如果采用等距差分，并且有hx=hy=h，設(shè)步長為h，x=0:h:3/2;y=0:h:1/2;nx=length(x)-1;%為所求區(qū)域中x軸上內(nèi)點(diǎn)的個數(shù)ny=length(y)-1; %為所求區(qū)域中y軸上內(nèi)點(diǎn)的個數(shù)對于原來的梯形G1來說，有網(wǎng)格內(nèi)點(diǎn)nx*ny-(ny2-my)/2對于矩形區(qū)域S1+G1而言，網(wǎng)格內(nèi)點(diǎn)數(shù)為nx*ny，所以在后面的程序中要比在梯形G1中多計算

4、了(ny2-my)/2 個點(diǎn)的函數(shù)值，但對程序效率的影響并不是很大，可以接受。以下具體說明只需在G1上求解poisson方程的原因所求方程為：設(shè)直線l是經(jīng)過原點(diǎn)o的任意一條直線，其方程為：y=kx設(shè)p(x,y)是區(qū)域內(nèi)任一點(diǎn)，則其關(guān)于l對稱的點(diǎn)為p'(s,t)S=(k-1)2+|2y)/(k2+1) t=(2kx+(k2-1)y)/(k2+1)則同理可得：將u(s,t)代替u(x,y)得：Uxx+uyy=uss+utt=1其依然滿足上述方程。令=arctan(k)點(diǎn)p的橫坐標(biāo)x=rcos() y=rsin()則p關(guān)于直線l的對稱點(diǎn)為p'(rcos(2-),rsin(2-)由上述

5、證明可知u(p)=u(p')。由和p點(diǎn)的任意性知，對于函數(shù)u圖像上的任意一點(diǎn)p，其關(guān)于任意一條經(jīng)過原點(diǎn)直線l的對稱點(diǎn)p'都在u的圖像上，即u(+)=u(),即u關(guān)于原點(diǎn)是旋轉(zhuǎn)對稱的。當(dāng)l為x軸時，有u(x,y)=u(x,-y)l為y軸時，u(x,y)=u(-x,y)坐標(biāo)軸旋轉(zhuǎn)不改變方程的形式，所以函數(shù)在直角坐標(biāo)系中u關(guān)于原點(diǎn)是旋轉(zhuǎn)對稱的，又所求十字形區(qū)域關(guān)于x，y軸是軸對稱和關(guān)于原點(diǎn)中心對稱的，因此可通過直求解區(qū)域G1，就可以知道函數(shù)在整個十字形區(qū)域的圖像。（2）對S1+G1形成的矩形進(jìn)行正方形網(wǎng)格剖分，則求解Poisson方程的差分格式和化為如下形式：對于正則內(nèi)點(diǎn)其差分如下：

6、-uij=1/h2*(-u(i,j+1)-u(i-1,j)+4u(i,j)-u(i+1,j)-u(i,j-1)=fij(1)對于矩形區(qū)域S1+G1，nx=(xb-xa)/h-1ny=(yb-ya)/h-1按從左到右，從下到上的次序排列網(wǎng)點(diǎn)(ij)：j=1,1inx；j=2,1inx;,;j=ny,1inx，定義向量Uh=u11,u21,unx-1;u1,nx-1,u2,nx-1,uny-1,nx-1T利用邊界條件可以將（1）式寫成如下形式:其中A=為nx*ny階矩陣，I為nx階單位矩陣，B為nx階單位矩陣。B=右端向量g的元素，依賴于邊值a和右端項(xiàng)f，顯然A是對稱正定矩陣，也是稀疏矩陣因此可以

7、采用逐次超松弛法，共軛梯度法和交替方向法萊求解，但為了方便程序設(shè)計采用了matlab的'運(yùn)算符來求解u。針對本題的Poisson方程，對S1+G1形成的矩形網(wǎng)格的五點(diǎn)差分的具體分析。對S1+G1形成的矩形五點(diǎn)差分的具體分析。ny+1ny3211 2 3 4 i nx-1 nx nx+1G1S2S11/2O3/21/2紅色線條圍成的區(qū)域?yàn)镚1，L為紅色斜邊，S1為L左側(cè)的三角形，S2為L右側(cè)的三角形。由對稱性知，S1和S2中關(guān)于L相互對稱的網(wǎng)格點(diǎn)其上U的函數(shù)值是相同的。按從左到右，從下到上的次序排列網(wǎng)點(diǎn)(ij)。（1）對于三角形S1中的網(wǎng)點(diǎn)有u(i,j)， nyi>j1，有u(i,

8、j)-u(j,i)=0 所以S1中網(wǎng)點(diǎn)的差分就為：u(i,j)-u(j,i)=0（2）由于原點(diǎn)o點(diǎn)的特殊性，其鄰點(diǎn)中u(1,2)=u(1,-1)=u(-1,1)=u(2,1)所以其差分為：u(1,1)-4u(1,2)= h2*f(i,j)程序中語句：A(1:nx,nx+1:2*nx)=2*I; 和A(1,2)=-2;保障了上面差分方程的系數(shù)。（3）對于網(wǎng)格上最下面一行上除了原點(diǎn)o的所有正則網(wǎng)格內(nèi)點(diǎn)，由對稱性得u(1,i)的鄰點(diǎn)中u(1,i-1)=u(1,i+1)所以其差分為：4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)=h2*f(i,j)對于右下角的非正則內(nèi)點(diǎn)u(1

9、,nx)其差分為：4u(i,j)-u(i-1,j)-2*u(i,j+1)=h2*f(i,j)程序中的相關(guān)語句為：A(1:nx,nx+1:2*nx)=2*I; A(nx+1:2*nx,nx+1:2*nx)=B;(4)對于G1中的所有正則內(nèi)點(diǎn)，其差分為：4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j-1)-u(i,j+1)=h2*f(i,j)程序中相關(guān)語句：A(i*nx+1:(i+1)*nx,i*nx+1:(i+1)*nx)=B;A(i-1)*nx+1:i*nx,i*nx+1:(i+1)*nx)=I;A(i*nx+1:(i+1)*nx,(i-1)*nx+1:i*nx)=I;(5)對

10、于網(wǎng)格最后一列的所有非正則內(nèi)點(diǎn)u(:,nx),其差分為：4u(nx,j)-u(nx,j-1)-u(nx,j+1)-u(nx-1,j)=h2*f(i,j)(6)在求出了所有的內(nèi)點(diǎn)后，對于剩下的邊界點(diǎn)賦值有：U(ny+1,ny+1:nx)=0%最上面一行上的邊界點(diǎn)U(1:ny+1,nx+1)=0%最右側(cè)一列的邊界點(diǎn)u(ny+1,1:ny)=u(1:ny,ny+1);%三角形S1和S2邊界上的網(wǎng)點(diǎn)，它們是S1的邊界點(diǎn)，但是整個求解區(qū)域的內(nèi)點(diǎn)。三、程序設(shè)計及分析function poisson5spot(h)if nargin<1 %默認(rèn)步長h=0.02 h=0.02;endx=0:h:3/2;

11、y=0:h:1/2;nx=length(x)-1;%取區(qū)域的內(nèi)點(diǎn)ny=length(y)-1;%=構(gòu)造矩陣B=B=eye(nx,nx)*4;for i=1:nx-1 B(i,i+1)=-1;endfor i=2:nx B(i,i-1)=-1;endI=-eye(nx,nx);%=構(gòu)造系數(shù)矩陣A=A=zeros(nx*ny,nx*ny);A(1:nx,1:nx)=B;%由區(qū)域的對稱性，正方形網(wǎng)格最下面一行的差分形式%改為4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)%因?yàn)閡(i,j+1)=u(i,j-1)A(1:nx,nx+1:2*nx)=2*I; A(nx+1:2*n

12、x,nx+1:2*nx)=B;%矩形網(wǎng)格左下角第一個點(diǎn)的差分形式改為：%4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j+1)-u(i,j-1)%=4u(i,j)-2u(i,j+1)-2u(i+1,j)A(1,2)=-2;%為了方便，本程序往梯形中增加了一個三角形，以方便編程A(nx+1:2*nx,1:nx)=I;for i=2:ny-1 A(i*nx+1:(i+1)*nx,i*nx+1:(i+1)*nx)=B; A(i-1)*nx+1:i*nx,i*nx+1:(i+1)*nx)=I; A(i*nx+1:(i+1)*nx,(i-1)*nx+1:i*nx)=I;end%=%bi=h

13、*h*f;右端向量b=h*h*ones(nx*ny,1);%由于對稱性，左側(cè)三角形中有u(i,j)-u(j,i)=0 i>j%因此令A(yù)(i,j)=1 A(j,i)=-1%所以在本程序中多計算了(ny2-my)/2 個點(diǎn)的函數(shù)值%但對程序的影響并不是很大for i=2:ny A(i-1)*nx+1:(i-1)*nx+i-1,:)=0; for j=1:i-1 A(i-1)*nx+j,(i-1)*nx+j)=1; A(i-1)*nx+j,(j-1)*nx+i)=-1; b(i-1)*nx+j)=0; endendx=Ab;%為了方便采用左除求解網(wǎng)格點(diǎn)上的函數(shù)值%x=gmres(A,b,100

14、);%按順序?qū)賦值給uu=zeros(ny+1,nx+1);for i=1:ny for j=1:nx u(i,j)=x(i-1)*nx+j); endend%根據(jù)對稱性，給網(wǎng)格最上面一行的點(diǎn)賦值u(ny+1,1:ny)=u(1:ny,ny+1); %=作Poisson方程在區(qū)域上的圖形=x,y=meshgrid(0:h:3/2,0:h:1/2);hold onsurf(x,y,u)%11第一象限的第一塊區(qū)域，下面的以此類推surf(y,x,u)%12surf(-y,x,u)%21surf(-x,y,u);%22surf(-x,-y,u)%31surf(-y,-x,u)%32surf(y,-x,u)%41surf(x,-y,u);%42 四、實(shí)驗(yàn)結(jié)果1在區(qū)域G1上求解后的圖形顯示:圖形表示了在S1+G1區(qū)域上的函數(shù)圖像，而不是單純的G1上的函數(shù)圖像。2通過拼接后圖形顯示如下：由上圖可知，除了邊界點(diǎn)外網(wǎng)格點(diǎn)上的函數(shù)值都有u(i,j)>0.Lhuij>0，對任意(xi,yj)Gh，uij不可能在內(nèi)點(diǎn)取得負(fù)的極小，與極值定理相符合。3、翻轉(zhuǎn)后圖形如下：4 網(wǎng)格間距h=0.01時得到