MATLAB编程求解二维泊松方程doc资料

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

M A T L A B编程求解二维

泊松方程

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

相关文档
最新文档