用五点有限差分格式求解椭圆型方程(偏微分方程) 程序3.doc

用五点有限差分格式求解椭圆型方程（偏微分方程） 程序32010-06-16 06:55function varargout=liu(varargin)a=0;b=2;c=0;d=1;h1=1/32;h2=1/32;x=a:h1:b;y=c:h2:d;m=length(x);n=length(y);C=tri1(-1/h12,2*(1/h12+1/h22),-1/h12,n-2);D=-1/h12*eye(n-2);AA=tri2(D,C,D,m-2);BB=fc1t(a,b,c,d,h1,h2);XX=AABB;UU=zeros(m,n);for r=1:m for j=1:n UU(r,j)=fU(x(r),y(j); endendUL=UU;for r=2:m-1 UL(r,2:n-1)=XX(r-2)*(n-2)+1:(r-1)*(n-2);endY,X=meshgrid(y,x);Z=abs(UL-UU);mesh(Y,X,Z);shading flat;xlabel(X,FontSize,14);ylabel(Y,FontSize,14);zlabel(error,FontSize,14);title(误差图)%求解下问题%-(u_xx+u_yy)=f(x,y) x,y 在区域内 x in axb,cyd%u=g(x,y) x,y在边界上%u=g(a,y)=g1y u=g(b,y)=g2y c=y=d%u=g(x,c)=g1x u=g(x,d)=g2x axb%h1离散x方向的步长%h2离散y方向的步长function c1=fc1t(a,b,c,d,h1,h2)x=a:h1:b;y=c:h2:d;m=length(x);n=length(y);c1=zeros(n-2)*(m-2),1);c1(1)=-1/h12*g1y(y(2)-1/h22*g1x(x(2);for r=2:n-2 c1(r)=-1/h12*g1y(y(r+1);endfor r=2:m-3 c1(r-1)*(n-2)+1)=-1/h22*g1x(x(r+1); c1(r*(n-2)=-1/h22*g2x(x(r+1);endc1(m-3)*(n-2)+1)=-1/h22*g1x(x(m-1)-1/h12*g2y(y(2);for r=2:n-2 c1(m-3)*(n-2)+r)=-1/h12*g2y(y(r+1);endc1(n-2)*(m-2)=-1/h12*g2y(y(n-1)-1/h22*g2x(x(m-1);for r=2:m-1 for j=2:n-1 c1(r-2)*(n-2)+j-1)=c1(r-2)*(n-2)+j-1)-f(x(r),y(j); endendc1=-c1;function y=tri1(A,B,C,n)y=diag(B*ones(n,1)+diag(C*ones(n-1,1),1)+diag(A*ones(n-1,1),-1);function y=tri2(A,B,C,n)n1,n2=size(A);y=zeros(n*n1,n*n1);for j=2:n-1 y(j-1)*n1+1:j*n1,(j-2)*n1+1:(j-1)*n1)=A; y(j-1)*n1+1:j*n1,(j-1)*n1+1:j*n1)=B; y(j-1)*n1+1:j*n1,j*n1+1:(j+1)*n1)=C;endy(1:n1,1:n1)=B;y(1:n1,n1+1:2*n1)=C;y(n-1)*n1+1:n*n1,(n-1)*n1+1:n*n1)=B;y(n-1)*n1+1:n*n1,(n-2)*n1+1:(n-1)*n1)=A;function z=fU(x,y)%精确解函数z=exp(x)*sin(pi*y);function z=f(x,y)z=(pi2-1)*exp(x)*sin(pi*y);function z=g1x(x)z=0;