资源描述:
《有限差分法的Matlab程序(椭圆型方程).doc》由会员上传分享,免费在线阅读,更多相关内容在学术论文-天天文库。
1、有限差分法的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-1 x(i)=a+i*h; end%定义网格点坐标 forj=1:m-1 y(j)=c+j*l; end%定
2、义网格点坐标 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))
3、; u(1,m-1)=z; %对靠近上边界的除第一点和最后点外网格点进行处理 fori=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; ifabs(u(i,m-1)-z)>norm; norm=abs(u(i,m-1)-z); end
4、 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; ifabs(u(n-1,m-1)-z)>norm norm=abs(u(n-1,m-1)-z); end u(n-1,m-1)=z; %对不靠近上下边界的网格点进行处
5、理 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)>norm norm=abs(u(1,j)-z); end u(1,j)=z; %对不靠近左右边界的网格
6、点进行处理 fori=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; ifabs(u(i,j)-z)>norm norm=abs(u(i,j)-z); end u(i,j)=z; end
7、%对靠近右边界的网格点进行处理 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)>norm norm=abs(u(n-1,j)-z); end u(n-1,j)=z; end %对靠近下边界的网格点进行处理 %对左下角的
8、网格点进行处理 z=(-h^2*fun(x(1),y(1))+gun(a,y(1))+r