## Operator Splitting Monodomain Model

In Section 2.2.3 we derived a simplified version called the monodomain model, given by (2.35) and (2.36). If we combine this model with one of the models derived in Section 2.4 for the ionic current term, we have a complete model given by

xGH,

V • (MiVv) = dV + Iion(v, s) X G H, n • (MiVv) = 0

Compared with the models discussed in Section 3.1, both the monodomain model and the bidomain model are significantly more complex. They are time-dependent systems, which must be discretized in time as well as in space. They also involve coupled systems of nonlinear ODEs and PDEs, which makes it much more challenging to derive efficient numerical methods.

In this section we will present numerical methods for solving both the monodomain model and the bidomain model. The methods are based on a technique known as operator splitting, in combination with a finite difference discretization in time. The spatial discretization of the equations will be based on a finite element procedure, similar to the one we used for the simpler equations considered above.

### 3.2.1 Operator Splitting

Operator splitting is an attractive technique for solving coupled systems of PDEs, since complex equation systems may be split into smaller parts that are easier to solve. Several operator splitting techniques exist, but we will apply a class of methods often referred to as fractional step methods; see, e.g., [82]. To introduce this technique, consider an initial value problem of the form

where L1 and L2 are operators on u, and u0 is a given initial condition. If we choose a small time step At, an approximate solution at t = At may be computed by first solving the problem

for t e [0, At]. Thereafter, we solve the problem dw . .

for t e [0, At]. Note that the final value of the solution of (3.41), v(At), is used as an initial condition in this step.

Although it may now seem that we have found an approximate solution after a time interval 2At, we have only included parts of the right hand side in each integration step. To see that the result w(At) is in fact a consistent approximation to u(At), we perform a Taylor series expansion of both the original solution u, and the approximation w obtained by the operator splitting. We have

— = (L1 + L2)u, and if L1 and L2 do not depend explicitly on t, we obtain by direct differentiation ddtuU = (Li + L2 )(Li + L2)u, for which we introduce the shorter notation d2 u dt2

Repeating these steps n times gives the general result

^ = <L' + L2)'u, where the notation (L1 + L2)n simply means that the operator (L1 + L2) is applied n times to u.

Inserted into the Taylor series, this gives

u(At) = uo + At(Li + L2)uo + —(Li + L2)2uo + O(At3). (3.45)

A similar Taylor expansion can be made for the solution v of the simplified equation (3.41). We get v(At) = uo + AtL1uo + At2L21uo + O(At3).

We now use the same series expansion for the solution of (3.43), with v(At) as the initial condition. We get w(At) = v(At) + AtL2v(At) + At2L2v(At) + O(At3 ), and inserting the series expansion for v(At) gives

w(At) = uo + At(Li + L2)uq + —(A + 2L2L1 + L2,)u0 + O(At3). (3.46)

The splitting error at t = At is the difference between the operator splitting solution w(At) and the solution u(At) of the original problem. Inserting the series expansions (3.45) and (3.46), we get

We see that the error after one time step is proportional to At2, and we expect this error to accumulate to nAt2 after n time steps. When solving the equations over a fixed time interval, e.g. t e [0, b], the number of time steps n is proportional to At-1, and the error at t = b is therefore proportional to At. The splitting method, commonly referred to as Godunov splitting, is hence a first-order method.

By a small modification it is possible to make the splitting algorithm second-order accurate. The idea is that instead of first solving (3.41) for a full time step of length At, we solve the problem for a time step of length At/2. We then solve the problem (3.43) for a full step of length At, and finally (3.41) once more, again for a time interval of length At/2. The details of the three steps of the algorithm are as follows. We first solve the problem

is solved for t e [0, At]. Finally, we solve dv = Li(v), (3.51)

for t e [At/2, At]. This three-step algorithm is called Strang splitting.

To see that the Strang splitting gives second-order accuracy, we compute the Taylor series expansion of the approximate solution v at t = At. Following the steps outlined above for the Godunov splitting, we first find a Taylor expansion of the solution v of (3.47)-(3.48), at t = At/2;

At At2

Using this as the initial condition for a Taylor expansion of the solution w(At) from the second step, we get

At At2 At2 At2

w(At) = uo+—- Liuo+AtL2UoH—— L\u0+—— L2Liuo+—— L2>uo+O(At3). 2 8 2 2

And finally, by a Taylor expansion of the third step, we find

v(At) = uo + At(Li + L2)uo + — (Li + L1L2 + L2L1 + L2)uo + O(At3).

Comparing this with the Taylor expansion (3.45) of the solution of (3.39), we see that the second-order terms in the local splitting error cancel, and we have v(At) - u(At) = O(At3).

With the local error proportional to At3, the accumulated error after n ~ At-1 steps is hence proportional to At2.

We can derive a more compact notation for the operator splitting methods by introducing the concept of solution operators. We introduce operators V and W, defined so that for a given time step At, V(At) applied to uo gives the solution of (3.41) at t = At, with uo as initial condition. Similarly, W(At)uo gives the solution of (3.43) at t = At, with uo as initial condition. With this notation, the solution un obtained by Godunov splitting at t = tn = nAt is given by un = (W (At)V (At))nuo.

As above, this is just a short notation for applying the combined operator (W(At) V(At)) to uo n times. The solution obtained with Strang splitting may be written as un = (V (At/2)W (At)V (At/2))nuo.

By introducing a parameter 6 e [0,1], we can formulate a more general operatorsplitting algorithm, for which Godunov and Strang splitting are obtained as special cases. One step of the general formulation may be written as u(At) = V((1 - 6)At)W(At)V(6At)uo.

Noting that V(0) is simply the identity operator, we see that choosing 6 = 0 or 6 =1 gives two different versions of the Godunov splitting. The Strang splitting is obtained by choosing 6 = 1/2. By a Taylor series expansion similar to that above, it can be shown that all choices of 6 = 1/2 give a first-order splitting scheme. We are mostly interested in choosing 6 so that we get either Strang splitting or Godunov splitting, but other choices are possible. The greatest advantage of the general 6-formulation of the operator splitting is that it leads to a very convenient implementation of the algorithm, where the splitting method used can be varied by adjusting a single parameter.

### 3.2.2 Operator Splitting for the Monodomain Model

So far, we have introduced the ideas of Godunov and Strang splitting for the general initial value problem (3.39)-(3.40). To concretize the ideas, let us now apply the algorithm to the monodomain model (3.37)-(3.38), as proposed by Qu and Garfinkel [117].

Equation (3.37) can be written in the form (3.39), with the two operators defined by

0 0