Poisson方程九點差分格式-米瑞琪

數(shù)值實驗報告實驗名稱Poisson方程九點差分格式實驗時間實驗名稱Poisson方程九點差分格式實驗時間2016年4月15日姓名米瑞琪班級信息1303學(xué)號1309010304成績一、實驗?zāi)康模瑑?nèi)容1、理解Poisson方程九點差分格式的構(gòu)造原理；2、理解因為網(wǎng)格點的不同排序方式造成的系數(shù)矩陣格式的差異；3、學(xué)會利用matlab的spdiags(),kron()函數(shù)生成系數(shù)矩陣；二、算法描述針對一個Poisson方程問題：在Poisson方程五點差分格式的基礎(chǔ)上，采用Taylor展開分析五點差分算子的截斷誤差，可以得到：為了提高算子截斷誤差的精度，在(1)式中配湊出了差分算子的形式，將原Poisson方程代入(1)式有：考慮，有:將（3）代回（2）可得得到Poisson方程的九點差分格式:在計算機上實現(xiàn)（4）式，需要在五點差分格式的基礎(chǔ)上在等式兩端分別增加一部分，將等式左側(cè)新增的部分寫成緊湊格式，有:對于該矩陣，可以看成是兩個矩陣的組合:以及則生成這兩個矩陣可以采用Kroncker生成，方法類似于五點差分格式。對于右端添加的關(guān)于f(x,y)的二階導(dǎo)數(shù)，可以采用中心差分格式進行近似代替，即:寫成相應(yīng)的緊湊格式有:該式中的矩陣又可以分解為兩個矩陣的和:三.程序代碼根據(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;%生成網(wǎng)格點與邊界點處f的函數(shù)值，其中f_b1是中間全部為0的向量，f_b2是每一個分塊中間為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);fori=1:N1-1f_b2((i-1)*(N2-1)+1)=f(x(i),ya);f_b2(i*(N2-1))=f(x(i),yb);forj=1:N2-1f_in((i-1)*(N2-1)+j)=f(x(i),y(j));ifi==1f_b1(j)=f(xa,y(j));endifi==N1-1

f_b1((i-1)*(N2-1)+j尸f(xb,y(j));endendendf_in=f_in';%生或迭代矩陣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;%生成迭代矩陣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*h1A2*h2A2)*B;B=B/h2A2;%生成迭代矩陣DS1=spdiags([e2e2],[-11],N1-1,N1-1);S2=spdiags(e3,0,N2-1,N2-1);D=(h1A2+h2A2)/(12*h1A2*h2A2)*kron(S1,K2);%生成f的中心差分格式對應(yīng)矩陣E和FE=kron(S1,S2);K3=spdiags([e3,-4*e3,e3],[-101],N2-1,N2-1);F=kron(I2,K3);f_vec=-f_in-(1/12)*((E+F)*f_in+f_b1+f_b2);%生成網(wǎng)格點[X,Y]=meshgrid(x,y);%計算網(wǎng)格函數(shù)值u_d=(A+B+C+D)\f_vec;u0=zeros(length(f_vec),1);%計算誤差u_real=@(x,y)exp(pi*(x+y))*sin(pi*x).*sin(pi*y);fori=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é)果針對課本P93給出的問題，分別采用步長，將計算出的誤差列表如下：Jk五點差分格式誤差九點差分格式誤差—0.02861330711346713.59173530739554e-06—0.007156344071610482.24435417806035e-07可見米用九點差分格式可以進一步縮小誤差，達到更高階的精度。

五.計算中出現(xiàn)的問題，解決方法及體會在生成九點差分格式的時候，等號右端涉及到了對f的二階偏導(dǎo)，我最初利用符號函數(shù)定義了f,隨后求出具二