Strongly Implicit Procedure (SIP)

We want to solve the following linear system of equations

  $$A \, \vec{N} = \vec{b}$$ (6.1)

where $A$ is some non-symmetric penta-diagonal matrix, $\vec{N}$ is the wave action vector to be solved and $\vec{b}$ contains source terms and boundary values.


The basis for the SIP method (Stone, 1968; Ferziger and Perić, 1999) lies in the observation that an LU decomposition is an excellent general purpose solver, which unfortunately can not take advantage of the sparseness of a matrix. Secondly, in an iterative method, if the matrix $M = LU$ is a good approximation to the matrix $A$, rapid convergence results. These observations lead to the idea of using an approximate LU factorization of $A$ as the iteration matrix $M$, i.e.:

  $$M = L\,U = A + K$$ (6.2)

where $L$ and $U$ are both sparse and $K$ is small. For non-symmetric matrices the incomplete LU (ILU) factorisation gives such an decomposition but unfortunately converges rather slowly. In the ILU method one proceeds as in a standard LU decomposition. However, for every element of the original matrix $A$ that is zero the corresponding elements in $L$ or $U$ is set to zero. This means that the product of $LU$ will contain more nonzero diagonals than the original matrix $A$. Therefore the matrix $K$ must contain these extra diagonals as well if Eq. (6.2) is to hold.


Stone reasoned that if the equations approximate an elliptic partial differential equation the solution can be expected to be smooth. This means that the unknowns corresponding to the extra diagonals can be approximated by interpolation of the surrounding points. By allowing $K$ to have more non zero entries on all seven diagonals and using the interpolation mentioned above the SIP method constructs an LU factorization with the property that for a given approximate solution $\phi$ the product $K\phi \approx 0$ and thus the iteration matrix $M$ is close to $A$ by relation (6.2).


To solve the system of equations the following iterations is performed, starting with an initial guess for the wave action vector ${\vec{N}}^0$ an iteration is performed solving:

  $$U\,{\vec{N}}^{s+1} = L^{-1}\,K\,{\vec{N}}^s + L^{-1}\,\vec{b}$$ (6.3)

Since the matrix $U$ is upper triangular this equation is efficiently solved by back substitution. An essential property which makes the method feasible is that the matrix $L$ is easily invertible. This iterative process is repeated $s=0,1,2,\cdots$ until convergence is reached.