Environmental Engineering Reference
In-Depth Information
At all nodes at which we require a Neumann condition, the 4 entry in the main
diagonal needs to be changed, according to formula (
21.18
). If we assume a
Neumann condition on the left and right boundary, we thus have to change the
corresponding diagonal entries to 3. This is done in the last two statements of the
first loop.
for
i = 1:ny
B(i*nx,2) = 0;
B(i*nx+1,4) = 0;
B((i-1)*nx+1,3) = -3;
B(i*nx,3) = -3;
end
At all nodes with Dirichlet condition, according to (
21.17
) we have to add an
additional term on the right hand side, i.e. in the vector
b
. This is done in the next
loop, before the matrix
A
finally is set up as intended.
for
i = 1:nx
b(i) = b(i)-btop;
b(N+1-i) = b(N+1-i)-bbottom;
end
A = spdiags(B,d,N,N);
After all these preparations we can solve the linear system. The solution is
obtained in the vector
U
using the backslash operator (see chapter *). In order to
rearrange calculated (before) unknown values into a matrix that corresponds to the
format of the grid, we use the
reshape
command.
% processing: solution
U = A\b;
U = reshape(U,nx,ny);
command
del2
computes the
outcome of the right hand side of (
21.23
). In our example this is a constant value
at all nodes. The last command delivers graphical output, as shown in Fig.
21.3
.
Finally we check the result. The MATLAB
®
% check & visualize
4*del2(U)
surf(U)
The complete code is included in the accompanying software under the name
'Poisson1.m'
Note that the Dirichlet boundary values, in the example 1 on the left and 0 on the
right side, are not taken at the outer nodes. In the given example all nodes are lying
in the interior of the model region. The boundaries, for which the conditions are
valid, have a distance of
h
from the outer nodes, as the m-file terms for the boundary
blocks correspond with the discrete formulation (
21.17
). In order to show the
solution in the entire model region, a vectors with boundary values would have to
be combined with
U
, for example by using:
surf([btop*ones(nx,1) U bbottom*ones(nx,1)])
