dcmri.conc_ncomp#

dcmri.conc_ncomp(J, T, E, t=None, dt=1.0, solver='diag', dt_prop=None)[source]#

Concentration in a linear and stationary n-compartment system.

See section N-compartment system for more detail.

Parameters:

  • J (array_like) – the indicator flux entering the system, as a rectangular 2D array with dimensions (n,k), where n is the number of compartments and k is the number of time points in J.

  • T (array_like) – n-element array with mean transit times of each compartment.

  • E (array_like) – dimensionless and square n x n matrix. An off-diagonal element E[j,i] is the extraction fraction from compartment i to compartment j. A diagonal element E[i,i] is the extraction fraction from compartment i to the outside.

  • t (array_like, optional) – the time points of the indicator flux J, in the same units as T. If t is not provided, the time points are assumed to be uniformly spaced with spacing dt. Defaults to None.

  • dt (float, optional) – spacing between time points for uniformly spaced time points, in the same units as T. This parameter is ignored if t is explicity provided. Defaults to 1.0.

  • solver (str, optional) – A string specifying the numerical method for solving the system. Two options are available: with solver = 'diag' the system is solved by diagonalising the system matrix, with solver = 'prop' the system is solved by forward propagation. The default is 'diag'.

  • dt_prop (float, optional) – internal time resolution for the forward propagation when solver = 'prop'. This must be in the same units as T. If dt_prop is not provided, it defaults to the sampling interval, or the smallest time step needed for stable results (whichever is smaller). This argument is ignored when solver = 'diag'. Defaults to None.

Returns:

Concentration in each compartment, and at each time point, as a 2D array with dimensions (n,k), where n is the number of compartments and k is the number of time points in J.

Return type:

numpy.ndarray

See also

res_ncomp, prop_ncomp, flux_ncomp

Note

The default solver 'diag' should be most accurate and fastest, but currently does not allow for compartments that trap the tracer. It relies on matrix diagonalization which may be more problematic in very large systems, such as spatiotemporal models. The alternative solver 'prop' is simple and robust and is a suitable alternative in such cases. It is slower and less accurate, though the accuracy can be improved at the cost of larger computation times by setting a smaller dt_prop.

Example

>>> import numpy as np
>>> import dcmri as dc

Consider a measurement with 10 time points from 0 to 20s, and a 2-compartment system with a constant influx in each compartment. The influx in compartment 1 is twice a large than in compartment 0:

>>> t = np.linspace(0, 20, 10)
>>> J = np.zeros((2, t.size))
>>> J[0,:] = 1
>>> J[1,:] = 2

The transit times are 6s for compartment 0 and 12s for compartment 1.

>>> T = [6,12]

The extraction fraction from compartment 0 to compartment 1 is 0.3 and the extraction fraction from 1 to 0 is 0.8. These are the off-diagonal elements of E. No indicator is trapped or created inside the system so the extraction fractions for each compartment must add up to 1. The extraction fractions to the outside are therefore 0.7 and 0.2 for compartment 0 and 1, respectively. These are the diagonal elements of E:

>>> E = [
...  [0.7, 0.8],
...  [0.3, 0.2]]

Calculate the concentrations in both compartments of the system:

>>> C = dc.conc_ncomp(J, T, E, t)

The concentrations in compartment 0 are:

>>> C[0,:]
array([ 0.        ,  2.13668993,  4.09491578,  5.87276879,  7.47633644,
8.91605167, 10.20442515, 11.3546615 , 12.37983769, 13.29243667])

The concentrations in compartment 1 are:

>>> C[1,:]
array([ 0.        ,  4.170364  ,  7.84318653, 11.0842876 , 13.94862323,
16.48272778, 18.72645063, 20.71421877, 22.47597717, 24.03790679])

Solving by forward propagation produces a different result because of the relatively low time resolution:

>>> C = dc.conc_ncomp(J, T, E, t, solver='prop')
>>> C[1,:]
array([ 0.        ,  4.44444444,  8.3127572 , 11.69333943, 14.65551209,
17.25550803, 19.54012974, 21.54905722, 23.31636527, 24.87156916])

But the difference can be made arbitrarily small by choosing a smaller dt_prop (at the cost of some computation time). In this case the results become very close with dt_prop = 0.01:

>>> C = dc.conc_ncomp(J, T, E, t, solver='prop', dt_prop=0.01)
>>> C[1,:]
array([ 0.        ,  4.17147736,  7.84511918, 11.08681805, 13.95158088,
16.48597905, 18.72988986, 20.71776196, 22.47955758, 24.04147164])