Environmental Engineering Reference
In-Depth Information
% boundary values (Dirichlet only)
btop = ones(1,nx); % top
bbottom = zeros(1,nx); % bottom
bleft = 2*ones(ny,1); % left
bright = zeros(ny,1); % right
The values given in the four vectors above will be taken only if the
corresponding position has a Dirichlet boundary condition, i.e. when the
corresponding entry in the logical vector is zero. In the example here only the first
half of the vector bleft and the second half of the bright are relevant for the
problem. All other values are maintained in order to allow M-file to be capable to
treat the general case.
As a generalization of the right hand side we write
q = [ones(nx/2,ny);zeros(nx/2,ny)]; % right hand side (source term)
The entries in the matrix q represent positions in the model region. If an element
is non-zero there is a source (or sink if the value is negative) in the corresponding
block of the model region. In the example we have a source term in the upper half of
the model region, while there is no source in the lower part - for whatever reason,
we don't have to discuss here. The next commands are identical to the former
program.
N = nx*ny;
d = [-nx,-1,0,1,nx];
B = [ones(N,2) -4*ones(N,1) ones(N,2)];
b = -h*h*reshape(q,N,1);
The last command only is modified. One has to take into account that the matrix
q has to be changed into a vector for the further computations. In the following we
modify matrix B and right hand side b according to the general boundary conditions.
In the first loop we treat top and bottom boundary conditions. For the top boundary
(indices 1 to n x ) an additional term has to be added to the right hand side in case of a
Dirichlet condition or the entry in the main diagonal has to be reset to 3 in case of
Neumann condition. These settings were outlined already in the previous subchap-
ter, but for the different boundaries as a whole. Now the distinction has to be made
for each location along the boundaries. The bottom boundary blocks then (indices
N - n x +1to N ) are treated similarly in the second command block.
for i = 1:nx
if ltop(i)
b(i) = b(i)-btop(i);
else
B(i,3) = -3;
end
Search WWH ::




Custom Search