Laplace九點(diǎn)差分格式

1、laplace九點(diǎn)差分格式 中南林業(yè)科技高校 本科課程設(shè)計(jì)說(shuō)明書 學(xué) 院： 理學(xué)院 專業(yè)班級(jí)： 2021級(jí)信息與計(jì)算科學(xué)二班 課 程： 科學(xué)計(jì)算課程設(shè)計(jì) 論文題目： laplace方程九點(diǎn)差分格式 指導(dǎo)老師： 陳宏斌 2021年6月 二維橢圓邊值問題的九點(diǎn)差分格式 1問題：laplace方程 uxx uyy 0, x,y g, g是xy平面上一有界區(qū)域，其邊界 為分段光滑曲線. 在 上u滿意下列邊值條件： u (x,y)(drichlet邊值條件). 在此考慮g為正方形區(qū)域，g=(x,y) | axb, ayb. 背景：拉普拉斯方程（laplace'sequation），又名調(diào)和方程、

2、位勢(shì)方程，是一種偏微分方程。由于由法國(guó)數(shù)學(xué)家拉普拉斯首先提出而得名。求解拉普拉斯方程是電磁學(xué)、天文學(xué)和流體力學(xué)等領(lǐng)域常常遇到的一類重要的數(shù)學(xué)問題，由于這種方程以勢(shì)函數(shù)的形式描寫了電場(chǎng)、引力場(chǎng)和流場(chǎng)等物理對(duì)象（一般統(tǒng)稱為“保守場(chǎng)”或“有勢(shì)場(chǎng)”）的性質(zhì)。 2區(qū)域剖分 區(qū)域g是一個(gè)正方形區(qū)域，邊界axb,ayb. 分別沿x，y軸，在a,b上取n+1個(gè)節(jié)點(diǎn)：a x0 x1 xn b,a y0 y1 yn b. 則步長(zhǎng) h=(b-a)/n; 3九點(diǎn)差分格式 3.1向量形式 利用taylor展式，有 （3.1） u(xi 1,yj) 2u(xi,yj) u(xi 1,yj) 2 h x2 68 h4 u(

3、xi,yj)2h6 u(xi,yj)8 o(h),68360 x8! x （3.2） 2u(xi,yj) 4 h2 u(xi,yj) 12 x4 4 h2 u(xi,yj) 12 y4 u(xi,yj 1) 2u(xi,yj) u(xi,yj 1) h 4 6 2 6 8 2u(xi,yj) y2 h u(xi,yj)2h u(xi,yj) o(h8),68360 y8! y 將（3.1），（3.2）兩式相加，則得 4u(xi,yj) 4u(xi,yj) h4 6u(xi,yj) 6u(xi,yj) h2 hu(xi,yj) u(xi,yj) 4466 12 x y x y 360 8u(xi

4、,yj) 8u(xi,yj) 2h6 o(h8), 88 8! x y 其中 hu(xi,yj) u(xi 1,yj) 2u xi,yj u(xi 1,yj) h 2 u(xi,yj 1) 2u(xi,yj) u(xi,yj 1) h 2 , u(xi,yj) 又 2u(xi,yj) x2 2u(xi,yj) y2 0. 6u(xi,yj) x6 6u(xi,yj) y6 2u xi,yj 2u(xi,yj) 6u(xi,yj) 6u(xi,yj) 4 4 224224 x4 y4 x y x y x y 2u xi,yj 2u xi,yj 4 0, 22 22 x y x y 6u(xi,y

5、j) x y 4 2 6u(xi,yj) x2 y4 4u xi,yj 4u xi,yj 2 2u xi,yj 2u xi,yj 4u xi,yj 4u xi,yj 2 2 2 2,442222222 x y x y x y x y x y 而 4u xi,yj x2 y2 u''xx(xi,yj 1) 2u''xx(xi,yj) u''xx(xi,yj 1) h2 2u xi,yj h4 8u xi,yj h2 4 o(h6)4 262 360 y x12 y x 1 u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2

6、(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4 4 4u(xi,yj) 6u(xi,yj 1) h2 6u(xi,yj 1) 6u(xi,yj) 6u(xi,yj 1) h2 6u(xi,yj)1 u(xi,yj 1) 2 2 4446662412 x x x360 x x x 12 x y 8 h4 u(xi,yj) o(h6)26360 x y 1 u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj) u

7、(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4 66888 h2 u(xi,yj)h2 u(xi,yj)1h4 u(xi,yj)h4 u(xi,yj)h4 u(xi,yj) * o(h6)422444622612 x y12 x y1212 x y360 x y360 x y 1 u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4 888 1h4 u(xi,yj)h4 u(xi,yj)h4 u

8、(xi,yj) * o(h6)4462261212 x y360 x y360 x y 1 u(xi 1,yj 1) 2u xi,yj 1 u(xi 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) h4 8 h4 u(xi,yj) o(h6)26720 x y 因此 4 8u(xi,yj) 8u(xi,yj) h2 u xi,yj 2h6 8 o(h) hu(xi,yj) u(xi,yj) 8822 8! x y6 x y 1 u(xi 1,yj 1) 2u xi,yj 1 u(xi

9、 1,yj 1) 2(u(xi 1,yj) 2u(xi,yj) u(xi 1,yj) u(xi 1,yj 1) 2u(xi,yj 1) u(xi 1,yj 1) 6h2 8 8u(xi,yj) 8u(xi,yj) h2h4 u(xi,yj)2h6 8 o(h) u(xi,yj) 8826 8! x y6720 x y 舍去截?cái)嗾`差項(xiàng)，可得laplace 方程的九點(diǎn)差分格式 1 20ui,j 4ui 1,j 4ui 1,j 4ui,j 1 4ui,j 1 ui 1,j 1 ui 1,j 1 ui 1,j 1 ui 1,j 1 0 6h2 從而可得差分算子 h的截?cái)嗾`差 8 8u(xi,yj) 8

10、u(xi,yj) h2h4 u(xi,yj)2h6 o(h8) ri,j(u) 88 8! x y6720 x2 y6 40h6 ri,j(u) m8, 3*8! 其中m8是u的8階偏導(dǎo)數(shù)的肯定值于考慮區(qū)域g 1,且u是laplace 方程的 的光滑解. 3.2 矩陣形式 au b其中 a1 a 2 a a2 a1 204 4 204 , a1 , a2 4 204 a1 4 20 (n 1)*(n 1) (n 1)*(n 1) a2 a2 a1a2 a2 4 1 1 4 1 1 41 ; 1 4 (n 1)*(n 1) u u1,u2, ,un 1 '， ui ui1,ui2,ui3

11、, ,ui,n 1 ',i 1,2,3, ,n 1; b b1,b2,b3, ,bn 1 ', b1 u00 4u01 u02 4u10 u20, u01 4u02 u03, u02 4u03 u04, , u0,n 3 4u0,n 2 u0,n 1, u0,n 2 4u0,n 1 u0,n 4u1,n u2,n ,bi ui 1,0 4ui,0 ui 1,0,0, ,0, ui 1,n 4ui,n ui 1,n ,i 2,3,4, n 2; bn 1 un 2,0 4un 1,0 un0 4un1 un2, un1 4un2 un3, un2 4un3 un4, , un,n

12、 3 4un,n 2 un,n 1, un,n 2 4un,n 1 un 2,n 4un 1,n un,n , 4數(shù)值試驗(yàn) 2u 2u x2 y2 0; u(x,y) ex(siny cosy); 0 x 1,0 y 1. 邊值條件： u(0,y) siny cosy, u(1,y) e(siny cosy), xu(x,0) e, x u(x,1) e(sin1 cos1). 5編程實(shí)現(xiàn) 5.1算法 方程au b中，a是大型稀疏矩陣矩陣，用塊gauss-seidel迭代法解； 步驟：1輸入?yún)^(qū)間上限、下限，節(jié)點(diǎn)； 2分解矩陣a，a d l u;則迭代矩陣：b (d l) 1u； 3由分塊矩陣乘

13、法得到塊gauss-seidel 迭代法的詳細(xì)形式： u(0) 0; (k 1)(k) a2u2 b1; a1u1 (k 1)k 1)k) a2ui( a2ui( 11 bi,i 2,3, ,n 2; a1ui (k 1)(k 1) au au1n 12n 2 bn 1; k 0,1,2, . 4輸出. 5.2流程圖 5.4編寫代碼 5.5測(cè)試 例1 步長(zhǎng) h=1/6 數(shù)值解： 1.3362 1.4653 1.5609 1.6172 1.6393 1.5423 1.6926 1.8067 1.8649 1.8897 1.7216 1.9659 2.1132 2.1590 2.0944 1.98

14、89 2.3380 2.5190 2.5661 2.4093 2.5776 2.8509 3.0488 3.1403 3.1553 精確解： 1.3610 1.5029 1.6031 1.6589 1.6688 1.6078 1.7754 1.8939 1.9598 1.9714 1.8994 2.0974 2.2373 2.3152 2.3290 2.2439 2.4778 2.6431 2.7351 2.7513 2.6508 2.9272 3.1224 3.2312 3.2503 例2 步長(zhǎng)h=1/8 數(shù)值解： 1.2377 1.3340 1.4182 1.4831 1.5260 1.5

15、498 1.3635 1.4637 1.5598 1.6322 1.6762 1.6973 1.4509 1.6028 1.7269 1.8113 1.8517 1.8508 1.5485 1.7708 1.9349 2.0355 2.0702 2.0351 1.6894 1.9993 2.2064 2.3254 2.3588 2.2917 1.9420 2.3330 2.5632 2.6982 2.7441 2.6806 2.5580 2.7946 3.0092 3.1568 3.2312 3.2330 精確解： 1.2656 1.3783 1.4694 1.5377 1.5819 1.60

16、15 1.4341 1.5618 1.6651 1.7424 1.7926 1.8147 1.6250 1.7697 1.8868 1.9744 2.0313 2.0564 1.8414 2.0054 2.1380 2.2373 2.3017 2.3302 2.0866 2.2724 2.4227 2.5352 2.6082 2.6404 2.3644 2.5749 2.7453 2.8728 2.9555 2.9920 2.6792 2.9178 3.1108 3.2553 3.3490 3.3904 取步長(zhǎng)h1=h2=1/10，其數(shù)值解和精確解的曲面如下： 1.5621 1.7197 1.

17、8206 1.9311 2.0957 2.4086 3.2172 1.5961 1.8086 2.0494 2.3223 2.6315 2.9819 3.3789 5.6附源程序代碼 %laplace方程九點(diǎn)差分格式 clear; %a=input('輸入?yún)^(qū)間下限：'); %c=input('輸入?yún)^(qū)間上限：'); n=input('輸入節(jié)點(diǎn)數(shù)：'); a=0;c=1;%n=6; n=10;%迭代次數(shù) %網(wǎng)格剖分 h=(c-a)/n; x=a+h:h:c; y=x; %解au=b %對(duì)矩陣b賦值 for p=1:n-1 for q=1:n-1 b(

18、p,q)=0; if p=1 if q=1 b(p,q)=-4*(sin(y(q)+cos(y(q)-1-(sin(y(2)+cos(y(2)-4*exp(x(1)-exp(x(2);%邊值條件 end if q=n-1 b(p,q)=-lap9_u(a,y(n-2)-4*lap9_u(a,y(n-1)-lap9_u(a,c)-4*lap9_u(x(1),c)-lap9_u(x(2),c); end if q1qn-1 b(p,q)=-lap9_u(a,y(q-1)-4*lap9_u(a,y(q)-lap9_u(a,y(q+1); end end if p1pn-1 if q=1 b(p,q)

19、=-lap9_u(x(1),a)-4*lap9_u(x(2),a)-lap9_u(x(3),a); end if q=n-1 b(p,q)=-lap9_u(x(1),c)-4*lap9_u(x(2),c)-lap9_u(x(3),c); end end if p=n-1 if q=1 b(p,q)=-lap9_u(x(n-2),a)-4*lap9_u(x(n-1),a)-lap9_u(c,a)-4*lap9_u(c,y(1)-lap9_u(c,y(2); end if q=n-1 b(p,q)=-lap9_u(x(n-2),c)-4*lap9_u(x(n-1),c)-lap9_u(c,c)-4*lap9_u(c,y(n-1)-lap9_u(c,y(n-2); end if q1qn-1 b(p,q)=-lap9_u(c,y(q-1)-4*lap9_u(c,y(q)-lap9_u(c,y(q+1); end end end end %對(duì)矩陣u賦初值 for p=1:n-1 for q=1:n-1 u0(p,q)=1; end end %賦值給