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

%%%% 真解 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

本文来源：https://www.bwwdw.com/article/5jze.html