《有限差分法的Matlab程序(椭圆型方程).doc》由会员分享，可在线阅读，更多相关《有限差分法的Matlab程序(椭圆型方程).doc（3页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精品文档，仅供学习与交流，如有侵权请联系网站删除有限差分法的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-
2、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; % 对不靠近上
3、下边界的网格点进行处理 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); e
4、nd 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(
5、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
6、=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); end【精品文档】第 3 页