Possion方程的差分法求解(共10頁)

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

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

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

4、my)/2 個點的函數(shù)值，但對程序效率的影響并不是很大，可以接受。以下具體說明只需在G1上求解poisson方程的原因所求方程為：設直線l是經(jīng)過原點o的任意一條直線，其方程為：y=kx設p(x,y)是區(qū)域內(nèi)任一點，則其關(guān)于l對稱的點為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)點p的橫坐標x=rcos() y=rsin()則p關(guān)于直線l的對稱點為p(rcos(2-),rsin(2-)由上述證明可知u(p)=u(p)。由和

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

6、,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)點(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和右端項f，顯然A是對稱正定矩陣，也是稀疏矩陣因此可以采用逐次超松弛法，共軛梯度法和交替方向法萊求解，但為

7、了方便程序設計采用了matlab的運算符來求解u。針對本題的Poisson方程，對S1+G1形成的矩形網(wǎng)格的五點差分的具體分析。對S1+G1形成的矩形五點差分的具體分析。ny+1ny3211 2 3 4 i nx-1 nx nx+1G1S2S11/2O3/21/2紅色線條圍成的區(qū)域為G1，L為紅色斜邊，S1為L左側(cè)的三角形，S2為L右側(cè)的三角形。由對稱性知，S1和S2中關(guān)于L相互對稱的網(wǎng)格點其上U的函數(shù)值是相同的。按從左到右，從下到上的次序排列網(wǎng)點(ij)。（1）對于三角形S1中的網(wǎng)點有u(i,j)， nyij1，有u(i,j)-u(j,i)=0 所以S1中網(wǎng)點的差分就為：u(i,j)-u(j

8、,i)=0（2）由于原點o點的特殊性，其鄰點中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)格上最下面一行上除了原點o的所有正則網(wǎng)格內(nèi)點，由對稱性得u(1,i)的鄰點中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)點u(1,nx)其差分為：4u(i,j)-u(i-1,j)-2*u(i,j+1

9、)=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)點，其差分為：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)對于網(wǎng)格最后一列的所有非正則內(nèi)點u(:,nx),其差分為：4u(nx,j

10、)-u(nx,j-1)-u(nx,j+1)-u(nx-1,j)=h2*f(i,j)(6)在求出了所有的內(nèi)點后，對于剩下的邊界點賦值有：U(ny+1,ny+1:nx)=0%最上面一行上的邊界點U(1:ny+1,nx+1)=0%最右側(cè)一列的邊界點u(ny+1,1:ny)=u(1:ny,ny+1);%三角形S1和S2邊界上的網(wǎng)點，它們是S1的邊界點，但是整個求解區(qū)域的內(nèi)點。三、程序設計及分析function poisson5spot(h)if narginj%因此令A(i,j)=1 A(j,i)=-1%所以在本程序中多計算了(ny2-my)/2 個點的函數(shù)值%但對程序的影響并不是很大for i=2:

11、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)格點上的函數(shù)值%x=gmres(A,b,100);%按順序?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)格最上面一行的點賦值u(ny+1,1:ny)=u(1:ny,ny+1); %=作Poisson方

12、程在區(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 四、實驗結(jié)果1在區(qū)域G1上求解后的圖形顯示:圖形表示了在S1+G1區(qū)域上的函數(shù)圖像，而不是單純的G1上的函數(shù)圖像。2通過拼接后圖形顯示如下：由上圖可知，除了邊界點外網(wǎng)格點上的函數(shù)值都有u(i,j)0.Lhuij0，對任意(xi,yj)Gh，uij不可能在內(nèi)點取得負的極小，與極值定理相符合。3、翻轉(zhuǎn)后圖形如下：4 網(wǎng)格間距h=0.01時得到的圖形： 五、實驗分析本次實驗將求解區(qū)域先利用對稱性將其縮小為原區(qū)域的1/