有限差分法的Matlab程序(橢圓型方程)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上有限差分法的Matlab程序（橢圓型方程）function FD_PDE(fun,gun,a,b,c,d)    % 用有限差分法求解矩形域上的Poisson方程    tol=10(-6);  % 誤差界    N=1000;  % 最大迭代次數(shù)    n=20;  % x軸方向的網(wǎng)格數(shù)    m=20;  % y軸方向的網(wǎng)格數(shù)    h=(b-a)

2、/n; % x軸方向的步長(zhǎng)    l=(d-c)/m; % y軸方向的步長(zhǎng)    for i=1:n-1        x(i)=a+i*h;    end % 定義網(wǎng)格點(diǎn)坐標(biāo)    for j=1:m-1        y(j)=c+j*l;    end % 定義網(wǎng)格點(diǎn)坐標(biāo)  

3、0; u=zeros(n-1,m-1); %對(duì)u賦初值    % 下面定義幾個(gè)參數(shù)    r=h2/l2;    s=2*(1+r);    k=1;    % 應(yīng)用Gauss-Seidel法求解差分方程    while k<=N        % 對(duì)靠近上邊界的網(wǎng)格點(diǎn)進(jìn)行處理     &#

4、160;      % 對(duì)左上角的網(wǎng)格點(diǎn)進(jìn)行處理            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;            % 對(duì)靠近上邊界的除第一點(diǎn)和最后點(diǎn)外網(wǎng)格點(diǎn)進(jìn)行處理            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、            % 對(duì)右上角的網(wǎng)格點(diǎn)進(jìn)行處理            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 ab

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

10、#160;    % 對(duì)不靠近上下邊界的網(wǎng)格點(diǎn)進(jìn)行處理            for j=m-2:-1:2                % 對(duì)靠近左邊界的網(wǎng)格點(diǎn)進(jìn)行處理           

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

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

13、#160;  % 對(duì)不靠近左右邊界的網(wǎng)格點(diǎn)進(jìn)行處理                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

14、(i+1,j)/s;                    if abs(u(i,j)-z)>norm                       norm=abs(u(i,j)-z)

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

16、60;      end                % 對(duì)靠近右邊界的網(wǎng)格點(diǎn)進(jìn)行處理                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

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

18、;        end                u(n-1,j)=z;            end        % 對(duì)靠近下邊界的網(wǎng)格點(diǎn)進(jìn)行處理   

19、0;        % 對(duì)左下角的網(wǎng)格點(diǎn)進(jìn)行處理            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、對(duì)靠近下邊界的除第一點(diǎn)和最后點(diǎn)外網(wǎng)格點(diǎn)進(jìn)行處理            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、160;      if abs(u(i,1)-z)>norm                  norm=abs(u(i,1)-z);               end     &

23、#160;         u(i,1)=z;            end            % 對(duì)右下角的網(wǎng)格點(diǎn)進(jìn)行處理            z=(-h2*fun(x(n-1),y

24、(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,1)-z);        

25、    end            u(n-1,1)=z;         % 結(jié)果輸出         if norm<=tol             fid = fo

26、pen('FDresult.txt', 'wt');             fprintf(fid,'n*用有限差分法求解矩形域上Poisson方程的輸出結(jié)果*nn');             fprintf(fid,'迭代次數(shù): %d次nn',k);   &#

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