有限差分法Matlab程序[椭圆型方程].docVIP

有限差分法的Matlab程序（椭圆型方程） function FD_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轴方向的步长??? for i=1:n-1??????? x(i)=a+i*h;??? end % 定义网格点坐标??? for j=1:m-1??????? y(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法求解差分方程??? while k=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;??????????? % 对靠近上边界的除第一点和最后点外网格点进行处理??????????? for i=2:n-2??????????????? z=(-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;??????????????? if abs(u(i,m-1)-z)norm;?????????????????? norm=abs(u(i,m-1)-z);??????????????? end??????????????? u(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;??????????? if abs(u(n-1,m-1)-z)norm?????????????? norm=abs(u(n-1,m-1)-z);??????????? end??????????? u(n-1,m-1)=z;??????? % 对不靠近上下边界的网格点进行处理??????????? for j=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;??????????????? if abs(u(1,j)-z)norm?????????????????? norm=abs(u(1,j)-z);??????????????? end??????????????? u(1,j)=z;??????????????? % 对不靠近左右边界的网格点进行处理??????????????? for i=2:n-2??????????????????? z=(-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;??????????????????? if abs(u(i,j)-z)norm?????????????????????? norm=abs(u(i,j)-z);??????????????????? end??????????????????? u(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;??????????????? if abs(u(n-1,j)-z)norm?????????????????? norm=abs(u(n-1,j)-