Solving a nonlinear algebraic system that includes a linear term

Question

I am trying to solve a particular system of non linear equations written as $F(x) = 0$ in an efficient way.

More specifically, $$F(x) = (I - \gamma A)x - g(x) + C$$ where $\gamma$ is a scalarconstant, $C$ is a vector constant, $A$ is a constant matrix, and $g$ is a non-linear function of $R^n \rightarrow R^n$. I am using Newton's method and I am relatively satisfied with it. It involves computing the Jacobian $F'(x) = I - \gamma A - g'(x)$ and solving a linear system at each Newton iteration of the form $ F'(x^k) \Delta x^k = -F(x^k)$ and then set $x^{k+1} = x^k + \Delta x^k $. In this specific case it converges quickly (2-3 iterations at most) and everything works well.

However, as the problem becomes more complex and the size of my system of equations increases, computing the Jacobian and solving that linear system becomes a bottleneck. To overcome this, I use a simplified Newton's method which keeps F' constant over the whole Newton iteration process, and some parallel computing since the Jacobian can be computed in parallel efficiently in my case. Unfortunately, solving the linear system still remains a bottleneck (despite use of diverse parallel linear algebra techniques).

My question is then: Is there a way to take advantage of the specific form of $F$ ($i.e$ the fact that it is the sum of a linear operator and a non linear operator) ? In particular, I am especially interested in methods which would require me to solve linear systems with either $A$ (or $I - \gamma A$) or $g'$, but not their sum. Indeed, I have specific algorithms to solve linear systems with these matrices, but solving linear system with the sum of $I - \gamma A$ and $g'$ is very inefficient ($g'$ turns out to be block diagonal so adding $I - \gamma A$ messes it up, and solving $(I - \gamma A)X = B$ by itself is not an issue thanks to the specific structure of A)

Thanks

Answer 1

Please state up-front in your question that you are solving a reaction-diffusion problem. You should also be aware that time splitting methods, such as Strang Splitting, are very popular for these problems, though they can have serious accuracy and stability problems (e.g., Ropp, Shadid, and Ober, 2004). That said, you can use a multiplicative iterative method, preferably as a preconditioner to a linear or nonlinear Krylov method. To solve $(A + B) x = f$, compute

$$ \tilde x \gets x_0 + \omega_1 A^{-1} \big(f - (A+B) x_0\big) $$ $$ x_1 \gets \tilde x + \omega_2 B^{-1} \big(f - (A+B) \tilde x \big) $$

where $\omega_1,\omega_2$ are weighting parameters. The iteration matrix of this method is

$$ T = \big(1 - \omega_1 A^{-1}(A+B)\big) \big(1 - \omega_2 B^{-1}(A+B)\big) = (1-\omega_1 - \omega_1 A^{-1} B)(1 - \omega_2 - \omega_2 B^{-1} A). $$

The map is contractive if you choose $\omega_1,\omega_2$ such that $\lvert T \rvert < 1$, which may not be possible in case of indefinite reaction. One would typically use a Krylov method (and perhaps $\omega_1=\omega_2=1$) instead of tuning $\omega_1$ and $\omega_2$.

Answer 2

Is $g'$ by any change a low-rank correction to the remainder of the matrix? If so, you could think about using the Sherman-Morrison formula:

http://en.wikipedia.org/wiki/Sherman–Morrison_formula