Poisson方程九點(diǎn)差分格式米瑞琪

1、數(shù)值實(shí)驗(yàn)報(bào)告I實(shí)驗(yàn)名稱Poisson方程九點(diǎn)差分格式實(shí)驗(yàn)時(shí)間2016年4月15日姓名米瑞琪班級(jí)信息1303學(xué)號(hào)04成績(jī)、實(shí)驗(yàn)?zāi)康?，?nèi)容1、理解Poisson方程九點(diǎn)差分格式的構(gòu)造原理；2、理解因?yàn)榫W(wǎng)格點(diǎn)的不同排序方式造成的系數(shù)矩陣格式的差異;3、學(xué)會(huì)利用matlab的spdiags(),kron()函數(shù)生成系數(shù)矩陣;二、算法描述針對(duì)一個(gè)Poisson方程問(wèn)題：Au 二-f(X, y)在Poisson方程五點(diǎn)差分格式的基礎(chǔ)上，采用 Taylor展開分析五點(diǎn)差分算子的截?cái)嗾`差，可以 得到：AhU(肛，yj)1( 7/(“力)/u&.y)1二 Au(x“Vj) + - hi + + 0(h

2、); AuGi,力)+ 14 dxdy / J n Waufxi, yj) a'uxirv) h/ + 丸(wj hi2 ： h22 -J + 0(力取 的叭&a/ /12a/羽(1)為了提高算子截?cái)嗾`差的精度，在(1)式中配湊出了差分算子的形式，將原Poisson方程代入(1) 式有:、/、 1 / 47 鏟儼u(Xi，Vj)+ h/血Vj)"¥)=而"力卜視 .+加匐+卜獲落十0")，有:(/u(xi yj)1 i ”=-UxK(Ki.Vj + 1)- 2uxx(Xi- ¥j)+ uxx(xir Yj - i) + 0(h?

3、2)=-h7h?u(xj + i, yj + i) - 2u(xlr /j + i) + u(xj - n yj + JhiVu(x! h yj 4 i)hj12 d乂4u(xj + i. yj) - 2u(xj. yj) + u(xi i, yj)2h,丸(x#),5 天u(xj + Ip Vj J _ 2u(x(F Vj - 1) + U(Xi - 1, Vj - 1) hm(Xi+ yj - J+hi212 dJLhU(Xj( yj i) - N/Xxj , ¥j) + LhuUih yj _ J -1212h/h/Mx - Vj * i)Lhu(xl±yj i J

4、- SUutxi, yJ)+ Lhu(xn Yj _ 1)- c，血（xi將（3）代回（2）可得h, + h/AhU(Xi. v。+ -y(4u(xir Vj) - 2(u(xj, yj + 1) + u(xj,力i) + uGi + i ¥j) - u(xj - b yj)12hiT h?+ ufxj + uYj i) + u(xj 4 b Vj - i) + u(xj - n Vj i) + 口(肛-i"-i)/千Vj)得到Poisson方程的九點(diǎn)差分格式:I2h/h/(4叫-2(Uj j + 1 + u(, j - 1 + u( + i5 j + Uj 1, j) +

5、 5疊+曲a? Jy在計(jì)算機(jī)上實(shí)現(xiàn)（4）式，需要在五點(diǎn)差分格式Al】u 1=_fi j的基礎(chǔ)上在等式兩端分別增加一部分，將等式左側(cè)新增的部分寫成緊湊格式，有:l-2U11U12tll.Ni - 1U22UN, - 1 1U* 1 2對(duì)于該矩陣，可以看成是兩個(gè)矩陣的組合:4- 2r2 1 v-24C =以及-211- 2 1X-211- 2 1 *11- 2-211- 2 1 » «* *-21-211- 2 1""a"«i*1- 2則生成這兩個(gè)矩陣可以采用 Kroncker生成，方法類似于五點(diǎn)差分格式c對(duì)于右端添加的關(guān)于f(x,y)的

6、二階導(dǎo)數(shù)，可以采用中心差分格式進(jìn)行近似代替，即:寫成相應(yīng)的緊湊格式有:11-4 11«I-«4 .-4I 111121- 4-4111-41備+ '1111該式中的矩陣又可以分解為兩個(gè)矩陣的和:11"111119«* 11111111三.程序代碼根據(jù)上述分析，可以寫出程序代碼為：clcclear網(wǎng)格剖分xa=0;xb=1;N1=64;h1=(xb-xa)/N1;x=xa+1:(N1-1)*h1;ya=0;yb=1;N2=64;h2=(yb-ya)/N2;y=ya+1:(N2-1)*h2;9往成網(wǎng)格點(diǎn)與邊界點(diǎn)處f的函數(shù)值，其中f_b1是中間全部為0

7、的向量，f_b2是每一個(gè)分塊中間為 0的向量f=(x,y)-(2*piA2)*exp(pi*(x+y)*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y); f_b1=zeros(N1-1)*(N2-1),1);f_b2=zeros(N1-1)*(N2-1),1);for i=1:N1-1f_b2(i-1)*(N2-1)+1)=f(x(i),ya);f_b2(i*(N2-1)=f(x(i),yb);for j=1:N2-1f_in(i-1)*(N2-1)+j)=f(x(i),y(j); if i=1f_b1(j)=f(xa,y(j); endif i=N1-1f_

8、b1(i-1)*(N2-1)+j)=f(xb,y(j); endend end f_in=f_in'熱成迭41矩陣Ae2=ones(N1-1,1);K1=spdiags(e2,-2*e2,e2,-1,0,1,N1-1,N1-1);e3=ones(N2-1,1);I1=spdiags(e3,0,N2-1,N2-1);A=kron(K1,I1);A=A/h1A2;9往成迭41矩陣B和CI2=spdiags(e2,0,N1-1,N1-1);K2=spdiags(e3,-2*e3,e3,-1,0,1,N2-1,N2-1);B=kron(I2,K2);C=-2*(h1A2+h2A2)/(12*h

9、1A2*h2A2)*B;B=B/h2A2;9往成迭41矩陣DS1=spdiags(e2 e2,-1 1,N1-1,N1-1);S2=spdiags(e3,0,N2-1,N2-1);D=(h1A2+h2A2)/(12*h1A2*h2A2)*kron(S1,K2);9往成f的中心差分格式對(duì)應(yīng)矩陣 E和FE=kron(S1,S2);K3=spdiags(e3,-4*e3,e3,-1 0 1,N2-1,N2-1);F=kron(I2,K3);f_vec=-f_in-(1/12)*(E+F)*f_in+f_b1+f_b2);熱成網(wǎng)麻點(diǎn)X,Y=meshgrid(x,y);%+算網(wǎng)格函數(shù)值u_d=(A+B+

10、C+D)f_vec;u0=zeros(length(f vec),1);%+算誤差u_real=(x,y)exp(pi*(x+y)*sin(pi*x).*sin(pi*y);for i=1:N1-1u_m(i-1)*(N2-1)+1:i*(N2-1)=u_real(x(i),y);endu_v=u_m'err_d=max(abs(u_d-u_v);sol=reshape(u_d,N2-1,N1-1);mesh(X,Y,sol)四.數(shù)值結(jié)果,將計(jì)算出的誤差列表如下:針對(duì)課本P93給出的問(wèn)題，分別采用步長(zhǎng) :高J k五點(diǎn)差分格式誤差九點(diǎn)差分格式誤差1641128可見采用九點(diǎn)差分格式可以進(jìn)一步縮小誤差，達(dá)到更高階的精度。五.計(jì)算中出現(xiàn)的問(wèn)題，解