Efficient Eigenvector Calculation method

The carbon and climate models both use about 40 ocean layers (x2 oceans for climate model). If using a simple integration method, the timestep would have to be quite small (e.g. 1 month) in order that the fluxes are much smaller than the box contents, and so the calculation over hundreds of years is rather slow. Therefore in order to get an instant response in the plots as the parameters are dragged by the mouse, a more efficient matrix-eigenvector method has been applied.

For examples see:

  • Carbon Module
  • Climate Module

    Note also:

  • How does JCM work so fast?:

    Features of this method

    :
  • An exact analytical solution
  • The only approximation is that the change in non-linear fluxes is linear within one timestep.
  • The timestep can be any size (depending on input data and plot resolution)
  • The ramp function must be iterated when a non-linear flux is a function of box contents.
  • Only matrix cells needed for data output or non-linear input are calculated.
  • Usually more efficient than pulse response function especially with many timesteps.
  • Unlike PRF, preserves knowledge of all box contents (needed for thermal expansion or depth profiles etc.)

    The eigenvectors only need to be recalculated when you change a parameter that alters the linear fluxes (such as internal mixing rates), not when the external emissions change.

    Eigenvectors, inverses, etc. are calculated using the convenient Java Matrix Package ("JAMA"), which was developed at MIT, and has simply been bolted onto Java Climate Model.
  • See: http://math.nist.gov/javanumerics/jama/

    Mathematical Principle

    Differential equation

    Start with a simple differential equation:

    dq/dt = -λq + x

    If x changes linearly from xt to xt+dt, it can be shown that:

    qt+dt = prop qt + step xt + ramp (xt+dt - xt)

    Where:
    prop = e-λdt
    step = (1 - prop)/λ
    ramp = (dt - step)/(λ.dt)

    System of boxes (n equations)

    A system of n boxes with diffusive fluxes between them, plus extra non-linear inputs, can be written as:

    dQ/dt = R Q + X

    Where Q is the vector of box contents, R is an nxn matrix for the fluxes, and X is the vector of extra inputs.

    We can diagonalise R:

    R = S diag(E) S-1

    Where S is the matrix of eigenvectors, and diag(E) the diagonal matrix of eigenvalues.

    Premultiplying all by S-1 gives:

    d/dt(S-1Q) = diag(E) (S-1Q) + S-1X

    S-1Q and S-1X are just vectors, with each element of S-1Q multiplied by just one eigenvalue, so this is just like n simple differential equations as above, with diag(E) providing the λ terms

    Now, if we have a lot of equal timesteps, we only need to calculate the n prop, step and ramp functions once, at the beginning.

    If we preserve the state information in the form S-1Q, then the diffusion for each timestep is given just by multiplying each element by it's prop function.

    To add the step and ramp functions for extra inputs, we need to calculate S-1X, but if we have only y boxes with non-linear inputs (for example due to carbonate chemistry at the surface of an upwelling-diffusion ocean), we only have to multiply out y columns of the nxn matrix S-1 (which we can also premultiply by the step and ramp functions just once at the beginning).

    To plot a result, we also need to convert S-1Q back to Q by premultiplying with S, but if we are plotting only z elements of Q (for example, the contents of the surface layer), we only have to multiply z rows of the nxn matrix S.

    So if there are t timesteps, we have altogether
    t*n*(1+2*y+z) multiplications
    (the 2 is for step + ramp functions)

    (Note that if (1+2*y+z)>n, it might be quicker preserving the state in Q rather than S-1Q).

    Iteration

    If the X are not just external inputs, but are also dependent on the contents of the box (as with the carbonate chemistry), then we have to iterate. A good first guess is usually that the change in X is the same as in the previous timestep. This can still be fast and accurate, providing the assumption that X changes linearly within one timestep is reasonable.

    If X is a non-linear function of box contents, it helps to separate out an arbitrary linear part of X and include this within R, so that the remaining non-linear perturbation is as small as possible.

    Compare with PRF

    Note that (im)pulse response functions are derived from the same principle: an exact PRF for a linear system is derived from the n prop functions as above. The difference lies the way of adding it up:

    With PRFs, the fate of each input (Xt) from every timestep must be calculated seperately for all subsequent timesteps. So if there are n boxes and t timesteps, then we have altogether n * t2 /2 multiplications.

    So this method seems to be slower, so long as t > 8!
    (note for comparison y and z must be 1, since PRF can only cope with input and output in one box).