有限差分法的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 轴方向的步长 fori=1:n-1 x(i)=a+i*h;

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-1,m-1); %对u赋初值

% 下面定义几个参数

r=h^2/l^2;

s=2*(1+r);

k=1;

% 应用Gauss-Seidel法求解差分方程

while k%26lt;=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));

u(1,m-1)=z;

% 对靠近上边界的除第一点和最后点外网格点进行处理

for i=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;

if abs(u(i,m-1)-z)%26gt;norm;

norm=abs(u(i,m-1)-z);

end

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;

if abs(u(n-1,m-1)-z)%26gt;norm

norm=abs(u(n-1,m-1)-z);

end

u(n-1,m-1)=z;