I have taken a course on undergrad scientific computing which discussed nonlinear algebraic equations about half-way in, and PDEs at the very end, but never discussed nonlinear PDEs.
However in my research work, we have to use a solver that deals with a nonlinear coupled system of time-independent PDEs and I'm trying to have a very high level conceptual understanding of how Newton's method might have been applied in that solver internally, as well as a practical question of dealing with an unstable (stiff would be a better word?) scheme (Note: I'm familiar with FD/FV/FE and our solver uses FV for discretization, but I would use FD stencils for my discussion below). Kindly confirm or correct my line of thinking as follows:
For Newton's method applied to two nonlinear algebraic equations in two variables, the Jacobian matrix would be 2 by 2. However, if we have a coupled system of 2 nonlinear PDEs, let's say in 1-D (think a Poisson equation, and a continuity equation derived from Maxwell's equations, both nonlinearly coupled to each other in space, but otherwise time-independent),
First we would decide on the discretization and pick N points lying along the 1-D simulation domain. This will result in 2 times N algebraic equations (N algebraic equations for each of the 2 differential equations).
Each algebraic equation would have M independent variables, the value depending on the method of discretization, e.g., using 3-point FD stencil for 1-D, we have each algebraic equation in 3 variables? (except at the boundaries where the variables will be less than 3)
What would be the size of Jacobian? I guess the number of rows would be 2 times N, but the number of columns would depend on discretization? Again for 3-point FD, we'll have a 2-times-N by 3 Jacobian matrix? Or would it be 2-times-N by 2-times-N sparse Jacobian matrix? or something completely different? (edit 0: on further thought, I think it would be 2-times-N by N since we have 2-times-N algebraic equations but only N variables?)
If dimensionality of the system goes from 1-D to 2-D, and let's say we now use 5-point FD stencil, with N^2 discretization points. Would the Jacobian matrix be 2-times-N^2 by 5? or sparse 2-times-N^2 by 2-times-N^2?
Secondly, (somewhat separate from above), from a practical point of view, such a Newton scheme would be very unstable, especially if the PDEs are highly nonlinear in character, and therefore one way to mitigate this is to use, what our solver is using, a ramp-up Newton's method? where "ramp-up" is described as follows:
In a first step we ignore the given external Dirichlet BC, and use zero-everywhere BC (think Poisson equation). (As an aside, we could further split this into two steps by first doing, not only zero BC, but zero RHS, which turns it into a Laplace equation, then, in the second step, we iteratively ramp-up the RHS 1% at every step, similar to what I describe below).
And we solve one of the two equations first, e.g., only the Poisson equation.
Then we use the solution in a coupled solve, still zero-everywhere DBC.
Now we crank up the BC to 1% of target value, and do a coupled solve.
And we keep going, 1% at every step (call them steps of an outer loop), until we hit target DBC.
Is there a better name for this kind of ramp-up Newton scheme, or "guided Newton scheme", which feels like a forward Euler even though there is no time involved (it's not an initial value problem)?
Finally, a third question, if I try to solve the two PDEs in a "de-coupled" way, i.e., solve Poisson equation, use the solution as input to the continuity equation and solve it, and use that solution back into the next Poisson solve, that would make it a self-consistent scheme, which is different from Newton scheme?
Also a recommended reading on this topic that gives such a rundown would be highly appreciated. Thanks in advance.
Edit 1: Came across a very useful pdf shared here. I realized a better name for "guided Newton's method" is Newton's method with continuation (or numerical continuation). Plus there is other very useful information in the pdf. However I would still highly appreciate further comments.
Edit 2: I guess another question related to the second one above is, which of the two is a better answer for the reasoning behind using Newton's method with continuation instead of direct Newton's method:
Because high nonlinearity makes Newton's method unstable/stiff? (the typical explanation being, Newton's method works best if it starts very close to the solution, which means the initial guess for NM probably needs to be a solution of some non-NM solver that does not have that problem.)
Because high nonlinearity means solutions are not unique, and in a space of solutions, 99% are non-physical, and in order to reach the 1 physically meaningful solution, we have to "guide" the method we're using, through continuation.
(in other words, is it an issue with the method, unstable/stiff, or is it an issue with the system of PDEs, many non-unique solutions most of which are physically meaningless?)