九点差分格式.docx

1、九点差分格式二阶椭圆型方程的差分格式题目：九点差分格式1.考虑问题 考虑Poisson 方程: (1)G是xy平面上一有界区域，其边界为分段光滑曲线。在上满足下列边值条件之一： (第一边值条件、强制边值条件)， (第二边值条件)， (第三边值条件、混合边值条件)，及都是连续函数，。2.九点差分格式 2.1公式推导因为有 将（3）、（4）相加得：又因为因此舍去截断误差，便得到逼近Poisson方程的九点差分格式：其截断误差的阶为化其形式为 =同样将上述方程化成如下形式：其中：另外，。我们用来表示数值解满足的线性方程组，则有：(具有三对角矩阵的特征)，其中：，。2.2差分方程数值例子的求解考虑问题

2、 其精确解2.3差分方程数值例子求解考虑问题 其精确解2.4求解图像3.附录f(x,y)函数function f=ff(x,y)f=(pi*pi-1)*exp(x)*sin(pi*y);u(x,y)函数function f=fu(x,y)f=sin(pi*y)*exp(x);clear;clc;format short;%数据准备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);%对自变量进行赋值for i=1:M x(i)=x0+i*Hx; ux0(i)=

3、fu(x(i),0); ux1(i)=fu(x(i),1);endfor j=1:N y(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*

4、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(dia

5、g(ones(1,M-1),-1),J);max=maxA+maxB+maxC;for ii=1:M for jj=1:N f(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); end if(jj=N) f(ii-1)*N+jj)= f(ii-1)*N+jj)-a*u0y(jj-1)-

6、b*u0y(jj)-c*fu(x0,yn)-f1*ux1(ii)-m*ux1(ii+1); %jj end if(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); end end if(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); end if(jj=N) f(ii-1)*N+jj)= f(ii-1)*N+jj)-c*ux1(ii-1)-f1

7、*ux1(ii)-g*u2y(jj-1)-h*u2y(jj)-m*fu(xn,yn); end if(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); end end if(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); end if(jj=N) f(ii-1)*N+jj)= f(ii-1)*N+jj)-c*ux1(ii-1)-f1*ux1(ii)-m*ux1(ii+1); end end endendjg=(inv(max)*f);jg;jjgg=reshape(jg,M,N);for iii=1:M for jjj=1:N JJGG(iii,jjj)=fu(x(iii),y(jjj); endendX,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);