next up previous contents
Next: Metabolism Up: The Standard Model Previous: Mass Conservation and Nonnegativity

Numerics

Note that all entries in the system matrix of (3.3) except those on the diagonal are nonnegative, while the diagonal entries are nonpositive. Moreover, the inhomogeneous part is nonnegative (or, more precisely, strictly positive for compartment tex2html_wrap_inline2749) if and only if the air compartment contains a positive pollutant mass and there is an exchange between the air compartment and compartment k. Consider for the moment the case that the inhomogeneous part is zero, i. e. tex2html_wrap_inline2935 or tex2html_wrap_inline2937, and that no backflow to the air compartment takes place, i. e. tex2html_wrap_inline2939. As a consequence, all column sums of the matrix A (t) are zero for all t, and any consistent multistep or Runge-Kutta method for solving (3.3) will lead to mass consistency, regardless what step size is used and no matter if the method is explicit or implicit.

Due to (3.1) and the theorem of Gershgorin, the eigenvalues of A (t) have only nonpositive real parts (for all tex2html_wrap_inline2787), and the only possible eigenvalue with real part equal to zero is 0. This makes the system stiff, while a lower bound for the real part of the eigenvalues at time t is
displaymath3787
At first glance an implicit numerical solver should be used to solve (3.3). But implicit methods require the solution of a linear system in n-1 variables per time step. In contrast to this, an explicit solver just needs one matrix-vector multiplication to compute the next approximation. Moreover, the accuracy demand in the application discussed here is rather moderate, since the inhomogeneous part of (3.3) cannot be calculated with high accuracy. As a consequence, we use Euler's method in the present implementation of the compartment model within OLAF. Nonnegativity of the iterates is achieved by the following step size control. In the homogeneous case, we have
displaymath3788
as the iteration procedure. Here, the superscripts denote the time steps for the approximations tex2html_wrap_inline2955 that we compute with step length h > 0. Given a (componentwise) nonnegative vector tex2html_wrap_inline2959, the nonnegativity of tex2html_wrap_inline2961 can be assured if tex2html_wrap_inline2963 for tex2html_wrap_inline2785. The step size chosen by OLAF is then
 equation309
the largest one possible under this constraint. Note that larger step sizes can also be allowed, e. g. if
displaymath3789
for a tex2html_wrap_inline2967, tex2html_wrap_inline2969, then use a step size h with
 equation323
This rule depends then on the values of tex2html_wrap_inline2959. Even larger steps can be made if the inhomogeneous term is used in the above estimate:
 equation328
for all tex2html_wrap_inline2785 with tex2html_wrap_inline2977.

With the implicit Euler method,
displaymath3790
i. e.
 equation336
one may use arbitrary step sizes h > 0 without violating the nonnegativity of tex2html_wrap_inline2961. This can be seen in the following way. The unique solution tex2html_wrap_inline2983 of the homogeneous version of (3.3),
 equation345
is nonnegative. Now write A(t) = N(t) - D(t) with the positive diagonal matrix tex2html_wrap_inline2987 and define the operator
displaymath3791
for arbitrary vectors tex2html_wrap_inline2989 and tex2html_wrap_inline2787. Note that tex2html_wrap_inline2993 for all y and all tex2html_wrap_inline2787. With this, the system of differential equations
 equation358
has a unique solution tex2html_wrap_inline2879. But then tex2html_wrap_inline3001, and consequentially tex2html_wrap_inline3003 for all tex2html_wrap_inline2787. Therefore, employing the implicit Euler method on (3.10) leads to the same approximations tex2html_wrap_inline3007 as if the method has been used to solve (3.9). But then
displaymath3792
and when the right hand side is nonnegative (by induction), tex2html_wrap_inline3009 will also be nonnegative. Note that this reasoning works only if the system of linear equations (3.8) is solved exactly, i. e. with a direct method and not an iterative one. Therefore, a nonnegativity enforcing version of Euler's implicit method needs tex2html_wrap_inline3011 floating point operations per time step, as compared to the tex2html_wrap_inline3013 floating point operations for the explicit version. More precise estimates depend on the sparsity structure of the system matrix A.

Note, however, that computing an approximative solution of this system of ordinary differential equations is necessary for each point of the computational grid of the air dispersion code, and that any step size rule depending on tex2html_wrap_inline2959 decouples the time grid of the different systems to be solved. This poses a problem when the actual objective function has to be evaluated, because each mass in each compartment is then also equipped with a time interval representing the length of stay of the mass inside the compartment. Due to this, further tests with step size rules like (3.6) or (3.7) have been postponed.

While atmospheric dispersion models usually use a rather large time scale, chemokinetic models do not. Typical time steps in atmospheric dispersion are one hour or several hours, while the scale needed in the chemokinetic model is, according to what is said above, only minutes. As a consequence, the concentration of the pollutant in the air as used in the ecosystem dispersion and the chemokinetic model is approximated by a piecewise constant function. While it is possible to compute out of this data a smooth approximation, it is not clear how mass balance can be achieved. On the other hand, nonsmooth right-hand sides in (3.3) pose serious numerical problems to the corresponding solver. Due to this, Euler's method is used only from one time step of the atmospheric dispersion model to the next. The last step size used in Euler's method is computed in such a way that the step taken coincides with the end point of the interval in which the current constant approximation of pollutant concentration in the air holds. The next step of Euler's method uses the new pollutant concentration from the next time step of the atmospheric dispersion model.


next up previous contents
Next: Metabolism Up: The Standard Model Previous: Mass Conservation and Nonnegativity

Joerg Fliege
Wed Dec 22 12:25:31 CET 1999