MATLAB编程求解二维泊松方程

%%%%真解u=sin(pi*x)*sin(pi*y)%%%%%%%方程-Laplace(u)=f%%%%%%%%%%f=2*pi^2*sin(pi*x)*sin(pi*y)%%%%%%%%%%differencecodeforellipticequationswithconstantcoefficient%%%%%%clearall%clcN=20;h=1/N;S=h^2;x=0:h:1;y=0:h:1;%%%StiffmatrixA=zeros((N-1)^2,(N-1)^2);fori=1A(i,i)=4/h^2;A(i,i+1)=-1/h^2;A(i,i+(N-1))=-1/h^2;endfori=N-1A(i,i-1)=-1/h^2;A(i,i)=4/h^2;A(i,2*i)=-1/h^2;%A(i,i+(N-1))=-1/h^2endfori=(N-2)*(N-1)+1A(i,i-(N-1))=-1/h^2;A(i,i)=4/h^2;A(i,i+1)=-1/h^2;endfori=(N-1)^2A(i,i-(N-1))=-1/h^2;A(i,i)=4/h^2;A(i,i-1)=-1/h^2;endforn=2:N-2i=(N-2)*(N-1)+n;A(i,i-(N-1))=-1/h^2;A(i,i-1)=-1/h^2;A(i,i)=4/h^2;A(i,i+1)=-1/h^2;endfori=2:N-2A(i,i-1)=-1/h^2;A(i,i)=4/h^2;A(i,i+1)=-1/h^2;A(i,i+(N-1))=-1/h^2;endform=1:N-3i=m*(N-1)+1;A(i,i-(N-1))=-1/h^2;A(i,i)=4/h^2;A(i,i+1)=-1/h^2;A(i,i+(N-1))=-1/h^2;endform=2:N-2i=m*(N-1);A(i,i-(N-1))=-1/h^2;A(i,i-1)=-1/h^2;A(i,i)=4/h^2;A(i,i+(N-1))=-1/h^2;end%form=1:N-3%i=m*(N-1)+(N-1);%A(i,i-(N-1))=-1/h^2;%A(i,i-1)=-1/h^2;%A(i,i)=4/h^2;%A(i,i+(N-1))=-1/h^2;%endform=1:N-3forn=2:N-2i=m*(N-1)+n;A(i,i-(N-1))=-1/h^2;A(i,i-1)=-1/h^2;A(i,i)=4/h^2;A(i,i+1)=-1/h^2;A(i,i+(N-1))=-1/h^2;endend%%%RighttermF=zeros((N-1)^2,1);form=0:N-2forn=1:N-1i=m*(N-1)+n;F(i)=2*pi^2*sin(pi*n*h)*sin(pi*(m+1)*h);endend%U=zeros((N-1)^2,1);u1=A\F;u=zeros((N+1)^2,1);form=1:N-1u(m*(N+1)+2:m*(N+1)+N)=u1((m-1)*(N-1)+1:m*(N-1));endU=zeros(N+1,N+1);form=1:N+1U(m,:)=u((m-1)*(N+1)+1:m*(N+1));endsurf(x,y,U)u_exact=zeros((N+1)^2,1);form=0:Nforn=1:N+1i=m*(N+1)+n;u_exact(i)=sin(pi*n*h)*sin(pi*m*h);endendU_exact=reshape(u_exact,N+1,N+1);subplot(1,2,)err=max(abs(u-u_exact));l2_err=norm(u-u_exact)*h;errl2_err