PDA

查看完整版本 : [MATLAB数学相关] 求助差分算法求解椭圆方程,程序没出错但结果有误差


ivanxia
2012-04-29, 11:27
程序是用5点差分求解椭圆方程,程序可以运行,但画出的三维图跟实际有差别。程序如下:
function Y=zploy(L1,L2,N,M)
dx=L1/(N-1);
dy=L2/(M-1);
% Initialize matrix and right hand side
A=sparse(N*M,N*M);
F=zeros(N*M,1);
% Mapping between (i,j) and position in solution vector
map=zeros(N,M);
map( : ) =1:N*M;
% 二阶差分和一阶差分
mol_xi=[0,-1,0;0,0,0;0,1,0]/2/dx;
mol_xixi=[0,1,0;0,-2,0;0,1,0]/dx^2;
mol_yiyi=[0,0,0;1,-2,1;0,0,0]/dy^2;
mol_yi=[0,0,0;-1,0,1;0,0,0]/dy;
% 内部点
for i=2:N-1
for j=2:M-1
mol=-(mol_xixi+mol_yiyi);
% Insert in A
for i1=-1:1
for j1=-1:1
A(map(i,j),map(i+i1,j+j1))=mol(i1+2,j1+2);
end
end
% RHS is 1
F(map(i,j))=1;
end
end
% Bottom boundary (T=0)
for i=1:N
A(map(i,1),map(i,1))=1;
F(map(i,1))=0;
end
% Left boundary (T=0)
for j=1:M
A(map(1,j),map(1,j))=1;
F(map(1,j))=0;
end
% Top boundary (T=0)
for i=1:N
A(map(i,M),map(i,M))=1;
F(map(i,M))=0;
end
for j=1:M
A(map(N,j),map(N,j))=1;
F(map(N,j))=0;
end

Y=A\F;
Y=reshape(Y,N,M);

再做画三维图检验的时候发现
L1=10;
L2=10;
N=21;
M=21;
zz=tfield(L1,L2,N,M);
h=L1/(N-1);
h2=L2/(M-1);
x=0:h:L1;
y=0:h2:L2;
[xx,yy]=meshgrid(x,y);
surf(xx,yy,zz);

画出的图不在中心,不知是哪里出了问题?
是矩阵初始化的问题还是标记?求指教,谢谢。