Two coupled linear systems whose coefficient matrices depend on the solution of each other Consider the following coupled linear systems:
$A(y)x=p  \quad\quad\quad$ (1)
$B(x)y=q  \quad\quad\quad$ (2)
Here $p$ and $q$ are constant vectors, and $A$ and $B$ are symmetric and positive definite for all $x$ and $y$.
In practice, I found that I could solve this problem alternately, and the error seemed decreasing strictly. By solving alternately, I meant starting from $x_0$, solve $y_0$ by (2) with $B(x_0)$, and then solve $x_1$ by (1) with $A(y_0)$, then $y_1$ by (2)  with $B(x_1)$...
My question: does there exist any theory about a problem of this  kind?
Thanks in advance.
 A: Your question can be reformulated as a nonlinear system like this:
$$\begin{bmatrix}
    A(y) & 0\\
    0 & B(x) \\
\end{bmatrix} \begin{bmatrix}
    x \\
    y \\
\end{bmatrix}=\begin{bmatrix}
    p \\
    q \\
\end{bmatrix}$$
Now we can set:
$$ C=\begin{bmatrix}
    A(y) & 0\\
    0 & B(x) \\
\end{bmatrix}
 $$ $$
\mathbf{x}=\begin{bmatrix}
    x \\
    y \\
\end{bmatrix}$$
and
$$\mathbf{b}=\begin{bmatrix}
    p \\
    q \\
\end{bmatrix} $$
So your system can be written as:
$$C\mathbf{x}=\mathbf{b}$$
where
$$C=C(\mathbf{x})$$
or as
$$F(\mathbf{x})=C\mathbf{x}-\mathbf{b}=\mathbf{0}$$
The derivative, or rather Jacobian, of this function is the matrix:
$$J(\mathbf{x})=F'(\mathbf{x})$$
which in the linear case corresponds to
$$F'(\mathbf{x})=C$$
but will be a bit more complicated (chain rule) in your case.
You can solve this type of system iteratively by means of Newton's method:
$$\mathbf{x}_{n+1}=\mathbf{x}_{n}-J(\mathbf{x}_n)^{-1} F(\mathbf{x}_x)$$
However, it is in practice not a good idea to invert matrices and it is better to resort to stable numerical linear algebra methods. See also:
https://en.wikipedia.org/wiki/Newton%27s_method
Newton's method will converge quickly if you are close enough to the final solution with your starting guess.
The method you have used is instead a form of fixed point iteration method. See
https://en.wikipedia.org/wiki/Fixed-point_iteration
Such methods will only converge if the corresponding function/mapping is a contraction. See also:
https://www.cs.cornell.edu/~bindel/class/cs4220-s16/lec/2016-04-06-notes.pdf
There are usually many different ways of formulating a fixed point iteration. Some of them will lead to a contraction and converge to a solution.
In your case you can formulate a fixed-point iteration by writing:
$$\mathbf{x}=C(\mathbf{x})^{-1}\mathbf{b}$$
or iteratively
$$\mathbf{x}_{n+1}=C(\mathbf{x}_{n})^{-1}\mathbf{b}$$
However, again this is not the best way of solving. Better to resort to numerical linear algebra methods and instead view this as the linear system:
$$C(\mathbf{x}_{n})\mathbf{x}_{n+1}=\mathbf{b}$$
