资源描述
有限差分法的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-1x(i)=a+i*h;end %定义网格点坐标for j=1:m-1 y(j)=c+j*l;end %定义网格点坐标u=zeros(nT,m-l); %对 u 赋初值%下面定义几个参数r二h2/l2;s=2*(1+r);k=1;%应用Gauss-Seidel法求解差分方程while knorm; norm=abs(u(i,m1)z);end u(i,m1)=z;end%对右上角的网格点进行处理z=(-h“2*fun(x(n-l),y(m-l)+gun(b,y(m-l)+r*gun(x(n-l),d)+r*u(nT,m- 2)+u(n2,m1)/s;if abs(u(n1,m1)z)normnorm=abs(u(n1,m1)z);endu(n1,m1)=z;% 对不靠近上下边界的网格点进行处理for j=m-2:-1:2% 对靠近左边界的网格点进行处理z=(-h2*fun(x(l),y(j)+gun(a,y(j)+r*u(l,j+l)+r*u(l,j-l)+u(2,j)/s;if abs(u(1,j)-z)norm norm=abs(u(1,j)-z);endu(1,j)=z;% 对不靠近左右边界的网格点进行处理for i=2:n-2z=(-h2*fun(x(i),y(j)+u(i-l,j)+r*u(i,j+l)+r*u(i,j-1)+u(i+1,j)/s;if abs(u(i,j)-z)norm norm=abs(u(i,j)-z);endu(i,j)=z;end% 对靠近右边界的网格点进行处理z=(-h“2*fun(x(n-l),y(j)+gun(b,y(j)+r*u(nT,j+l)+r*u(nT,j-1)+u(n-2,j)/s;Jif abs(u(n-1,j)-z)norm norm=abs(u(n-1,j)-z);endu(n-1,j)=z;end% 对靠近下边界的网格点进行处理%对左下角的网格点进行处理z=(-h2*fun(x(l),y(l)+gun(a,y(l)+r*gun(x(l),c)+r*u(l,2)+u(2,l)/sif abs(u(1,1)-z)norm norm=abs(u(1,1)-z);endu(1,1)=z; %对靠近下边界的除第一点和最后点外网格点进行处理 for i=2:n-2z=(h2*fun(x(i),y(l)+r*gun(x(i),c)+r*u(i,2)+u(i+l,l)+u(iT,1)/s;if abs(u(i,1)-z)normnorm=abs(u(i,1)z);endu(i,1)=z;end%对右下角的网格点进行处理页脚内容2z=(-h“2*fun(x(n-l),y(l)+gun(b,y(l)+r*gun(x(n-l),c)+r*u(nT,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=tolfid = 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-1for j=1：m-1fprintf(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);endend fclose(fid);break;%用来结束while循环endk=k+1;end if k=N+1fid = fopen(FDresult.txt, wt);fprin tf(fid,超过最大迭代次数，求解失败！)； fclose(fid);end clca1 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;f=f;pf = tramnmx(f,mina,maxa);an2 = sim(net,pf);g =postmnmx(an2,minb,maxb);g= g;E2 = g- f4;mm =sse(E2)nn = mse(E2)
展开阅读全文