二阶椭圆偏微分方程的差分算法-适用于求解矩形区域-任意边界条件-MATLAB源程序

functionlh2(p,q,f,alpha,xspan,yspan,M,N)%¶þ½×ÍÖÔ²ÐÍ·½³Ì±ßÖµÎÊÌâµÄ²î·Ö¸ñʽ%-[d(p*du/dx)/dx+d(p*du/dy)/dy]+q*u=f;%u(x,y)=alpha(x,y)every(x,y)belongstod(xspan*yspan)%p,q,f¾ùΪ¹ØÓÚx,yµÄ¶þÔªº¯Êý%test:u(x,y)=x^2+y^2;f(x,y)=-4+x^2+y^2;p=@(x,y)1;q=@(x,y)0;f=@(x,y)4-x.^2-y.^2;xspan=[-1,1];yspan=[-1,1];alpha=@(x,y)0;M=30;N=10;h1=(xspan(2)-xspan(1))/M;x=xspan(1):h1:xspan(2);h2=(yspan(2)-yspan(1))/N;y=yspan(1):h2:yspan(2);u=zeros(M+1,N+1);H=zeros((M-1)*(N-1),(M-1)*(N-1));g=zeros(M-1,N-1);u(1,:)=alpha(x(1),y(:));u(M+1,:)=alpha(x(M+1),y(:));u(:,1)=alpha(x(:),y(1));u(:,N+1)=alpha(x(:),y(N+1));uh=u([2:M],[2:N]);fori=1:M-1,g(i,:)=f(x(i),y(2:N));end;[x,y]=meshgrid(x,y);x=x';y=y';u=u';uh=uh';g=g';fori=1:(M-1)*(N-1),alpha1=p(x(i)+h1/2,y(i))/h1^2;%alpha1alpha2=p(x(i),y(i)+h2/2)/h2^2;%alpha2alpha3=p(x(i)-h1/2,y(i))/h1^2;%alpha3alpha4=p(x(i),y(i)-h2/2)/h2^2;%alpha4alpha0=(alpha1+alpha2+alpha3+alpha4)+q(round(i/(M-1)),mod(i,M-1));%alpha0ifi1,H(i,i-1)=-alpha3;elseg(i)=g(i)+u(i);end;ifiM-1,H(i,i-M+1)=-alpha4;elseg(i)=g(i)+alpha3*u(i);end;H(i,i)=alpha0;ifi(M-1)*(N-2),H(i,i+M-1)=-alpha2;elseg(i)=g(i)+u(i);end;ifi(M-1)*(N-1)-1,H(i,i+1)=-alpha1;elseg(i)=g(i)+u(i);end;end;uh(:)=H^-1*g(:);x=x';y=y';u=u';uh=uh';g=g';%uh(:)=in_guass_seidel(H,g(:),eye(length(g(:)),1))%Óøß˹ÈûµÂ¶ûµü´úÇó½â%uh(:)=in_jacobi(H,g(:),eye(length(g(:)),1))%ÓÃÑſ˱ȵü´úÇó½âu([2:M],[2:N])=uh;figure;mesh(x,y,u);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Ñſ˱ȵø´úfunctionx=in_jacobi(A,b,x1)exps=1.0e-6;count=0;n=length(x1);forrow=1:n,A(row,[1:row-1,row+1:n])=-A(row,[1:row-1,row+1:n])/A(row,row);b(row,1)=b(row,1)/A(row,row);A(row,row)=0;endx=x1+2*exps;whilemax(abs(x-x1))exps&count100,x=A*x1+b;temp=x1;x1=x;x=temp;count=count+1;end;in_jacobi=x;disp('µþ´ú´ÎÊýÊÇ:');count%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%¸ß˹ÈûµÃ¶ûµü´úfunctionx=in_guass_seidel(A,b,x)exps=1.0e-6;dx=x;count=0;n=length(x);whilemax(abs(dx))exps&count=100,forrow=1:n,ifA(row,row)==0,return;enddx(row)=(b(row)-A(row,[1:row-1])*x([1:row-1])-A(row,[row:n])*x([row:n])).../A(row,row);x(row)=x(row)+dx(row);endcount=count+1;endin_guass_seidel=x;disp('µþ´ú´ÎÊý:');count