九点差分格式.doc

二阶椭圆型方程的差分格式题目:: (1)G是xy平面上一有界区域,其边界为分段光滑曲线。在上满足下列边值条件之一: (第一边值条件、强制边值条件), (第二边值条件), (第三边值条件、混合边值条件),,,,,及都是连续函数,。(3)、(4)相加得:又因为因此舍去截断误差,便得到逼近Poisson方程的九点差分格式:其截断误差的阶为化其形式为=同样将上述方程化成如下形式:其中:另外,。我们用来表示数值解满足的线性方程组,则有:(具有三对角矩阵的特征),其中:,,,,。(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);en