有限差分法的Matlab程序(椭圆型方程)3页

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轴方向的步长

2、    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-

3、1,m-1); %对u赋初值    % 下面定义几个参数    r=h2/l2;    s=2*(1+r);    k=1;    % 应用Gauss-Seidel法求解差分方程    while k<=N        % 对靠近上边界的网格点进行处理       

4、60;    % 对左上角的网格点进行处理            z=(-h2*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);      &#

5、160;     u(1,m-1)=z;            % 对靠近上边界的除第一点和最后点外网格点进行处理            for i=2:n-2               

6、 z=(-h2*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;                if abs(u(i,m-1)-z)>norm;                  

7、norm=abs(u(i,m-1)-z);                end                u(i,m-1)=z;            end  &

8、#160;         % 对右上角的网格点进行处理            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)-

9、z)>norm               norm=abs(u(n-1,m-1)-z);            end            u(n-1,m-1)=z;     &#

10、160;  % 对不靠近上下边界的网格点进行处理            for j=m-2:-1:2                % 对靠近左边界的网格点进行处理              

11、;  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=

12、abs(u(1,j)-z);                end                u(1,j)=z;                %

13、 对不靠近左右边界的网格点进行处理                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;

14、60;                   if abs(u(i,j)-z)>norm                       norm=abs(u(i,j)-z);  

15、                  end                    u(i,j)=z;           

16、0;    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; &

17、#160;              if abs(u(n-1,j)-z)>norm                   norm=abs(u(n-1,j)-z);          

18、      end                u(n-1,j)=z;            end        % 对靠近下边界的网格点进行处理      

19、;      % 对左下角的网格点进行处理            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     &#

20、160;         norm=abs(u(1,1)-z);            end            u(1,1)=z;            % 对靠近下边界的除第一点和最

21、后点外网格点进行处理            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;          

22、60;    if abs(u(i,1)-z)>norm                  norm=abs(u(i,1)-z);               end       &#

23、160;       u(i,1)=z;            end            % 对右下角的网格点进行处理            z=(-h2*fun(x(n-1),y(1)+gun(b,y(1

24、)+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);          &

25、#160; end            u(n-1,1)=z;         % 结果输出         if norm<=tol             fid = fopen('FDre

26、sult.txt', 'wt');             fprintf(fid,'n*用有限差分法求解矩形域上Poisson方程的输出结果*nn');             fprintf(fid,'迭代次数: %d次nn',k);     

27、60;       fprintf(fid,'    x的值    y的值       u的值           u的真实值      |u-u(x,y)|n');             for i=1:n-1                 for j=1:m-1