Is it possible to solve nonlinear PDEs without using Newton-Raphson iteration?

Question

I am trying to understand some results and would appreciate some general comments on tackling nonlinear problems.

Fisher's equation (a nonlinear reaction-diffusion PDE),

$$ u_t = du_{xx} + \beta u (1 - u) = F(u) $$

in discretised form,

$$ u_j^{\prime} = \boldsymbol{L}\boldsymbol{u} + \beta u_j (1 - u_j) = F(\boldsymbol{u}) $$

where $\boldsymbol{L}$ is the differential operator and $\boldsymbol{u}=(u_{j-1}, u_j, u_{j+1}) $ is the discretisation stencil.

Method

I wish to apply a implicit scheme because I require stability and unrestricted time step. For this purpose I am using the $\theta$-method, (note that $\theta=1$ gives a fully implicit scheme and $\theta=0.5$ gives the trapezoidal or "Crank-Nicolson" scheme),

$$ u_{j}^{\prime} = \theta F(\boldsymbol{u}^{n+1}) + (1-\theta) F(\boldsymbol{u}^{n}) $$

However, for nonlinear problems this cannot be done because the equation cannot be written in a linear form.

To get around this problem I have been exploring two numerical approaches,

1. IMEX method

   $$ u_j^{\prime} = \underbrace{\theta\boldsymbol{L}\boldsymbol{u^{n+1}} + (1-\theta)\boldsymbol{L}\boldsymbol{u^{n}}}_{\theta-\text{method diffusion term}} + \underbrace{\beta u_j^{n} (1 - u_j^{n})}_{\text{Fully explicit reaction term}} $$

   The most obvious route is to ignore the nonlinear part of the reaction term and just update the reaction term with the best possible value, i.e. that from the previous time step. This results in the IMEX method.

2. Newton solver

$$ \nu^{k+1} = \nu^{k} - (I - \theta\tau A^{n})^{-1} \left( \nu^{k} - u_{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) \right) $$

The full $\theta$-method equation can be solved using the a Newton-Raphson iteration to find the future solution variable. Where $k$ is the iteration index ($k\geq0$) and $A^{n}$ is the Jacobian matrix of $F(w^n)$. Here I use the symbols $\nu^{k}$ for iteration variables such that they are distinguished from solution of the equation at a real time point $u^n$. This is actually a modified Newton solver because the Jacobian is not updated with every iteration.

Results

The results above are calculated for a reasonably large time step and they show the difference between the time stepping approach and a full Newton iteration solver.

Things I don't understand:

1. I am surprised that the time-stepping method does "OK" but it eventually lags behind the analytical solution as time goes by. (NB if I had chosen a smaller time-step then the time-stepping approach gives results closed to the analytical model). Why does the time-stepping approach give reasonable results to a nonlinear equation?

2. The Newton model does much better, but starts to lead the analytical model as time goes forward. Why does the accuracy of the Newton approach decrease with time? Can accuracy be improved?

3. Why is there a general feature that after many iterations then numerical model and the analytical model begin to diverge? Is this just because the time step is too large or will this always happen?

Answer 1

I assume, that you have conducted a space discretization, so that you are about solving the (vector-valued) ODE $$ \dot u_h(t) = F_h(t,u_h(t)), \text{ on [0,T] }, u_h(0) = \alpha. $$ via a numerical scheme $\Phi$ that advances the approximation $u_h^n$ at the current time instance $t=t^n$ to the next value $u_h^{n+1}$ at $t=t^{n+1}:=t^n+\tau$.

Then your questions refer to properties of explicit, where the update writes as $$ u_h^{n+1} = u_h^n + \Phi_e(t^n,\tau,u_h^n),\quad \quad \quad $$

implicit, written like $$ u_h^{n+1} = u_h^n + \Phi_i(t^n,\tau,u_h^{n+1},u_h^n), \quad \quad \quad (*) $$

or a combination of both ('IMEX', see @Jed Brown's answer) single-step time-stepping schemes.

In this setup, the Newton method simply is an approach to solve the possibly nonlinear in $u_h^{n+1}$ systems resulting from $(*)$.

And my answers base on results from the numerical analysis of single-step methods.

1. If you use convergent schemes, in terms of the convergence order, there is no general advantage of using implicit schemes (see. 2.). However, for stiff systems, e.g. your system containing a Laplacian, there are implicit schemes that are stable without time-step restrictions. Nevertheless, in theory, for the explicit scheme, you get better results with smaller time-steps, as long as your equation itself is stable (e.g., referring to Picard-Lindelof Theorem, if $F_h$ is Lipshitz in the second argument) and your time-step is not too small.
2. You can find examples, where explicit schemes perform better. (Theoretically, you can reverse the time in your example, start from the terminal value, and find implicit and explicit interchanged.) If you make the Newton error sufficiently small, you can still improve accuracy by decreasing the time-step or by using time-stepping schemes of higher order.
3. The constant $C$ in the error estimate for the global error grows exponentially with the length of the time-interval. See, e.g., here for the explicit Euler scheme. This is true for every single-step method. As the estimate is of type $err \approx C \tau^p$, $p>0$, a smaller time-step $\tau$ only postpones this effect.

---

Some more remarks and the final answer:

• IMEX schemes can be used to treat only the linear part implicitly what avoids the nonlinear solves. See Jed Brown's answer.
• Crank-Nicolson is a single-step method. The 'multi' in multi-step methods refers to the use of a number of preceding timesteps to define the current update. E.g. like $$ u_h^{n+1} = \Phi_m(t^n,\tau,u_h^{n+1},u_h^n,u_h^{n-1}). $$ This is very different from single-step and also split-step or IMEX methods, where the update is defined not recurring to previous values.

So, my answer is: Yes, you can solve nonlinear PDEs without Newton's method. You can use explicit schemes, 'IMEX' schemes, or socalled linearly implicit methods (e.g. the Rosenbrock methods). Also, you can employ other approaches to solve the systems from $(*)$ like fixed-point iteration or, in particular cases, algebraic solvers.

Answer 2

Short answer

If you only want second order accuracy and no embedded error estimation, chances are that you'll be happy with Strang splitting: half-step of reaction, full step of diffusion, half step of reaction.

Long answer

Reaction-diffusion, even with linear reaction, is famous for demonstrating splitting error. Indeed, it can be much worse, including "converging" to incorrect steady-states, mistaking steady-states for limit cycles, confuse stable and unstable configurations, and more. See Ropp, Shadid, and Ober (2004) and Knoll, Chacon, Margolin, and Mousseau (2003) for the computational physicists perspective on this. For the mathematician's analysis in terms of order conditions, see Hairer and Wanner's book on stiff ODE (Rosenbrock-W methods are a linearly-implicit IMEX method), Kennedy and Carpenter (2003) for nonlinearly-implicit IMEX "additive" Runge-Kutta, and Emil Constantinescu's page for more recent IMEX methods.

In general, IMEX methods have more order conditions than the underlying implicit and explicit methods alone. IMEX method pairs can be designed with desired linear and nonlinear stability and so that they satisfy all order conditions up to the design order of the method. Satisfying all order conditions will keep the asymptotic splitting error of the same scale as the error in each scheme separately. It says nothing about the pre-asymptotic regime (large time steps/low accuracy requirement), but it is rarely more stringent than resolution of each part separately. In any case, the splitting error is visible to the embedded error estimator (when using adaptive error control).

PETSc has many IMEX methods of the Rosenbrock-W and additive Runge-Kutta families, and will have extrapolation and linear multistep IMEX in our next release.

Disclaimer: I wrote much of the PETSc time integration support and collaborate with Emil (linked above).