《有限差分法的Matlab程序(椭圆型方程)说课讲解.pdf》由会员分享，可在线阅读，更多相关《有限差分法的Matlab程序(椭圆型方程)说课讲解.pdf（5页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、 有限差分法的 Matlab程序(椭圆型方程)精品文档 收集于网络，如有侵权请联系管理员删除 有限差分法的 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=z
2、eros(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;%对不靠近上下边界的网格
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);end u(i,j)=z;en
4、d%对靠近右边界的网格点进行处理 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
5、,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
6、,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