MATLAB编程求解二维泊松方程
3页
1、% 真解 u=sin(pi*x)*sin(pi*y) % 方程 -Laplace(u)=f % f=2*pi2*sin(pi*x)*sin(pi*y) %difference code for elliptic equations with constant coefficient %clear all%clcN=20;h=1/N;S=h2;x=0:h:1;y=0:h:1;% Stiff matrix A=zeros(N-1)2,(N-1)2);for i=1A(i,i)=4/h2;A(i,i+1)=-1/h2;A(i,i+(N-1)=-1/h2;endfor i=N-1A(i,i-1)=-1/h2;A(i,i)=4/h2;A(i,2*i)=-1/h2; %A(i,i+(N-1)=-1/h2endfor i=(N-2)*(N-1)+1A(i,i-(N-1)=-1/h2;A(i,i)=4/h2;A(i,i+1)=-1/h2;endfor i=(N-1)2A(i,i-(N-1)=-1/h2;A(i,i)=4/h2;A(i,i-1)=-1/h2;endfor n=2:N-2i=(N-2)*(
2、N-1)+n;A(i,i-(N-1)=-1/h2;A(i,i-1)=-1/h2;A(i,i)=4/h2;A(i,i+1)=-1/h2;endfor i=2:N-2A(i,i-1)=-1/h2;A(i,i)=4/h2;A(i,i+1)=-1/h2;A(i,i+(N-1)=-1/h2;endfor m=1:N-3i=m*(N-1)+1;A(i,i-(N-1)=-1/h2;A(i,i)=4/h2;A(i,i+1)=-1/h2;A(i,i+(N-1)=-1/h2;endfor m=2:N-2i=m*(N-1);A(i,i-(N-1)=-1/h2;A(i,i-1)=-1/h2;A(i,i)=4/h2;A(i,i+(N-1)=-1/h2;end% for m=1:N-3% i=m*(N-1)+(N-1);% A(i,i-(N-1)=-1/h2;% A(i,i-1)=-1/h2;% A(i,i)=4/h2;% A(i,i+(N-1)=-1/h2;% endfor m=1:N-3for n=2:N-2i=m*(N-1)+n;A(i,i-(N-1)=-1/h2;A(i,i-1)=-1/h2;A(i,i)
3、=4/h2;A(i,i+1)=-1/h2;A(i,i+(N-1)=-1/h2;endend% Right termF=zeros(N-1)2,1);for m=0:N-2for n=1:N-1i=m*(N-1)+n;F(i)=2*pi2*sin(pi*n*h)*sin(pi*(m+1)*h);endend%U=zeros(N-1)2,1);u1=AF; u=zeros(N+1)2,1);for m=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);for m=1:N+1U(m,:)=u(m-1)*(N+1)+1:m*(N+1);endsurf(x,y,U)u_exact=zeros(N+1)2,1);for m=0:Nfor n=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
《MATLAB编程求解二维泊松方程》由会员飞***分享，可在线阅读，更多相关《MATLAB编程求解二维泊松方程》请在金锄头文库上搜索。
