程序有错：二阶椭圆型偏微分方程的五点差分格式

chentaocuc

我想实现二阶椭圆偏微分方程的五点差分格式的程序，但是程序编好了以后总是和真解有出入。我把程序贴在这里，请大家帮我看看问题出在哪里？
function lh2(p,q,f,alpha,xspan,yspan,M,N)
% 二阶椭圆型方程边值问题的差分格式
% -[d(p*du/dx)/dx+d(p*du/dy)/dy]+q*u=f;
% u(x,y)=alpha(x,y) every (x,y) belongs to [xspan*yspan]
% p,q,f均为关于x,y的二元函数
% 2007-7-3
% test: u(x,y)=(x^2-1)*(y^2-1);
% 真解: [x,y]=meshgrid(-1:0.1:1,-1:0.1:1);z=(x.^2-1).*(y.^2-1);mesh(x,y,z);
%调试部分的初始化
p=@(x,y) 1;%偏系数
q=@(x,y) 0; %u系数
f=@(x,y) 4-2*x.^2-2*y.^2;%方程右边
xspan=[-1,1];yspan=[-1,1];%矩形区域
alpha=@(x,y) 0;%边界
M=10;N=10;%剖分可以不同

% 程序部分
%分别对x和y方向进行剖分
tic;
h1=(xspan(2)-xspan(1))/M;
x=xspan(1):h1:xspan(2);
h2=(yspan(2)-yspan(1))/N;
y=yspan(1):h2:yspan(2);
u=zeros(M+1,N+1);%初始化函数函数值矩阵
H=zeros((M-1)*(N-1));%初始化H矩阵
g=zeros(M-1,N-1);%初始化方程右部g
%根据边界函数α(x,y)计算u的边界值
u(1,:)=alpha(x(1),y(:));
u(M+1,:)=alpha(x(M+1),y(:));
u(:,1)=alpha(x(:),y(1));
u(:,N+1)=alpha(x(:),y(N+1));
uh=u([2:M],[2:N]);%取u的内点,即网格内点
for i=1:M-1,
g(i,:)=f(x(i),y([2:N]));%第一遍初始化g
end;
[x,y]=meshgrid(x,y);%真正网格化x,y区域
x=x'';
y=y'';
u=u'';
uh=uh'';
g=g'';%将画图的节点排序转化为算法的排序
for i=1:(M-1)*(N-1),
alpha1=p(x(i)+h1/2,y(i))/h1^2;%alpha1
alpha2=p(x(i),y(i)+h2/2)/h2^2;%alpha2
alpha3=p(x(i)-h1/2,y(i))/h1^2;%alpha3
alpha4=p(x(i),y(i)-h2/2)/h2^2;%alpha4
alpha0=(alpha1+alpha2+alpha3+alpha4)+q(x(i),y(i));%alpha0
if i>1,
H(i,i-1)=-alpha3;
else
g(i)=g(i)+alpha3*u(i);
end;
if i>M-1,
H(i,i-M+1)=-alpha4;
else
g(i)=g(i)+alpha4*u(i);
end;
H(i,i)=alpha0;
if i<(M-1)*(N-2),
H(i,i+M-1)=-alpha2;
else
g(i)=g(i)+alpha2*u(i);
end;
if i<(M-1)*(N-1)-1,
H(i,i+1)=-alpha1;
else
g(i)=g(i)+alpha1*u(i);
end;
end;
uh(:)=H^(-1)*g(:);
%uh(:)=inner_guass_seidel(H,g(:),eye(length(g(:)),1));
%uh(:)=inner_jacobi(H,g(:),eye(length(g(:)),1))
u=u'';
uh=uh'';
u([2:M],[2:N])=uh;
x=x'';
y=y'';%将算法的节点排序转化为画图的排序
toc;
figure;
mesh(x,y,u'');
xlabel(''x'');ylabel(''y'');zlabel(''u'');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%雅克比跌代
function x=inner_jacobi(A,b,x1)
exps=1.0e-6;
count=0;
n=length(x1);
for row=1:n,
A(row,[1:row-1,row+1:n])=-A(row,[1:row-1,row+1:n])/A(row,row);
b(row,1)=b(row,1)/A(row,row);
A(row,row)=0;
end
x=x1+2*exps;
while max(abs(x-x1))>exps,
x=A*x1+b;
temp=x1;
x1=x;
x=temp;
count=count+1;
end;
in_jacobi=x;
disp(''叠代次数是:'');
count
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%高斯塞得尔迭代
function x=inner_guass_seidel(A,b,x)
exps=1.0e-6;
dx=x;
count=0;
n=length(x);
while max(abs(dx))>exps,
for row=1:n,
if A(row,row)==0,
return;
end
dx(row)=(b(row)-A(row,[1:row-1])*x([1:row-1])-A(row,[row:n])*x([row:n]))...
/A(row,row);
x(row)=x(row)+dx(row);
end
count=count+1;
end
in_guass_seidel=x;
disp(''叠代次数:'');
count