Environmental Engineering Reference
In-Depth Information
%----------------------execution-----------------------------------
solinit = bvpinit (linspace(0,L,N),[cin; 0]);
sol = bvp4c(@bw,@bwbc,solinit,odeset,D,v,k,KS,f,cin);
%---------------------- graphical output --------------------------
plotyy (sol.x,sol.y(1,:),sol.x,sol.y(2,:));
legend ('Corg','SO_4'); grid;
%----------------------functions-----------------------------------
function
dydx = bw(x,y,D,v,k,KS,f,cin)
monod = k*y(2)/(KS+y(2));
dydx = zeros (3,1);
dydx(1) = -monod*y(1)/v;
dydx(2) = y(3);
dydx(3) = (v*y(3)+f*monod*y(1))/D;
function
res = bwbc (ya,yb,D,v,k,KS,f,cin)
res = [ya(1:2)-cin; yb(3)];
The complete code is included in the accompanying software under the name
'boudreau_westrich.m'.
The initial commands, in which the parameters
L, v, D, k, KS, f, cin
and
N
are
specified, are omitted in the listing above. Before the call of the solution routine in
the second line, an initial guess of the unknown functions is required. This is done
by the
bvpinit
command. The first formal parameter in the command is the vector
of
x
-values, for which the variable-values are to be computed. The second formal
parameter consists of three values (
cin
is a two-component column vector). Each of
these three is a guess of a constant valued function.
In the
bvp4c
call there is a list of formal parameters. First in the list are the
function calls: (1) the function in which the system of differential equations is
specified, (2) the function in which the boundary conditions are specified.
solinit
is the initial guess obtained from the previous command.
odeset
yields the standard
options for the solution routine. What follows as formal parameters is the set of
parameters for the described example.
In the
bw
function the differential equations are specified, following (
9.19
) and
(
9.20
).
y
denotes the vector of unknown variables. Thus
y(1)
denotes the concen-
tration of organic matter,
y(2)
the sulphate concentration and
y(3)
the derivative of
sulphate concentration. The Monod-term for sulphate is calculated under
monod
.
The
bwbc
function specifies the boundary conditions. At the left side, i.e. at the
lowest
x
-value (
a
) given above, the condition is computed using the
ya
variable.
The variable
yb
is responsible for the condition at the highest
x
-value (
b
) . The
condition is formulated in a way making the residual
res
vanish. Thus the first
two variables at boundary
x ¼ a
take the values given in the
cin
vector, while the
third variable has a vanishing value at boundary
x ¼ b
. The physical meaning of
the latter condition is that there is no diffusive flux because of a vanishing
concentration gradient.
