Solving a system of nonlinear equations with an ODE solver is faster than with the Newton method?

Question

This is somehow unexpected, but my recent experience with solving a system of nonlinear equations is that treating them as the right hand side of a system of ordinary equations and then evolve the system with an ODE solver can be considerably faster than with the usual Newton-Raphson iteration.

Of course this only applies to some specific types of nonlinear equations. The equations I am dealing with is from a chemical equilibrium problem. The chemical evolution equation looks like $$\frac{d}{dt}x_i = f_i(\mathbf{x}) = -k_1x_i-k_2x_ix_j -k_3x_ix_k -\cdots + k_4x_j + k_5x_kx_h + \cdots,$$ $$\cdots$$ where the $k_i$s are all constants.

The goal is to solve the equations $$f_i(\mathbf{x}) = 0,\ \text{for }i=1,\ldots,n.$$ The approach I am talking about is to start from an initial guess $x_0$, and let it evolve until the $f_i$s are very close to $0$; an equilibrium solution is guaranteed to exist based on physical intuition (though not necessarily unique).

What I don't understand is that the ODE solver also has an Newton-Raphson iteration built-in, which is needed in an implicit scheme (which is the one I am using); how could the ODE approach be faster?

EDIT: The nonlinear equations have a trivial solution $\mathbf{x}=0$, but this is rarely the one we want. My understanding is that there is no guarantee the Newton iteration will converge to a ``physical'' solution starting from an initial estimate $\mathbf{x}_0$, while the ODE approach seems to have this property.

Update: I replaced one of the nonlinear equation with a conservation equation, and now the Newton method runs much faster (similar to roughly 4 times faster than the ODE approach), and it is likely that not doing this will yield the wrong result, since the system of nonlinear equations is not completely independent (which makes the Newton solver complain).

Answer 1

This is somehow unexpected, but my recent experience with solving a system of nonlinear equations is that treating them as the right hand side of a system of ordinary equations and then evolve the system with an ODE solver can be considerably faster than with the usual Newton-Raphson iteration.

It sounds like you're doing some sort of pseudotransient continuation (PtC). (See Is it well known that some optimization problems are equivalent to time-stepping?)

What I don't understand is that the ODE solver also has an Newton-Raphson iteration built-in, which is needed in an implicit scheme (which is the one I am using); how could the ODE approach be faster?

It could be that the sequence of initial guesses the PtC method goes through requires fewer iterations than a Newton-Raphson method for solving the algebraic equations you wish to solve. It would be helpful to give an example of the iterate history you're seeing.

My understanding is that there is no guarantee the Newton iteration will converge to a ``physical'' solution starting from an initial estimate x0, while the ODE approach seems to have this property.

Yes, there are no guarantees that the Newton iteration will converge to a "physical" solution. If the ODE approach converges, it should converge to a stable equilibrium point of the dynamical system, which, in your case, is the physical solution you are talking about.

I replaced one of the nonlinear equation with a conservation equation, and now the Newton method runs much faster (similar to the ODE approach), and it is likely that not doing this will yield the wrong result, since the system of nonlinear equations is not completely independent (which makes the Newton solver complain).

If the Jacobian is rank-deficient, then the linear solves may or may not return solutions. For instance, a Krylov subspace solver will typically calculate a Drazin inverse, which yields a least-squares solution. Gaussian elimination may throw an error due to a zero pivot or something like that; a least-squares solution could also be calculated using QR or SVD.

The Jacobian matrices for time steps in the PtC method will usually be full rank; they will be of the form $(I - \Delta{t}J)$, where $J$ is the Jacobian of the left-hand side of $f(x) = 0$. Rank-deficiency is probably another reason the convergence of the original Newton-Raphson method you used wasn't so great.

The conservation equation would yield a DAE, which you could also use in a PtC approach, if you like. See Coffey, Kelley, and Keyes (2003) for analysis of PtC for DAE compared to ODE formulations.