I am exploring the Method of Lines as a way of time stepping semi-discretised PDEs with ODE time-integration solvers. For an excellent introduction to this technique see the scholarpedia.org article. A popular ODE solver is VODE from http://www.netlib.org/ode/vode.f. It is also possible to access this solver by the scipy.integrate.ode set of solvers if you use Python.

I wish to solve an extremely stiff set of non-linear PDEs so I am interested in implicit methods because they offer stability and unrestricted time-step. I had planned to time-step the nonlinear system using a Newton iteration at every time step. However, I recently became aware of the Method of Lines (MOL).

The MOL technique seems a very powerful technique for time stepping semi-discretised PDEs because one can rely on robust and well tested ODE solver code that are well integrated with modern languages such as Python etc.

Question 1 How does the standard implicit Newton time-stepping method relate to the implicit Adams-Moulton technique used by the VODE solver?

To my current understanding the only difference is that Adams-Moulton is a generalisation of the implicit backwards Euler method which is include previous time points in the time integration up to a require order.

Question 2 Does the algorithm work like this?

  1. Solve the nonlinear system with $\theta$-method to give $u^{n+1}$, $$ w^{\prime}(u^{n+1}, u^{n}) = b_1 F(u^{n+1}) + b_0 F(u^{n}) $$

  2. Solve the system for the solution variable at the next time step $u^{n+2}$ but include the pervious value and appropriate coefficients to improve the time integration, $$ w^{\prime}(u^{n+2}, u^{n+1}) = b_2 F(u^{n+2}) + \underbrace{b_1 F(u^{n+1}) + b_0 F(u^{n})}_{\text{known}} $$

  3. Continue this processes up until the desired order ($N$) at which point we can start dropping the older terms, $$ w^{\prime}(u^{n+k}, u^{n+k-1}) = b_k F(u^{n+k}) + \underbrace{b_{k-1} F(u^{n+k-1}) + \,\,...\,\, +\, b_{k-N} F(u^{n+k-N})}_{\text{known}} $$

Question 3 This seems almost identical to the BDF approach, how does these approaches differ?

Question 4 For stiff nonlinear system of ODEs an implicit scheme (such as Newton iteration time-stepping) seems to be what people recommend. However, for a stiff ODE system why are Backwards Differentiation Formula (BDF) recommended and not Adams-Moulton?

Perhaps some of these equations don't make sense because I am confused regarding the basic technique, but I will be happy to be put straight.

  • $\begingroup$ Just as a warning: implicit does not imply A-stable, which is probably what you want for your method. In particular, no multistep method of order greater 2 is A-stable (Dahlquist barrier) $\endgroup$ Sep 2 '13 at 11:39
  • $\begingroup$ A-stability guarantees that the ODE tends to a solution at $t\rightarrow\infty$, so if the higher order method has problems I should use order 2. What about BDF? How do they compare? $\endgroup$
    – boyfarrell
    Sep 2 '13 at 12:36
  • $\begingroup$ If a variable-order method has good error control and heuristics, these concerns shouldn't matter. A-stability means you can take large time steps without worrying that the error will increase rapidly, but the error may still be unacceptably large with large time steps. Similarly, a high-order method may give you high accuracy, but stability will limit time step size. Good error heuristics should balance these two factors. $\endgroup$ Sep 2 '13 at 21:28
  • $\begingroup$ I see nothing in this question indicating whether $A$-stability is important and the above statements about $A$-stability are dubious/misleading. If anything "extremely stiff" usually refers to modes that decay rapidly due to eigenvalues near the negative real axis. In that case, a method need only be $A(\alpha)$-stable for small angle $\alpha$, in which case higher order BDFs are perfectly suitable. If you have eigenvalues near the imaginary axis, then you might want to check out RK-type methods, including Rosenbrock-W or ARK IMEX. $\endgroup$
    – Jed Brown
    Sep 5 '13 at 5:17

To answer your questions in order:

  1. Any implicit method for solving an ordinary differential equation involves solving a system of nonlinear equations. You can do this through variants of Newton's method, successive substitution, full approximation scheme, or any other approach that solves systems of nonlinear equations. (The caveat is, of course, depending on the problem, some methods work better than others, and some may not work at all. The case of successive substitution is a good example.)

  2. You are describing, more or less, how some implicit linear multistep methods work. Adams-Moulton and BDF methods are both linear multi step methods. For more background on linear multi step methods, check out the book by Ascher and Petzold, the two-volume work of Hairer and Wanner, Butcher's book, or the book by Lambert (all on numerical methods for ODEs).

  3. BDF methods have different coefficients than Adams-Moulton methods, and thus, they also have different stability and accuracy properties. The form of BDF methods is different than the equation you mentioned because, in standard form, many of its coefficients are zero. The Wikipedia articles on BDF methods and linear multistep methods illustrate this distinction nicely.

  4. If I remember correctly, BDF methods are recommended because they have large stability regions, but they're by no means the only choice for stiff systems. Rosenbrock methods, implicit Runge-Kutta methods, exponential integrators, and general linear methods can be effective as well, depending on the properties of your ODE right-hand side (and the spectrum of its Jacobian matrix), and whether or not you can effectively precondition any of the linear systems you solve during integration. If your stiff system can be decomposed into stiff and nonstiff components (like a reaction-diffusion equation), then implicit-explicit methods can also be worthwhile. The impression I get is that, due to the high-quality implementations of BDF methods out there, they've become the method of choice for people who just want something to work and don't really care about investigating alternatives or understanding their problem more deeply. They are reliable and accurate, but other algorithmic choices may be faster for a given problem.

Some final thoughts: VODE is really an old version of the CVODE solver in the SUNDIALS suite, and it might be worth using that instead of VODE to benefit from the algorithmic improvements made over the years (like incorporating inexact Newton methods into the nonlinear solves with good error control). Assimulo is a decent Python interface to SUNDIALS. PETSc also implements many of the methods I mentioned above and has a Python interface (called petsc4py), and Trilinos also implements some of those methods and has a Python interface (called PyTrilinos). If you don't think you need the new features, using scipy and VODE is probably good enough; lots of people still use it.

  • $\begingroup$ I disagree with the 2. point! The retrieval of starting values is a general problem in the application of multistep methods. If you do as the OP describes, the first step (2.1.; $\theta$-scheme) will limit the convergence order of the multi-step method by $2$ (if $\theta = 0.5$) or $1$ (else). It is recommended to compute the starting values with the same consistency order, e.g. via a single-step method. $\endgroup$
    – Jan
    Sep 2 '13 at 22:28
  • 2
    $\begingroup$ My understanding of linear multistep methods was that lower-order formulas within a family of methods could be used with small starting step sizes in adaptive schemes to generate values that bootstrap higher-order formulas in the same family of methods. I believe VODE uses this procedure; it would allow for significant code reuse, and would not introduce significant errors if the step size is small. VODE has starting step heuristics precisely for this purpose. $\endgroup$ Sep 3 '13 at 0:58
  • $\begingroup$ "In addition to adjusting the step size to meet the local error test, cvode periodically adjusts the order, with the goal of maximizing the step size. The integration starts out at order 1 and varies the order dynamically after that. The basic idea is to pick the order $q$ for which a polynomial of order $q$ best fits the discrete data involved in the multistep method.", page 8, CVODE user guide $\endgroup$ Sep 3 '13 at 1:54
  • $\begingroup$ I see, proper starting values can also be provided by means of a step size or order control. $\endgroup$
    – Jan
    Sep 3 '13 at 5:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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