Environmental Engineering Reference
In-Depth Information
nx = 5; ny = 20;
% dimensions in x- and y-direction
h = 1;
% grid spacing
btop = 1;
% boundary condition at left side
bbottom = 0;
% boundary condition at right side
q = .1;
% right hand side (source term)
In the main part we compute the matrix and the right hand side, as shown in
(
21.20
) according to the finite difference method.
N
is the total number of
unknowns. In the example we compute an approximate solution at all locations of
a grid of 5
20. For the exact set-up of the matrix (
21.20
) we use the
spdiags
command (see Sidebar 21.1) before starting the 'real' processing (the solution of
the linear system). All commands before that statement prepare the exact entries of
the diagonals.
The vector
d
determines the locations of the diagonal. To understand this we
have to be aware that the locations of the grid, the nodes, are numbered in canonical
order, i.e. starting with the 1 for the (1,1) position we number in
x
-direction first and
in
y
-direction subsequently. In that order their function values appear in the
unknown vector x. The left and right (
x
-direction!) neighbors of the node
i
then
have the indices
i
+ 1 and
i
-1, and the top and bottom (y-direction) neighbor nodes
have the indices
i+n
x
and
i-n
x
. This determines the locations of the off-diagonals,
and is taken into account by the given choice of
d
.
As mentioned above, annlogously to the derivation of formula (
21.19
) for the
Poisson equation we obtain the discretization:
4
uðx; yÞ¼h
2
q
(21.22)
uðx þ
D
x; yÞþuðx
D
x; yÞþuðx; y þ
D
yÞþuðx; y
D
yÞ
or
h
2
q
4
(21.23)
1
4
ð
uðx þ
D
x; yÞþuðx
D
x; yÞþuðx; y þ
D
yÞþuðx; y
D
yÞ
4
uðx; yÞ
Þ ¼
4 entries in the main diagonal and 1 is the standard entry in the
off-diagonal. This we have to consider in the assignment of matrix
B
, in which all
diagonals of the matrix
A
are stored. Thus we proceed with the commands:
Thus there are
N = nx*ny;
d = [-nx,-1,0,1,nx];
B = [ones(N,2) -4*ones(N,1) ones(N,2)];
b = -q*h*h*ones(N,1);
However, there are some exceptions that must be taken into account. For some
nodes not all four neighbor nodes are available. For example along the left bound-
ary the left neighbor is missing, and at all nodes on the right side the right neighbor
is missing. Thus in the matrix the corresponding off-diagonal entries should be
0 and not 1. In the example this is done in the first two commands of the first loop.
