I am trying to use Newton's method to get a stationary solution for a system of equations of the following form:

$$ \begin{Bmatrix} \frac{\partial x}{\partial t} \\ 0 \end{Bmatrix} = \begin{Bmatrix} f(x, y) \\ g(x, y) \end{Bmatrix} $$

For vectors $x$ and $y$.

Due to the nature of the equations I'm trying to solve, the Jacobian can be expressed in the form:

$$ J = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} $$

Such that $A$ and $C$ are dense matrices, while $B$ and $D$ are sparse. It should also be noted that the elements in the diagonals of $A$ and $C$ are much larger than those in $B$ and $D$ (by ~100 times).

Newton's method works perfectly when I solve the linear problem using Gaussian elimination, and the eigenvalues of the Jacobian are such that no severe numerical problems are faced with that algorithm.

I need to use an iterative method for larger problems, however, and the linear problem isn't appropriately solved with GMRES, IDR(s) or biconjugate gradient stabilized method when the matrix is taken as a whole, so I suppose another technique or precaution could be necessary.

I have tried preconditioning with an incomplete LU decomposition and with a diagonal preconditioner, and failed with both. When individually solving for blocks $A$ and $D$, however, one notices that the linear system $A\Delta x = - f(x, y)$ exhibits the same behaviour as the complete Jacobian, while $D \Delta y = - g(x, y)$ is easily solved by any of the Krylov subspace algorithms I cited above.

Is there any algorithm that could help solve the linear system?

  • $\begingroup$ Have you tried using Schur’s complement, the version involving the inverse of D? $\endgroup$ Commented Jun 15, 2021 at 19:33
  • $\begingroup$ @AmitHochman, in order to use Schur's complement, I'd have to compute the product $BD^{-1}C$, for which I'd have to invert D, right? I'd like to avoid that if possible, since I can't guarantee that D is tridiagonal (or of any other structure that allows for O(n) inversion) for all cases $\endgroup$ Commented Jun 15, 2021 at 20:07
  • $\begingroup$ You don’t have to invert D, but you would have to solve DX=B, or D^TX=C, whichever is easier. $\endgroup$ Commented Jun 15, 2021 at 20:24
  • 1
    $\begingroup$ There are methods, such as Krylov subspace recycling, which you could use to lower the cost of solving the same system with different right hand sides. $\endgroup$ Commented Jun 15, 2021 at 20:27
  • 1
    $\begingroup$ The form of your system is called in literature differential algebraic equation DAE. This might help in your literature search. $\endgroup$
    – Bort
    Commented Jun 16, 2021 at 12:25


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.