matlab求解二维泊松方程,MATLAB编程求解二维泊松方程

MATLAB编程求解二维泊松方程

%%%% 真解 u=sin(pi*x)*sin(pi*y) %%% %%%% 方程 -Laplace(u)=f %%%%%% %%%% f=2*pi^2*sin(pi*x)*sin(pi*y) %%%%%% %%%%difference code for elliptic equations with constant coefficient %%%%% %clear all %clc N=20; h=1/N; S=h^2; 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/h^2; A(i,i+1)=-1/h^2; A(i,i+(N-1))=-1/h^2; end for i=N-1 A(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^2 end for i=(N-2)*(N-1)+1 A(i,i-(N-1))=-1/h^2; A(i,i)=4/h^2; A(i,i+1)=-1/h^2; end for i=(N-1)^2 A(i,i-(N-1))=-1/h^2; A(i,i)=4/h^2; A(i,i-1)=-1/h^2; end for n=2:N-2 i=(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; end for i=2:N-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; end for m=1:N-3 i=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; end for m=2:N-2 i=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 % for m=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; % end for m=1:N-3 for n=2:N-2 i=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; end end %%% Right term F=zeros((N-1)^2,1); for m=0:N-2 for n=1:N-1 i=m*(N-1)+n; F(i)=2*pi^2*sin(pi*n*h)*sin(pi*(m+1)*h); end end %U=zeros((N-1)^2,1); u1=A\F; 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)); end U=zeros(N+1,N+1); for m=1:N+1 U(m,:)=u((m-1)*(N+1)+1:m*(N+1)); end surf(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); end end U_exact=reshape(u_exact,N+1,N+1); subplot(1,2,) err=max(abs(u-u_exact)); l2_err=norm(u-u_exact)*h; err l2_err