「MATLAB编程求解二维泊松方程」.doc

% 真解 u=sin(pi*x)*sin(pi*y) % 方程 -Laplace(u)=f % f=2*pi2*sin(pi*x)*sin(pi*y) %difference code for elliptic equations with constant coefficient %clear all%clcN=20;h=1/N;S=h2;x=0:h:1;y=0:h:1;% Stiff matrix A=zeros(N-1)2,(N-1)2);for i=1 A(i,i)=4/h2; A(i,i+1)=-1/h2; A(i,i+(N-1)=-1/h2;endfor i=N-1 A(i,i-1)=-1/h2; A(i,i)=4/h2; A(i,2*i)=-1/h2; %A(i,i+(N-1)=-1/h2endfor i=(N-2)*(N-1)+1 A(i,i-(N-1)=-1/h2; A(i,i)=4/h2; A(i,i+1)=-1/h2;endfor i=(N-1)2 A(i,i-(N-1)=-1/h2; A(i,i)=4/h2; A(i,i-1)=-1/h2;endfor n=2:N-2 i=(N-2)*(N-1)+n; A(i,i-(N-1)=-1/h2; A(i,i-1)=-1/h2; A(i,i)=4/h2; A(i,i+1)=-1/h2;endfor i=2:N-2 A(i,i-1)=-1/h2; A(i,i)=4/h2; A(i,i+1)=-1/h2; A(i,i+(N-1)=-1/h2;endfor m=1:N-3 i=m*(N-1)+1; A(i,i-(N-1)=-1/h2; A(i,i)=4/h2; A(i,i+1)=-1/h2; A(i,i+(N-1)=-1/h2;endfor m=2:N-2 i=m*(N-1); A(i,i-(N-1)=-1/h2; A(i,i-1)=-1/h2; A(i,i)=4/h2; A(i,i+(N-1)=-1/h2;end% for m=1:N-3% i=m*(N-1)+(N-1);% A(i,i-(N-1)=-1/h2;% A(i,i-1)=-1/h2;% A(i,i)=4/h2;% A(i,i+(N-1)=-1/h2;% endfor m=1:N-3 for n=2:N-2 i=m*(N-1)+n; A(i,i-(N-1)=-1/h2; A(i,i-1)=-1/h2; A(i,i)=4/h2; A(i,i+1)=-1/h2; A(i,i+(N-1)=-1/h2; endend% Right termF=zeros(N-1)2,1);for m=0:N-2 for n=1:N-1 i=m*(N-1)+n; F(i)=2*pi2*sin(pi*n*h)*sin(pi*(m+1)*h); endend%U=zeros(N-1)2,1);u1=AF; u=zeros(N+1)2,1);for m=1:N-1 u(m*(N+1)+2:m*(N+1)+N)=u1(m-1)*(N-1)+1:m*(N-1);endU=zeros(N+1,N+1);for m=1:N+1 U(m,:)=u(m-1)*(N+1)+1:m*(N+1);endsurf(x,y,U)u_exact=zeros(N+1)2,1);for m=0:N for n=1:N+1 i=m*(N+1)+n; u_exact(i)=sin(pi*n*h)*sin(pi*m*h); endendU_exact=reshape(u_exact,N+1,N