CardioWave | Dev / ImplicitTimeIntegrator browse

We will use Newton's Method to formulate a fully implicit time-integrator. We first recast the equations as:

We then compute the Jacobian of the system, that is, the set of partial derivative of all the functions in G() with respect to all of the variables in x. This is then a matrix of partial derivatives. We proceed with the Newton algorithm by solving the linear system:

where delta is a correction term. The superscripts in k indicate the Newton iteration count. After computing delta, we update the vector x:

where we choose the scalar s to force the norm of G() to decrease. If the norm of G() is below some threshold, we consider this time-step to be solved. Otherwise, we must recompute the Jacobian (since it is a function of x, which was just updated), solve the new linear system for delta, find s and update x, etc. This would constitute the second Newton iteration.

For the monodomain equations, we start with:

The Jacobian (linear) system looks like:

where the remaining partial derivative terms come from the membrane model being used. We define a matrix Q as the upper sub-block:

Note that due to the structure of the membrane functions, Q will be, at worst, block diagonal with a block size equal to the membrane model's patch size (ie. 3x3 blocks for the HH model with 3 state variables). This also means that inverting Q will not be extremely difficult, even on parallel machines. Following Nigel Hooke's work, we then use one step of block Gaussian elimination to reduce the linear system to:

where matrix T and vector z are defined as:

Note that due to the sizes of dFv/dq, Q, and dFq/dv, the T matrix has terms added only to it's main diagonal. To solve the linear system, then, we start by solving the smaller system in T:

Once we've computed delta-v, we can back solve to find delta-q:

For the bidomain equations, we again start with the G(x)=0 formulation:

The Jacobian (linear) system is now a little larger, in block form it is:

We define the matrix Q as above. Rewriting the Jacobian with the lower 2x2 sub-block as WXYZ:

We again do one step of block Gaussian elimination to reduce the size of the linear system. In the bidomain case, we must modify the WXYZ portion of the system:

where, similar to what we saw before, we define the matrix T and vectors zi and ze as:

In this case, T is strictly diagonal (however, it is added to off-diagonal blocks of the matrix). Note that those terms added to the main diagonal in the monodomain all show up here in the bidomain, we have simply computed them slightly differently. Once we have solved the linear system in delta-v-i and delta-v-e, we then backsolve to find delta-q:

Note that in both the implicit mono- and bi-domains, we will require the membrane model programmer to compute:

(where the scaling by beta and A can be handled by the time-stepper programmer), and also to do the backsolve, that is, to solve:

for given vectors x and y. Since these functions will require detailed knowledge of the membrane dynamics, the membrane module programmer is the only one in a position to do this (ie. this is not general-purpose functionality, it requires specific knowledge of the membrane module). Also, the structure of the membrane functions, which ionic flows interact with what state variables, affects the structure of the Q. In particular, many membrane modules may be re-ordered such that the Q matrix is lower triangular (Hooke refers to these as having "Property H") or even diagonal. Take the HH model, for example. There are no interactions between the state variables:

The rest of the Jacobian will be zeros:

Thus, there may be extremely efficient means of inverting Q, but only the membrane model programmer will be in a position to know that.