九点差分格式

二阶椭圆型方程的差分格式题目：九点差分格式1.考虑问题考虑Poisson方程:-(u+u)=(,),(,)xxyyufxyxyG(1)G是xy平面上一有界区域，其边界为分段光滑曲线。在上u满足下列边值条件之一：|(,)uxy(第一边值条件、强制边值条件)，|(,)uxyn(第二边值条件)，|(,)ukuxyn(第三边值条件、混合边值条件)，(,)fxy，(,)xy，(,)xy，(,)xy，及(,)kxy都是连续函数，0k。2.九点差分格式2.1公式推导因为有1,1,,1,1221222ijijijijijijhijijuuuuuuufhh将（3）、（4）相加得：44224124422422222241212222222221(,)(,)1(,)(,)()12(,)(,)(,)1(,)()1212,1,12ijijhijijijijijijijijuxyuxyuxyuxyhhOhxyuxyuxyuxyhhuxyhhOhxyxyxyfxyfxyhx2422241222222,(,)()12ijijfxyuxyhhhOhyxy又因为4''''''112222221111221211211111(,)(,)2(,)(,)()1[(,)2(,)(,)2((,)2(,)(,))(,)2(,)(,)]()ijxxijxxijxxijijijijijijijijijijuxyuxyuxyuxyOhxyhuxyuxyuxyhhuxyuxyuxyuxyuxyuxyOh因此111111111122221112122222412221(,)[4(,)2((,)(,)(,)12(,))(,)(,)(,)(,)]()/(,)(,)1(,)()()12hijijijijijijijijijijijijijuxyuxyuxyuxyuxyuxyuxyuxyuxyuxyhhhhfxyfxyfxyhhOhxy舍去截断误差，便得到逼近Poisson方程的九点差分格式：1,,11,,122221,11,11,11,112122''2''121[42()12]()/1[(,)(,)]12hijijijijijijijijijijijxxijyyijuuuuuuuuuuhhhhfhfxyhfxy其截断误差的阶为4()Oh化其形式为2221,1,1,1,1221212[(2)(2)]ijijijijijijhuuuhuuuhh1,,11,,11,11,11,11,1[42()]ijijijijijijijijijuuuuuuuuu=222212122212121[('''')]12ijxxyyhhfhfhfhh同样将上述方程化成如下形式：1,11,1,1,1,11,11,1,12''2''121((,)(,))12ijijijijijijijijijijxxijyyijaubucudueufuguhumufhfxyhfxy其中：222211121122,2,,2,4,,,,akbkckdkekhhhhfdgchbma另外，2312221212hhkhh。我们用Auf来表示数值解满足的线性方程组，则有：AAIJAIJAA(具有三对角矩阵的特征)，其中：efdAAfde，hmgImgh，bcaJcab，1111[,,,,,,]TNMMNuuuuu，1111[,,,,,,]TNMMNfffff。2.2差分方程数值例子的求解-(u+u)=(,),(,)xxyyufxyxyG考虑问题222222()(1)sin,02,01(0,)sin(2,)sin(,0)0(,1)0xuueyxyxyuyyuyeyuxux其精确解(,)sinxuxyey22222212122212121[()(1)sin]12ixijjhhfhheyhh2.3差分方程数值例子求解-(u+u)=(,),(,)xxyyufxyxyG考虑问题222222()(1)sin,02,01(0,)sin(2,)sin(,0)0(,1)0xuueyxyxyuyyuyeyuxux其精确解(,)sinxuxyey2.4求解图像3.附录f(x,y)函数functionf=ff(x,y)f=(pi*pi-1)*exp(x)*sin(pi*y);u(x,y)函数functionf=fu(x,y)f=sin(pi*y)*exp(x);clear;clc;formatshort;%数据准备M=40;%x轴划分M份N=40;%y轴划分N份x0=0;%x起点xn=2;%x终点y0=0;%y起点yn=1;%y终点Hx=(xn-x0)/(M+1);Hy=(yn-y0)/(N+1);%对自变量进行赋值fori=1:Mx(i)=x0+i*Hx;ux0(i)=fu(x(i),0);ux1(i)=fu(x(i),1);endforj=1:Ny(j)=y0+j*Hy;u0y(j)=fu(0,y(j));u2y(j)=fu(2,y(j));end%a是u(i-1,j)系数b是u(i,j-1)系数c是u(i,j)系数d是u(i,j+1)系数e是u(i+1,j)系数KK=-(Hx*Hx+Hy*Hy)/(12*Hx*Hx*Hy*Hy);a=KK;b=-1/(Hx*Hx)-2*KK;c=KK;d=-1/(Hy*Hy)-2*KK;e=2/(Hx*Hx)+2/(Hy*Hy)+4*KK;f1=-1/(Hy*Hy)-2*KK;g=KK;h=-1/(Hx*Hx)-2*KK;m=KK;A=diag(ones(1,N)*e)+diag(ones(1,N-1)*f1,1)+diag(ones(1,N-1)*d,-1);I=diag(ones(1,N)*h)+diag(ones(1,N-1)*m,1)+diag(ones(1,N-1)*g,-1);J=diag(ones(1,N)*b)+diag(ones(1,N-1)*c,1)+diag(ones(1,N-1)*a,-1);maxA=blkdiag(kron(eye(M),A));maxB=blkdiag(kron(diag(ones(1,M-1),1),I));maxC=blkdiag(kron(diag(ones(1,M-1),-1),J));max=maxA+maxB+maxC;forii=1:Mforjj=1:Nf((ii-1)*N+jj)=ff(x(ii),y(jj))+1/12*(Hx*Hx*uxx(x(ii),y(jj))+Hy*Hy*uyy(x(ii),y(jj)));if(ii==1)if(jj==1)f((ii-1)*N+jj)=f((ii-1)*N+jj)-a*fu(x0,y0)-b*u0y(jj)-c*u0y(jj+1)-d*ux0(ii)-g*ux0(ii+1);endif(jj==N)f((ii-1)*N+jj)=f((ii-1)*N+jj)-a*u0y(jj-1)-b*u0y(jj)-c*fu(x0,yn)-f1*ux1(ii)-m*ux1(ii+1);%jjendif(jj~=1&jj~=N)f((ii-1)*N+jj)=f((ii-1)*N+jj)-a*u0y(jj-1)-b*u0y(jj)-c*u0y(jj+1);endendif(ii==M)if(jj==1)f((ii-1)*N+jj)=f((ii-1)*N+jj)-a*ux0(ii-1)-d*ux0(ii)-g*fu(xn,y0)-h*u2y(jj)-m*u2y(jj+1);endif(jj==N)f((ii-1)*N+jj)=f((ii-1)*N+jj)-c*ux1(ii-1)-f1*ux1(ii)-g*u2y(jj-1)-h*u2y(jj)-m*fu(xn,yn);endif(jj~=1&jj~=N)f((ii-1)*N+jj)=f((ii-1)*N+jj)-g*u2y(jj-1)-h*u2y(jj)-m*u2y(jj+1);endendif(ii~=1&ii~=M)if(jj==1)f((ii-1)*N+jj)=f((ii-1)*N+jj)-a*ux0(ii-1)-d*ux0(ii)-g*ux0(ii+1);endif(jj==N)f((ii-1)*N+jj)=f((ii-1)*N+jj)-c*ux1(ii-1)-f1*ux1(ii)-m*ux1(ii+1);endendendendjg=(inv(max)*f')';jg;jjgg=reshape(jg,M,N);foriii=1:Mforjjj=1:NJJGG(iii,jjj)=fu(x(iii),y(jjj));endend[X,Y]=meshgrid(x,y);subplot(1,2,1);surf(X,Y,jjgg)title('九点格式');subplot(1,2,2);surf(X,Y,JJGG')title('u(x,y)=sin(pi*y)*exp(x)');