Newton method for a nonlinear system of time-independent PDEs

Question

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:

Conceptual question

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?

Practical question

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?)

Answer 1

Conceptually, forget about the fact that the problem is time dependent for a moment. That's because if you apply a time stepping scheme to your time dependent PDE, you're going to get a (in your case nonlinear) PDE or a system of PDEs that solves for the solution $u^n(x)$ of the PDE at time step $n$ and where $u^{n-1}$, the solution of the previous time step, appears as known data. So whether or not the PDE is time dependent, ultimately you'll have to solve a nonlinear equation that's not time dependent.

Now, how to approach this stationary problem is a different question. I could try to write a lot of words on this, but let me instead just link to a bunch of lectures I've already recorded -- see lectures 31.5 and following here: http://www.math.colostate.edu/~bangerth/videos.html

Answer 2

Let me propose you a PDE example: $$-L(x,u)u = f(x) \qquad + \quad BC\tag{1}$$ Where $L$ is an operator depending on the spatial coordinates $x$ and the variable $u$.

Clearly the problem $(1)$ is nonlinear. This nonlinearity may be high or low, it depends on the form of the operator $L$.

Newton's method, linearises equation $(1)$ with $u^{n+1} = u^{n}+ \delta u$, leading to the following scheme: $$-L(x,u^n+\delta u)(u^n+\delta u) = f(x)$$ Expanding the latter in powers of $\delta u$ one obtains the equation for $\delta u$, provided $u^n$ $$ -L_u(x,u^n)\delta u =L(x,u^n)u^n+ f(x) \tag{2}$$ Where $L_u$ is the derivative of the operator w.r.t. $u$.

The scheme presented in $(2)$, if highly nonlinear, strongly depends (as you say) on the seed $u^0$, which must be as close as possible to the true solution of $(1)$. One smart option is solving the problem with the hypothesis $u\ll 1$: $$\lim_{u\to 0}\left[-L(x,u)\right]u^0 = f(x)$$

Other methods consist on solving a time-dependent version of $(1)$, asserting (strong hypothesis) that $$\lim_{t\to \infty}\frac{\partial u}{\partial t} = 0\tag{3}$$

Hypothesis $(3)$ is not always true. For example in the Navier-Stokes equation for the cylinder problem, from one specific Reynolds number, $(3)$ does not hold any more.

Therefore, if you are sure of that your problem will tend to stationary solution you can safely add the time derivative term to $(1)$: $$\frac{\partial u}{\partial t}-L(x,u)u = f(x) \qquad + \quad BC\quad+ \quad IC\tag{4}$$

The set-up in $(4)$ converges to the true solution iff you provide an initial condition close enough to the actual stationary solution (as for the seed of Newton's method, but now instead of seed it is identified with initial condition).

The easiest way to obtain a solution of $(4)$ is discretising the time derivative with a first order accuracy scheme. The rest of the terms belonging to the operator $L$ deserve a deeper discussion:

The high order derivatives (more that second order) must be implicitly treated (otherwise the timestep would be prohibitively small in order to solve the system as a set of quasi-isentropic states), I mean, you must solve them for the next time step, for example $\triangle u^{n+1}$, while the first order derivatives may be explicitly treated (stability issues arise when diffusive effects do not dominate or Courant number is small enough) as for example $u^n\textrm{grad}u^n$, or implicitly treated (stability issues do not srise but you pay the cost of the nonlinearity of the system you have to solve by means of e.g. Newton Raphson because it will converge extremely fast within 2 iterations) for example $u^{n+1}\textrm{grad}u^{n+1}$. Other options involve a "linearisation" including, instead, the term $u^n\textrm{grad}u^{n+1}$.

The safest method (a large time-step is required because, recall that we only want to reach a stationary solution) is the implicit scheme, given by the discretised version of $(4)$ $$\frac{ u^{n+1}-u^n}{\delta t}-L(x,u^{n+1})u^{n+1} = f(x) \tag{5}$$

The scheme $(5)$ converges more easily to a physical solution than $(2)$, once provided that $(3)$ holds for your problem. This is so because the previous $u^n$ used for the obtainment of $u^{n+1}$ has a physical meaning, on the contrary of what occurred in $(2)$.