NodePy version 0.4 released

NodePy is a Python package for analyzing numerical integrators for initial value ODEs.  It's essentially a collection of all the kinds of analysis I've used in my time integrator research, collected in a single object-oriented package.  

If you have a new Runge-Kutta method and want to know all about it, NodePy can tell you most anything.  If you want to design new time integration methods, NodePy can help you.

Although I'm rather proud of it, it fills a very small niche in the world and I'm not aware of anyone using it outside of my group and close collaborators.  If you've used it, please let me know in the comments.

One of the thorniest issues in NodePy previously was that floating-point representations of method coefficients were sometimes insufficient, especially when studying very high order methods.  I've now updated NodePy to use Sympy Rationals (and radicals, etc.) wherever possible, allowing exact analysis of many properties.

That and much more awaits in NodePy version 0.4, now available via pip.

Can a region of absolute stability be rectangular?

When a one-step integrator is applied to the solution of the linear scalar ODE

$$u'(t) = \lambda u(t)$$

the resulting iteration takes the simple form

$$U^{n} = R(h\lambda) U^{n-1}$$

where $U^n$ is a numerical approximation to the solution $u(t_n)$ and $R(h\lambda)$ is called the stability function.  The details of the stability function depend on the choice of numerical method, but for any explicit Runge-Kutta method, $R(z)$ is a polynomial whose degree is at most the number of stages of the method.

The stability function completely characterizes the accuracy and stability of the method when applied to linear problems.  Consider the first order linear, autonomous ODE

$$u'(t) = L u(t)$$

where now $u$ is a vector and $L$ is a square matrix.  The numerical solution will be

$$U^{n} = R(hL) U^{n-1}.$$

The global error satisfies a similar recurrence; in particular, it gets multiplied by a factor $R(hL)$ at each step.  Let $\lambda$ denote any eigenvalue of $L$; then If $L$ is a normal matrix, the solution will be absolutely stable in the Euclidean norm if all values $h\lambda$ lie within the stability region $S$, defined as

$$S = \{ z\in\mathbb{C} : |R(z)|\le 1\}.$$

Thus the region of absolute stability defines the portion of the complex plane in which a given numerical integration method may appropriately be applied.

In our preprint on Runge-Kutta stability regions, Aron Ahmadia and I claim that we have an algorithm to generate a stability region appropriate for any spectrum.  By considering high-degree polynomials, we find that the resulting stability regions are tightly adapted to the shape of the imposed spectrum.

While this promises to be very useful for some problems, it also has an aspect that's just fun: we can generate stability regions with unusual shapes.  I haven't explored this much yet, but a first question that we ask in the preprint is how to generate a stability region for a spectrum of eigenvalues forming a rectangle in the left half of the complex plane.

Here is an example of a resulting stability region:

The gray region is the set $S$ for a certain degree-20 stability polynomial corresponding to a consistent twenty-stage Runge-Kutta method.  As one colleague told me when I showed it to him, "this seems too good to be true; is that rectangle really the stability region?"

Indeed it is.  Zooming in on the top edge we see the detailed structure of the boundary:

Zooming in even closer:

As is typical with optimal stability polynomials, we se that the boundary is tangent or nearly tangent to the desired region at a large number of points (about 20 in this case).

What other shapes can be approximated?  More on that later...

Support the stackexchange for computational science

There is a new site on stackexchange for discussing computational science!  Please come participate -- there are already some great discussions.  The site is now in public beta, so you don't need an invitation.

In case you don't know about stackexchange, you may also be interested in the folowing:

• http://math.stackexchange.com/
• http://physics.stackexchange.com/

What journals do you read?

As a scientist, one is defined by the kind of problems one works on, the conferences one attends, the journals one publishes in, and the journals one reads. All of these except the last are more or less publicly available information.

Only you know precisely which journals you choose to read, yet they're an essential part of your scientific identity. They determine the kind of new advances you're likely to be aware of and where your research may turn in the future.

I've made my Mendeley library public, so anyone can see in great detail not only what journals I read but which articles I read. But most of you are probably not interested in quite that level of detail, so here's a list of the journals I follow closely. I collect their RSS feeds (with Google Reader) and read at least the title of every article they publish. I've grouped them into 3 main categories, but otherwise they're in no particular order. Those listed in bold are journals where I have published; they also tend to be the journals most heavily represented in my Mendeley library.

Numerical Analysis and scientific computing:

• BIT Numerical Mathematics
• Journal of Scientific Computing
• SIAM Journal on Scientific Computing (SISC)
• SIAM Journal on Numerical Analysis (SINUM)
• Mathematics of Computation
• Numerische Mathematik
• Applied Numerical Mathematics (APNUM)
• Computational Science and Discovery
• Journal of Computational and Applied Mathematics
• Journal of Computational Physics (JCP)
• Computing in Science and Engineering (CiSE)
• International Journal of Numerical Methods in Fluids (IJNMF)
• IMA Journal of Numerical Analysis
• ACM Transactions on Mathematical Software (TOMS)
• Computer Physics Communications
• Acta Numerica
• math.NA on ArXiv

 

I would add Communications in Computational Physics and Computational Methods in Applied Mathematics, but as far as I know they have no RSS feed.

Nonlinear waves:

• Physica D: Nonlinear Phenomena
• Nonlinearity

 

Here I would add Communications in the Mathematical Sciences, which also has no RSS feed.

General applied math:

• SIAM Review
• SIAM Journal on Applied Math
• IMA Journal of Applied Mathematics

 

Although ArXiv isn't a journal, I've included it here. Indeed, I find useful articles in that feed much more often than for most of the listed journals.

What does your list look like?

The positivity pipe dream fulfilled?

This post is about my recently accepted SINUM paper. This is an attempt to provide a broad context for the paper and related work in a less formal way, which couldn't be included in a journal publication.

In 1979, Bolley and Crouzeix published a fairly astonishing result (in French here). Namely, they showed that, for linear differential equations whose solution is always positive, it is impossible to design numerical methods that always yield a positive solution under arbitrarily large step sizes, even if one considers the very broad class of ODE solvers known as general linear methods. The only exception to this statement is the backward Euler method, which has many nice properties but is, unfortunately, too inaccurate for most applications.

This result stands in stark contrast to corresponding results on, say, stability in inner-product norms, where use of implicit methods can get you unconditional stability even for nonlinear problems.

I'm speaking roughly here, and won't attempt to be more precise (go read the paper if you want details). But the motivation cited by Bolley and Crouzeix for looking at this question was the idea that one might be able to take large time steps and maintain positivity in the solution of hyperbolic problems. Unfortunately, their result showed that this was not possible unless one was willing to settle for first-order accuracy.

What this result didn't indicate is how large the positivity-preserving step size could be for an implicit method. This question was partially answered within a decade, by Lenferink and by van de Griend & Kraaijevanger, for linear multistep methods and for Runge-Kutta methods, respectively. In both cases they found that the largest positivity-preserving step size was no more than two times the step size allowed by using the forward Euler method.

Of course, an implicit solve costs significantly more than an explicit solve, so gaining only a factor of two in the step size just isn't worth it. I don't know of any more work that was done on this problem after 1991 until fairly recently, in two of my own papers. The results there are consistent with the apparent barrier of a 2X step size, although see my short note at the end.

In one attempt to circumvent this bound, Gottlieb, Macdonald & Ruuth investigated the class of diagonally split Runge-Kutta methods (which aren't general linear methods, so escape the implications of Bolley & Crouzeix). They did find higher order methods with unconditional strong stability properties, but the accuracy of these methods invariably reduced to first order when used with large timesteps! It seemed that any effort to find unconditionally positive methods would be thwarted one way or another.

This sets the stage for my recently accepted SINUM paper. In this paper I found that by using both upwind-biased and downwind-biased discretizations (an idea that goes all the way back to Shu's original paper on "TVD time discretizations") in implicit Runge-Kutta methods, one can obtain second-order accurate methods that preserve positivity under arbitrarily large step sizes -- and they have this property even when applied to nonlinear problems. Remarkably, the methods appear to give quite accurate results when applied to problems with shocks.

It seems that we have gotten around the 2X barrier at last! But important issues remain, most notably the efficient implementation of these "implicit downwind Runge-Kutta schemes" in combination with high order hyperbolic PDE discretizations (like WENO).

I must reiterate that I've glossed over several important technical points here. For more details, go read the paper.

Note 1: In my Ph.D. thesis, I did find a method that breaks the 2X barrier, but only just barely and only for linear problems.

Note 2: I've referred rather carelessly in this post to literature on positivity preservation, contractivity, and strong stability preservation (monotonicity) without distinguishing between the three, simply because the conditions on the method turn out to be the same. The articles mentioned generally focus on one property or the other, but the results almost always carry over to all three.