Code covered by the BSD License

Environmental Modeling - using MATLAB

Ekkehard Holzbecher (view profile)

Springer Book, 2nd Edition

DiffDecay2D.m
```% solution of 2D diffusion-decay equation using sparse matrix
%
%  E.Holzbecher,    19.3.2011
%
%--------------------------------------------------
nx = 12; ny = 12;                  % dimensions in x- and y-direction
h = 1/4;                           % grid spacing

% boundary type indicators (1=Dirichlet, 0=Neumann no-flow)
ltop = logical(zeros(1,nx));         % top
lbottom = logical(zeros(1,nx));      % bottom
lleft = logical([ones(6,1) zeros(6,1)]);       % left
lright = logical([zeros(6,1) ones(6,1)]);      % right

% boundary values (Dirichlet only)
btop = ones(1,nx);                  % top
bbottom = zeros(1,nx);              % bottom
bleft = ones(ny,1);                 % left
bright = zeros(ny,1);               % right

q = zeros(nx,ny);                   % right hand side (source term)
D = 1.e-5;                          % diffusivity
lambda = 2.e-7;                     % decay constant

N = nx*ny;
d = [-nx,-1,0,1,nx];
B = [ones(N,2) -(4+lambda/D)*ones(N,1) ones(N,2)];
q = reshape(q,N,1);
b = -h*h*q.*ones(N,1)/D;
for i = 1:nx
if ltop(i)
b(i) = b(i)-btop(i);
else
B(i,3) = -3;
%B(i-1,1) = 0;
end
if lbottom(i)
b(N-nx+i) = b(N-nx+i)-bbottom(i);
else
B(N-nx+i,3) = -3;
% B(N-nx+i,5) = 0;
end
end
for i = 1:ny
B(i*nx,2) = 0;
if i<ny B(i*nx+1,4) = 0; end
if lleft(i)
b((i-1)*nx+1) = b((i-1)*nx+1)-bleft(i);
else
B((i-1)*nx+1,3) = B((i-1)*nx+1,3)+1;
end
if lright(i)
b(i*nx) = b(i*nx)-bright(i);
else
B(i*nx,3) = B(i*nx,3)+1;
end
end

A = spdiags(B,d,N,N);

% processing: solution
U = A\b;
U = reshape(U,nx,ny);

% check & visualize
4*del2(U)
surf (U)

```