有限差分法的Matlab程序.doc

有限差分法的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=h2/l2; s=2*(1+r); k=1; % 应用Gauss-Seidel法求解差分方程 while knorm; norm=abs(u(i,m-1)-z); end u(i,m-1)=z; end % 对右上角的网格点进行处理 z=(-h2*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=(-h2*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=(-h2*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=(-h2*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)-z); end u(n-1,j)=z; end % 对靠近下边界的网格点进行处理 % 对左下角的网格点进行处理 z=(-h2*fun(x(1),y(1)+gun(a,y(1)+r*gun(x(1),c)+r*u(1,2)+u(2,1)/s; if abs(u(1,1)-z)norm norm=abs(u(1,1)-z); end u(1,1)=z; % 对靠近下边界的除第一点和最后点外网格点进行处理 for i=2:n-2 z=(-h2*fun(x(i),y(1)+r*gun(x(i),c)+r*u(i,2)+u(i+1,1)+u(i-1,1)/s; if abs(u(i,1)-z)norm norm=abs(u(i,1)-z); end u(i,1)=z; end % 对右下角的网格点进行处理 z=(-h2*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; if abs(u(n-1,1)-z)norm norm=abs(u(n-1,1)-z); end u(n-1,1)=z; % 结果输出 if norm=tol fid = fopen(FDresult.txt, wt); fprintf(fid,n*用有限差分法求解矩形域上Poisson方程的输出结果*nn); fprintf(fid,迭代次数: %d次nn,k); fprintf(fid, x的值 y的值 u的值 u的真实值 |u-u(x,y)|n); for i=1:n-1 for j=1:m-1 fprintf(fid, %8.3f %8.3f %14.8f %14.8f %14.8fn, x(i),y(j),u(i,j),gun(x(i),y(j),abs(u(i,j)-gun(x(i),y(j); end end fclose(fid); break; % 用来结束while循环 end k=k+1; end if k=N+1 fid = fopen(FDresult.txt, wt); fprintf(fid,超过最大迭代次数，求解失败！); fclose(fid); endclca1 a2 a3 a4 = textread(F:aa.txt,%f %f %f %f);a = a1 a2 a3;a=a;b=a4;pa,mina,maxa,pb,minb,maxb=premnmx(a,b);net =newrb(pa,pb,0,1.3,24,2);an =sim(net,pa);E = an - pb;m =sse(E)n = mse(E)f1 f2 f3 f4= textread(F:bb.txt,%f %f %f %f);f = f1 f2 f3