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.


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.


Fisher's equation comparison of numerical methods.

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?

  • $\begingroup$ I recommend reading the basic error analysis of ODE solvers, for instance in Hairer/Nørsett/Wanner, plus some of the stability analysis. Most of your questions will be answered then. $\endgroup$ Aug 27, 2013 at 14:14
  • 1
    $\begingroup$ @boyfarrell, to avoid confusion of fellow readers, you should put the terminology right where explaining your method: 1. IMEX - explicit in the nonlinearity and implicit in the linear part. 2. this is the standard $\theta$-scheme, that will typically require Newton's method to solve for the update $\endgroup$
    – Jan
    Sep 2, 2013 at 9:44
  • 1
    $\begingroup$ Hello @Jan I think I got everything. Thanks again for your help. $\endgroup$
    – boyfarrell
    Sep 2, 2013 at 10:09

2 Answers 2


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.

  • $\begingroup$ Yes, I have applied standard central difference stencil to the diffusion term. I cannot use an explicit scheme (for the real problem I want to solve) because the stable time step is unrealistically small. This is why I am exploring IMEX or implicit options. Regarding your third point, to avoid the error accumulation I must use a multistep methods. Is the Crank-Nicolson scheme I have used above (with the Newton solver) classed as a multistep method (it has two points in time)? I was surprised the error increased with time when using the Newton solver method. $\endgroup$
    – boyfarrell
    Aug 27, 2013 at 13:12
  • $\begingroup$ Crank-Nicolson is a single-step method, as it writes as $u_h^{n+1}=u_h^n + \Phi(t^n,\tau^n,u_h^n,u_h^{n+1})$. Also, I don't see why multi-step schemes should avoid error accumulation. $\endgroup$
    – Jan
    Aug 27, 2013 at 13:22
  • 1
    $\begingroup$ OK thanks for explaining about the CN method. Yes, it is interesting why the multistep methods seems to have lower error accumulation. The reason that the Newton solver has error build up is because it is a single step method, I understand now. By the way, I know you like Python. I did all the above using scipy, numpy and matplotlib, gist.github.com/danieljfarrell/6353776 $\endgroup$
    – boyfarrell
    Aug 27, 2013 at 13:54
  • $\begingroup$ I have removed the link to the paper by Trefethen et. al. on high-order IMEX integration from my answer as there are better references to learn about IMEX schemes. $\endgroup$
    – Jan
    Sep 2, 2013 at 9:36

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

  • $\begingroup$ I am certainly approaching this from a physics perspective so all the technical details take some time for me to follow because I am not familiar with many of the terms. I'm actually an experimentalist! Would you explain a little more about order conditions? IMEX are these multistep methods mentioned by Jan? $\endgroup$
    – boyfarrell
    Aug 28, 2013 at 7:59
  • $\begingroup$ Order conditions are relationships between coefficients of ODE methods (e.g., entries in a Butcher tableau for Runge-Kutta methods) that must be satisfied to have an order of accuracy. Order conditions are discussed in any book or paper that designs ODE integration methods, but it basically amounts to repeatedly applying derivatives and matching terms in a Taylor expansion. The number of order conditions grows rapidly for high-order methods, which is why it becomes difficult to design high-order methods. Barriers are established by showing that order conditions are mutually incompatible. $\endgroup$
    – Jed Brown
    Aug 28, 2013 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.