MATLAB编程求解二维泊松方程doc资料
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 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