有限差分法的Matlab程序(椭圆型方程)

有限差分法的Matlab程序（椭圆型方程）functionFD_PDE(fun,gun,a,b,c,d)%用有限差分法求解矩形域上的Poisson方程tol=10^(-6);%误差界N=1000;%最大迭代次数n=20;%x轴方向的网格数m=20;%y轴方向的网格数h=(b-a)/n;%x轴方向的步长l=(d-c)/m;%y轴方向的步长fori=1:n-1x(i)=a+i*h;end%定义网格点坐标forj=1:m-1y(j)=c+j*l;end%定义网格点坐标u=zeros(n-1,m-1);%对u赋初值%下面定义几个参数r=h^2/l^2;s=2*(1+r);k=1;%应用Gauss-Seidel法求解差分方程whilek=N%对靠近上边界的网格点进行处理%对左上角的网格点进行处理z=(-h^2*fun(x(1),y(m-1))+gun(a,y(m-1))+r*gun(x(1),d)+r*u(1,m-2)+u(2,m-1))/s;norm=abs(z-u(1,m-1));u(1,m-1)=z;%对靠近上边界的除第一点和最后点外网格点进行处理fori=2:n-2z=(-h^2*fun(x(i),y(m-1))+r*gun(x(i),d)+r*u(i,m-2)+u(i+1,m-1)+u(i-1,m-1))/s;ifabs(u(i,m-1)-z)norm;norm=abs(u(i,m-1)-z);endu(i,m-1)=z;end%对右上角的网格点进行处理z=(-h^2*fun(x(n-1),y(m-1))+gun(b,y(m-1))+r*gun(x(n-1),d)+r*u(n-1,m-2)+u(n-2,m-1))/s;ifabs(u(n-1,m-1)-z)normnorm=abs(u(n-1,m-1)-z);endu(n-1,m-1)=z;%对不靠近上下边界的网格点进行处理forj=m-2:-1:2%对靠近左边界的网格点进行处理z=(-h^2*fun(x(1),y(j))+gun(a,y(j))+r*u(1,j+1)+r*u(1,j-1)+u(2,j))/s;ifabs(u(1,j)-z)normnorm=abs(u(1,j)-z);endu(1,j)=z;%对不靠近左右边界的网格点进行处理fori=2:n-2z=(-h^2*fun(x(i),y(j))+u(i-1,j)+r*u(i,j+1)+r*u(i,j-1)+u(i+1,j))/s;ifabs(u(i,j)-z)normnorm=abs(u(i,j)-z);endu(i,j)=z;end%对靠近右边界的网格点进行处理z=(-h^2*fun(x(n-1),y(j))+gun(b,y(j))+r*u(n-1,j+1)+r*u(n-1,j-1)+u(n-2,j))/s;ifabs(u(n-1,j)-z)normnorm=abs(u(n-1,j)-z);endu(n-1,j)=z;end%对靠近下边界的网格点进行处理%对左下角的网格点进行处理z=(-h^2*fun(x(1),y(1))+gun(a,y(1))+r*gun(x(1),c)+r*u(1,2)+u(2,1))/s;ifabs(u(1,1)-z)normnorm=abs(u(1,1)-z);endu(1,1)=z;%对靠近下边界的除第一点和最后点外网格点进行处理fori=2:n-2z=(-h^2*fun(x(i),y(1))+r*gun(x(i),c)+r*u(i,2)+u(i+1,1)+u(i-1,1))/s;ifabs(u(i,1)-z)normnorm=abs(u(i,1)-z);endu(i,1)=z;end%对右下角的网格点进行处理z=(-h^2*fun(x(n-1),y(1))+gun(b,y(1))+r*gun(x(n-1),c)+r*u(n-1,2)+u(n-2,1))/s;ifabs(u(n-1,1)-z)normnorm=abs(u(n-1,1)-z);endu(n-1,1)=z;%结果输出ifnorm=tolfid=fopen('FDresult.txt','wt');fprintf(fid,'\n********用有限差分法求解矩形域上Poisson方程的输出结果********\n\n');fprintf(fid,'迭代次数:%d次\n\n',k);fprintf(fid,'x的值y的值u的值u的真实值|u-u(x,y)|\n');fori=1:n-1forj=1:m-1fprintf(fid,'%8.3f%8.3f%14.8f%14.8f%14.8f\n',[x(i),y(j),u(i,j),gun(x(i),y(j)),abs(u(i,j)-gun(x(i),y(j)))]);endendfclose(fid);break;%用来结束while循环endk=k+1;endifk==N+1fid=fopen('FDresult.txt','wt');fprintf(fid,'超过最大迭代次数，求解失败！');fclose(fid);end