poissonfd.m

sareli<matlab科学计算>配套程序

function [u,x,y,error]=poissonfd(a,c,b,d,nx,ny,f,g,uex)
%POISSONFD two-dimensional Poisson solver
%   [U,X,Y]=POISSONFD(A,C,B,D,NX,NY,FUN,BOUND) solves with
%   the 5-points finite difference scheme the problem
%   -LAPL(U) = FUN in the rectangle (A,C)X(B,D) with
%   Dirichlet boundary conditions U(X,Y)=BOUND(X,Y) for any
%   (X,Y) on the boundary on the rectangle. 
%
%   [U,X,Y,ERROR]=POISSONFD(A,C,B,D,NX,NY,FUN,BOUND,UEX) computes
%   also the maximum nodal error ERROR with respect to the exact
%   solution UEX.
%   FUN, BOUND and UEX can be online functions.
if nargin == 8 & nargout==4
    uex = '0 + 0.*x'; 
end 
nx=nx+1;         ny=ny+1;
hx = (b-a)/nx;   hy = (d-c)/ny; 
nx1 = nx+1;      hx2 = hx^2; 
hy2 = hy^2; 
kii = 2/hx2+2/hy2;    kix = -1/hx2;    kiy = -1/hy2;
dim = (nx+1)*(ny+1); 
K = speye(dim,dim); 
rhs = zeros(dim,1); 
y = c;
for m = 2:ny
    x = a; y = y + hy;
    for n = 2:nx
        i = n+(m-1)*(nx+1); 
        x = x + hx; 
        rhs(i) = feval(f,x,y);  
        K(i,i) = kii; 
        K(i,i-1) = kix; 
        K(i,i+1) = kix;
        K(i,i+nx1) = kiy; 
        K(i,i-nx1) = kiy;
    end
end 
rhs1 = zeros(dim,1);
x = [a:hx:b]; 
rhs1(1:nx1) = feval(g,x,c);
rhs1(dim-nx:dim) = feval(g,x,d);
y = [c:hy:d]; 
rhs1(1:nx1:dim-nx) = feval(g,a,y); 
rhs1(nx1:nx1:dim) = feval(g,b,y);
rhs = rhs - K*rhs1; 
nbound = [[1:nx1],[dim-nx:dim],[1:nx1:dim-nx],[nx1:nx1:dim]];
ninternal = setdiff([1:dim],nbound);
K = K(ninternal,ninternal); 
rhs = rhs(ninternal); 
utemp = K\rhs; 
uh = rhs1; 
uh (ninternal) = utemp; 
k = 1; y = c;
for j = 1:ny+1
    x = a;
    for i = 1:nx1
        u(i,j) = uh(k); 
        k = k + 1; 
        ue(i,j) = feval(uex,x,y); 
        x = x + hx; 
    end
    y = y + hy;
end
x = [a:hx:b]; 
y = [c:hy:d];
if nargout == 4, if nargin == 8
        error('Exact solution not available');
    else
        error = max(max(abs(u-ue)))/max(max(abs(ue)));
end, end
return