KX=F
Overview
Lattice equations showed how forcebalance can be applied to determine the position of each node. This process yields a system of equations that can be written as
\[K X = F\]where:
 \(K\) is a square matrix of size \(n\) x \(n\)
 \(X\) is a vector of size \(n\) x 1
 \(F\) is a vector of size \(n\) x 1
\(K\) and \(F\) can be calculated as described at lattice equations. The goal is to calculate \(X\). This process can be written as:
 \(X\) = \(K\) \ \(F\)
Size of the problem
The main challenge is the size of the problem. Consider the smallest possible lattice which has just 4 thick filaments.
 \(n\) = (no of thick filaments * no of nodes per thick filament) + (no of thin filaments * no of nodes per thin filament)
 \(n\) = (4 * 54) + (8 * 189) = 1728
Doing anything 1728 times will require some resources. In practice, simulations often need 100 or more thick filaments to produce smooth records at low levels of activation. With 100 filaments, \(n\) has increased to 43200.
Standard approaches to solve \(X = K \backslash F\) scale as \(n^3\) implying that calculating the node positions for 100 thick filaments will require 73 million operations. And that is for each timestep! The computational cost becomes prohibitive.
Strategy
As explained at lattice equations, the \(K\) matrix is sparse, meaning that most of the entries are 0. FiberSim takes advantage of this and uses special techniques that are designed for such matrices. In fact, the main codebase never stores the full \(K\) matrix in socalled dense form.
Principles

If a sparse matrix \(A\) is stored in Compressed Sparse Column (CSC) form, and \(B\) is a vector, the operation \(A.B\) can be performed efficiently.

In the special case where a matrix \(T\) is tridiagonal (values only on the leading diagonal, and diagonals immediately above and immediately below), there is an efficient algorithm to solve \(X = T \backslash F\)

As explained at lattice equations, the \(K\) matrix for FiberSim is dominated by the tridiagonals.
Approach
The goal is to solve \(X\) where \(K\) and \(F\) are known and
\[\begin{equation} K X = F \end{equation}\]Define \(K_0\) as the portion of \(K\) that correspond to the thin and thick filament backbones, noting that \(K_0\)’s nonzero elements all lie on the tridiagonals.
Define \(G(X)\) as a vector that depends on X and represents perturbations induced by crosslinks (in FiberSim’s case, titin, bound myosins, and bound cMyBPC molecules).
\[\begin{equation} K_0 X + G(X) = F \end{equation}\]Further define
\[\begin{equation} X_i = X_{g,i} + \Delta X_i \end{equation}\]where
 \(X_i\) is the ‘true’ \(X\)
 \(X_{g,i}\) is the i’th ‘guess’ at \(X_i\)
 \(\Delta X_i\) is the difference between them.
Substituting in, leads to
\[\begin{equation} K_0 (X_{g,i} + \Delta X_i) + G(X_{g,i} + \Delta X_i) = F \end{equation}\]This can be rearranged to yield
\[\begin{equation} \Delta X_i = K_0 ^{1} (F  K_0 X_{g,i}  G(X_{g,i} + \Delta X_i)) \end{equation}\]If \(G\) is dominated by \(X_{g}\), equation (5) can be rewritten as
\[\begin{equation} \Delta X_i \approx K_0 ^{1} (F  K_0 X_{g,i}  G(X_{g,i})) \end{equation}\]An improved guess for \(X\) can then be calculated from
\[\begin{equation} X_{g,i+1} = X_{g,i} + \Delta X_i \end{equation}\]Note that

\(K_0 X_{g,i}\) can be evaluated efficently if \(K_0\) is stored in compressed sparse column format

\((F  K_0 X_{g,i}  G(X_{g,i}))\) simplifies to \(some \ vector\)

\(\Delta X = K_0 ^{1} (some\ vector)\) can be calculated efficiently because \(K_0\) is tridiagonal
This leads to an iterative approach that runs as follows:
 Start with an initial guess for \(X_{g,i}\)
 FiberSim uses the tridiagonal solution for the first timestep and \(X\) from the last timestep thereafter.
 In a loop
 calculate \(G(X_{g,i})\)
 calculate \(\Delta X_i\)
 form a new \(X_{g,i+1}\) = \(X_{g,i} + \alpha \Delta X_i\)
 repeat until \(X_g\) converges
Since FiberSim Ver 2.3, \(\alpha\) is initialized at 1 but can be reduced to improve convergence. Specifically, if \(X\) starts to diverge, the current value of \(\alpha\) is reduced by 50% for the next step. This adds a smaller correction.
This approach
 maintains fast convergence for situations where the filaments are relatively stiff
 improves the likelihood of finding the correct \(X\) vector when the filaments are so compliant that the perturbations induced by crosslinks are comparable to the internode spacing.
An aside that didn’t work
FiberCpp uses the GSL Scientific Library for most of the scientific computing procedures.
GSL includes code that can solve \(X = K \backslash X\) when \(K\) is stored in a sparse format. Specifically, the algorithm implements a Generalized Minimum Residual (GMRES) method.
Pilot tests showed that the GMRES algorithm could calculate \(X\) under at least some conditions but the approach was so slow that the calculations became impractical.
Alternative methods may be more efficent.