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?
