The origin of the system

The system we are facing, is used in the process of solving a set of partial differential equations, describing how pollutants react with each other in the atmosphere.

It should be of no surprise to the reader, that chemical reactions between a number of species can be written in differential form. Let's look at a little example of this. Consider a system consisting of four species, each with concentrations $c_1$, $c_2$, $c_3$ and $c_4$. 1 and 2 react with each other, producing 3, while 4 decays producing 1. The concentration of species 1, $c_1$ is thus given as:

$$\dot{c}_1 = - c_1 c_2 + \alpha c_4$$ (1.1)

where $\alpha$ indicates how fast $c_4$ is decaying. This of course, is a very simple model of a reaction.

When modelling aerial pollution, we work with a fairly large number of species. This number $N$ depends on the model in use, but is often between 36 and 56, but can go as high as more than 100 species for the most complex models.

Working with many species, invites us to write this reaction scheme on a more compact form. For every species in our model, we are given the equation:

$$\dot{c}_i = \sum_{j,k} \beta_{ijk} c_j c_k + \sum_j \alpha_{ij} c_j.$$ (1.2)

In real applications however, the coefficients $\alpha$ and $\beta$ are often zero, since most chemicals only react (notably) with a small subset of the others.

The implicit Euler method is used for solving the resulting system of differential equations. The system of equations is written short as $\dot{c} = f(t,c)$. The Euler method lets us take small steps $h$ in time, thereby finding the concentrations of the species at any given time. A step in time, from some given initial condition $c_0$ is given by

$$c = c_0 + hf(t,c).$$ (1.3)

Or, put in another way, we find the concentrations after a time-step by finding a root in the equation

$$g(c) = c_0 - c - hf(t,c).$$ (1.4)

Simple Newton iteration yields

$$c^{i+1} = c^i - \left( \frac{\partial g}{\partial c}(c^i) \right)^{-1} g(c^i).$$ (1.5)

The Jacobi matrix $\mathbf{J}$ for $f$ is given by $\partial f / \partial c$. Thus,

$$\frac{\partial g}{\partial c} = -\mathbf{I} + h \frac{\partial f}{\partial c} = -\mathbf{I} + h\mathbf{J}.$$ (1.6)

Finding the concentrations of the species seems to be the ``simple'' matter of iterating thru equation 1.5. That is in a sense true, except from the fact, that iterating thru that equation involves inverting $-\mathbf{I}+h\mathbf{J}$. This is the reason we actually bother to discuss the efficient solution of a linear system of equations in the first place. It should also be clear to the reader, why it is so important that the solution of that system of equations run efficiently.

The reaction schemes doesn't change often. They are usually developed and implemented, and then used for years, for research and forecasts etc. Therefore, if we can only solve the model in an efficient way, it is perfectly reasonable to spend much time optimizing the solver, since an improvement there will benefit the users of the model for years.