Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a problem of the form

$$\left[\begin{array}{cc} -(\lambda+2\mu)\frac{d^2}{dx^2} & \alpha\frac{d}{dx} \\ \frac{\alpha}{\Delta t}\frac{d}{dx} & \frac{c_0}{\Delta t}I-\frac{\kappa}{\mu_f}\frac{d^2}{dx^2} \end{array}\right] \left[ \begin{array}{c} u\\ p\end{array}\right] =\left[\begin{array}{c} f\\g \end{array}\right]$$

When I apply a finite difference scheme to this system, I obtain a block matrix system of the form

$$\left[ \begin{array}{cc} m_{11} & m_{12}\\ m_{21} & m_{22}\end{array}\right]\left[ \begin{array}{c} \vec{u}\\ \vec{p}\end{array}\right] =\left[\begin{array}{c} f\\g \end{array}\right]$$

where $m_{ij}$ are all sparse symmetric matrices and $u_i$, $p_i$ are the discrete values of u & p at the discrete points $x_i$ in the domain.

I was thinking of solving this system of equations by a fixed point iteration by splitting the block matrices into the form

$$\left[ \begin{array}{cc} m_{11} & m_{12}\\ 0 & m_{22}\end{array}\right]\left[ \begin{array}{c} \vec{u}\\ \vec{p}\end{array}\right]^{k+1} +\left[ \begin{array}{cc} 0 & 0\\ m_{21} & 0\end{array}\right] \left[ \begin{array}{c} \vec{u}\\ \vec{p}\end{array}\right]^k =\left[\begin{array}{c} f\\g \end{array}\right]$$

I recall from basic numerical analysis that if $A=M+N$ and we iterate $Mx^{k+1}=-Nx^k+b$, then the method is guaranteed to converge from any starting vector $x^0$ if the eigenvalues of $M^{-1}N$ are $<1$ in absolute value. Is there an easy way to verify whether this is true for this particular matrix splitting? Can I apply an analogous proof as that of standard gauss-seidel or jacobi methods?

Note: I know that this system can be solved by other (possibly faster) means in the 1D case, but I'm particularly interested in this form because it allows me to decouple the multiphysics simulation and I can implement fast solvers for each equation in the more complicated multidimensional case.

share|improve this question

When you write down the product

$$M^{-1}N = \begin{pmatrix} m_{11}^{-1} & -m_{11}^{-1} m_{12} m_{22}^{-1} \\ 0 & m_{22}^{-1} \end{pmatrix} \begin{pmatrix} 0 & 0 \\ m_{21} & 0 \end{pmatrix} = \begin{pmatrix} -m_{11}^{-1} m_{12} m_{22}^{-1} m_{21} & 0 \\ m_{22}^{-1} m_{21} & 0 \end{pmatrix}$$

it is clear that the eigenvalues of $M^{-1} N$ depend on $\Delta t$. Note that you can consider the generalized eigenvalue problem $m_{12} m_{22}^{-1} m_{21} x = \lambda m_{11} x$. The stability of this iteration depends on $\Delta t$ and the relative sizes of the two terms in $m_{22}$.

If you want a robust solver based on splitting, you should look at approximate commutator preconditioners. The idea there is to approximate the Schur complement $m_{11} - m_{12} m_{22}^{-1} m_{21}$ or $m_{22} - m_{21} m_{11}^{-1} m_{12}$ using something that is feasible to apply. These are usually combined with a Krylov method to catch stray eigenvalues arising because the approximate commutator arguments are inexact, especially at boundaries.

share|improve this answer
I found some mistakes in my original question: the entire first row was written incorrectly. I have updated the question to reflect this. – Paul Feb 7 '13 at 20:04
Jed: I think your calculation of the product $M^{-1}N$ is incorrect, since $M^{-1}=\left[\begin{array}{cc} m_{11}^{-1} & -m_{11}^{-1}m_{12}m_{22}^{-1} \\ 0 & m_{22}^{-1}\end{array}\right]$ – Paul Feb 7 '13 at 21:42
The new version of the problem is better behaved, but the splitting is still only conditionally stable. Thanks for spotting my bug. – Jed Brown Feb 8 '13 at 0:04
So, the goal of the approximate commutator is to precondition $M^{-1}N$ such that the eigenvalues are closer to (or less than) 1 in absolute value? – Paul Feb 12 '13 at 16:04
The goal of any preconditioner for a stationary iteration is to make the eigenvalues of the iteration matrix $T = (1 - P^{-1} A)$, which, in this case, becomes $1 - M^{-1} (M+N) = -M^{-1}N$ as small as possible. For a Krylov method, we only need for the eigenvalues of $P^{-1}A$ to be clustered. Some approximate commutator methods, especially the "block diagonal" form typically used with MINRES for symmetric problems, intends to cluster the eigenvalues in three places. See this paper for further details. – Jed Brown Feb 12 '13 at 23:47

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.