Recurrence, Averaging and Predictability

Motivation and Background

Yet another installment in the Lorenz63 series. This time motivated by a commenter on Climate Etc. Tomas Milanovic claims that time averages are chaotic too in response to the oft repeated claim that the predictability limitations of nonlinear dynamical systems are not a problem in the case of climate prediction. Lorenz would seem to agree, “most climatic elements, and certainly climatic means, are not predictable in the first sense at infinite range, since a non-periodic series cannot be made periodic through averaging [1].” We’re not going to just take his word on it. We’ll see if we can demonstrate this with our toy model.

That’s the motivation, but before we get to toy model results a little background discussion is in order. In this previous entry I illustrated the different types of functionals that you might be interested in depending on whether you are doing weather prediction or climate prediction. I also made the remark, “A climate prediction is trying to provide a predictive distribution of a time-averaged atmospheric state which is (hopefully) independent of time far enough into the future.” It was pointed out to me that this is a testable hypothesis [2], and that the empirical evidence doesn’t seem to support the existence of time-averages (or other functionals) describing the Earth’s climate system that are independent of time [3]. In fact, the above assumption was critiqued by none other than Lorenz in 1968 [4]. In that paper he states,

Questions concerning the existence and uniqueness of long-term statistics fall into the realm of ergodic theory. [...] In the case of nonlinear equations, the uniqueness of long-term statistics is not assured. From the way in which the problem is formulated, the system of equations, expressed in deterministic form, together with a specified set of initial conditions, determines a time-dependent solution extending indefinitely into the future, and therefore determines a set of long-term statistics. The question remains as to whether such statistics are independent of the choice of initial conditions.

He goes on to define a system as transitive if the long-term statistics are independent of initial condition, and intransitive if there are “two or more sets of long-term statistics, each of which has a greater-than-zero probability of resulting from randomly chosen initial conditions.” Since the concept of climate change has no meaning for statistics over infinitely long intervals, he then defines a system as almost intransitive if the statistics at infinity are unique, but the statistics over finite intervals depend (perhaps even sensitively) on initial conditions. In the context of policy relevance we are generally interested in behavior over finite time-intervals.

In fact, from what I’ve been able to find, different large-scale spatial averages (or coherent structures, which you could track by suitable projections or filtering) of state for the climate system face similar limits to predictability as un-averaged states. The predictability just decays at a slower rate. So instead of predictive limitations for weather-like functionals on the order of a few weeks, the more climate-like functionals become unpredictable on slower time-scales. There’s no magic here, things don’t suddenly become predictable a couple decades or a century hence because you take an average. It’s just that averaging or filtering may change the rate that errors for that functional grow (because in spatio-temporal chaos different structures, or state vectors, will have different error growth rates and reach saturation at different times). Again Lorenz puts it well, “the theory which assures us of ultimate decay of atmospheric predictability says nothing about the rate of decay” [1]. Recent work shows that initialization matters for decadal prediction, and that the predictability of various functionals decay at different rates [5]. For instance, sea surface temperature anomalies are predictable at longer forecast horizons than surface temperatures over land. Hind-casts of large spatial averages on decadal time-scales have shown skill in the last two decades of the past century (though they had trouble beating a persistence forecast for much of the rest of the century) [6].

I’ve noticed in on-line discussions about climate science that some people think that the problem of establishing long term statistics for nonlinear systems is a solved one. That is not the case for the complex, nonlinear systems we are generally most interested in (there are results for our toy though [7, 8]). I think this snippet sums things up well,

Atmospheric and oceanic forcings are strongest at global equilibrium scales of 10⁷ m and seasons to millennia. Fluid mixing and dissipation occur at micorscales of 10^-3 m and 10^-3s, and cloud particulate transformations happen at 10^-6 m or smaller. Observed intrinsic variability is spectrally broad band across all intermediate scales. A full representation for all dynamical degrees of freedom in different quantities and scales is uncomputable even with optimistically foreseeable computer technology. No fundamentally reliable reduction of the size of the AOS [atmospheric oceanic simulation] dynamical system (i.e., a statistical mechanics analogous to the transition between molecular kinetics and fluid dynamics) is yet envisioned. [9]

Here McWilliams is making a point similar to that made by Lorenz in [4] about establishing a statistical mechanics for climate. This would be great if it happened, because that would mean that the problem of turbulence would be solved for us engineers too. Right now the best we have (engineers interested in turbulent flows and climate scientists too) is empirically adequate models that are calibrated to work well in specific corners of reality.

Lorenz was responsible for another useful concept concerning predictability, that is predictability of the first and second kind [1]. If you care about the time-accurate evolution of the order of states then you are interested in predictability of the first kind. If, however, you do not care about the order, but only the statistics, then you are concerned with predictability of the second kind. Unfortunately, Lorenz’s concepts of first and second kind predictability have been morphed in to a claim that first kind predictability is about solving initial value problem (IVP)s and second kind predictability is about solving boundary value problem (BVP)s. For example, “Predictability of the second kind focuses on the boundary value problem: how predictable changes in the boundary conditions that affect climate can provide predictive power [5].” This is unsound. If you read Lorenz closely, you’ll see that the important open question he was exploring about whether the climate is transitive, intransitive or almost intransitive has been assumed away by the spurious association of kinds of predictability with kinds of problems [1]. Lorenz never made this mistake, he was always clear that the difference in kinds of predictability depends on the functionals you are interested in, not whether it is appropriate to solve an IVP or a BVP (what reason could you have for expecting meaningful frequency statistics from a solution to a BVP?). Those considerations depend on the sort of system you have. In an intransitive or almost intransitive system even climate-like functionals depend on the initial conditions.

A good early paper on applying information theory concepts to climate predictability is by Leung and North [10], and there is a more recent review article that covers the basic concepts by DelSole and Tippett [11].

Recurrence Plots

Recurrence plots are useful for getting a quick qualitative feel for the type of response exhibited by a time-series [12, 13]. First we run a little initial condition (IC) ensemble with our toy model. The computer experiment we’ll run to explore this question will consist of perturbations to the initial conditions (I chose the size of the perturbation so the ensemble would blow-up around t = 12). Rather than sampling from a distribution for the members of the ensemble, I chose them according a stochastic collocation (this helps in getting the same results every time too).

---

(a)EnsembleTrajectories (b)EnsembleMean Figure 1: Initial Condition Ensemble

---

One thing that these two plots makes clear is that it doesn’t make much sense to compare individual trajectories with the ensemble mean. The mean is a parameter of a distribution describing a population of which the trajectories are members. While the trajectories are all orbits on the attractor, the mean is not.

---

(a)SingleTrajectory (b)EnsembleMean Figure 2: Chaotic Recurrence Plots

---

Comparing the chaotic recurrence plots with the plots below of a periodic series and a stochastic series illustrates the qualitative differences in appearance.

---

(a)PeriodicSeries (b)StochasticSeries Figure 3: Non-chaotic Recurrence Plots

---

Clearly, both the ensemble mean and the individual trajectory are chaotic series, sort of “between” periodic and stochastic in their appearance. Ensemble averaging doesn’t make our chaotic series non-chaotic, what about time averaging?

Predictability Decay

How does averaging affect the decay of predictability for the state of the Lorenz63 system, and can we measure this effect? We can track how the predictability of the future state decays given knowledge of the initial state by using the relative entropy. There are other choices for measures such as mutual information [10]. Since we’ve already got our ensemble though, we can just use entropy like we did before. Rather than just a simple moving average, I’ll be calculating an exponentially weighted one using an FFT-based approach, of course (there’s some edge effects we’d need to worry about if this were a serious analysis, but we’ll ignore that for now). The entropy for the ensemble is shown for three different smoothing levels in Figure 4 (the high entropy prior to t = 5 for the smoothed series is spurious because I didn’t pad the series and it’s calculated with the FFT).

---

Figure 4: Entropy of Exponentially Weighted Smoothed Series

---

While smoothing does lower the entropy of the ensemble (lower entropy for more smoothing / smaller λ), it still experiences the same sort of “blow-up” as the unsmoothed trajectory. This indicates problems for predictability even for our time-averaged functionals. Guess what? The recurrence plot indicates that our smoothed trajectory is still chaotic!

---

Figure 5: Smoothed Trajectory Recurrence Plot

---

This result shouldn't be too surprising, moving averages or smoothing (of whatever type you fancy) are linear operations. It would probably take a pretty clever nonlinear transformation to turn a chaotic series into a non-chaotic one (think about how the series in this case is generated in the first place). I wouldn't expect any combination of linear transformations to accomplish that.

Conclusions

I’ll begin the end with another great point from McWilliams (though I’ve not heard of sub-grid fluctuations referred to as “computational noise,” that term makes me think of round-off error) that should serve to temper our demands of predictive capability from climate models[9]:

Among their other roles, parametrizations regularize the solutions on the grid scale by limiting fine-scale variance (also known as computational noise). This practice makes the choices of discrete algorithms quite influential on the results, and it removes the simulation from the mathematically preferable realm of asymptotic convergence with resolution, in which the results are independent of resolution and all well conceived algorithms yield the same answer.

If I had read this earlier, I wouldn’t have spent so much time searching for something that doesn’t exist.

Regardless of my tortured learning process, what do the toy models tell us? Our ability to predict the future is fundamentally limited. Not really an earth-shattering discovery; it seems a whole lot like common sense. Does this have any implication for how we make decisions? I think it does. Our choices should be robust with respect to these inescapable limitations. In engineering we look for broad optimums that are insensitive to design or requirements uncertainties. The same sort of design thinking applies to strategic decision making or policy design. The fundamental truism for us to remember in trying to make good decisions under the uncertainty caused by practical and theoretical constraints is that limits on predictability do not imply impotence.

References

[1]   Lorenz, E. N., The Physical Basis of Climate and Climate Modeling, Vol. 16 of GARP publication series, chap. Climatic Predictability, World Meteorological Organization, 1975, pp. 132–136.

[2]   Pielke Sr, R. A., “your query,” September 2010, electronic mail to the author.

[3]   Rial, J. A., Pielke Sr, R. A., Beniston, M., Claussen, M., Canadell, J., Cox, P., Held, H., Noblet-Ducoudr, N. D., Prinn, R., Reynolds, J. F., and Salas, J. D., “Nonlinearities, Feedbacks And Critical Thresholds Within The EarthS Climate System,” Climatic Change, Vol. 65, No. 1-2, 2004, pp. 11–38.

[4]   Lorenz, E. N., “Climatic Determinism,” Meteorological Monographs, Vol. 8, No. 30, 1968.

[5]   Collins, M. and Allen, M. R., “Assessing The Relative Roles Of Initial And Boundary Conditions In Interannual To Decadal Climate Predictability,” Journal ofClimate, Vol. 15, No. 21, 2002, pp. 3104–3109.

[6]   Lee, T. C., Zwiers, F. W., Zhang, X., and Tsao, M., “Evidence of Decadal Climate Prediction Skill Resulting from Changes in Anthropogenic Forcing,” Journal of Climate, Vol. 19, 2006.

[7]   Tucker, W., The Lorenz Attractor Exists, Ph.D. thesis, Uppsala University, 1998.

[8]   Kehlet, B. and Logg, A., “Long-Time Computability of the Lorenz System,”http://lorenzsystem.net/.

[9]   McWilliams, J. C., “Irreducible Imprecision In Atmospheric And Oceanic Simulations,” Vol. 104 of National Academy of Sciences, National Academy of Sciences, pp. 8709 – 8713.

[10]   Leung, L.-Y. and North, G. R., “Information Theory and Climate Prediction,”Journal of Climate, Vol. 3, 1990, pp. 5–14.

[11]   DelSole, T. and Tippett, M. K., “Predictability: Recent insights from information theory,” Reviews of Geophysics, Vol. 45, 2007.

[12]   Eckmann, J.-P., Kamphorst, S. O., and Ruelle, D., “Recurrence Plots of Dynamical Systems,” EPL (Europhysics Letters), Vol. 4, No. 9, 1987, pp. 973.

[13]   Marwan, N., “A historical review of recurrence plots,” The European PhysicalJournal - Special Topics, Vol. 164, 2008, pp. 3–12, 10.1140/epjst/e2008-00829-1.

Predictions and Entropy in Ensembles

This post is another installment in the Lorenz63 series. In it I’ll try to address a common confusion among folks trying to understand what uncertainty in initial conditions means for our ability to predict the future state of a system (in particular as it pertains to climate change). The common pedagogical technique is to describe weather prediction as an initial value problem (IVP) and climate prediction as a boundary value problem (BVP) (for example see Serendipity). I don’t think that happens to be a very good teaching technique (it seems too hand-wavy to me [update: better reasoning given down in this comment and this comment]). I think there are probably better approaches which would give more insight into the problems. I’m a learn by doing / example kind of guy, so that’s what this post will focus on: using the Lorenz ’63 system as a useful toy to give insight into the problem.

In fact, predictions of both climate and weather often use the same models which approximately solve the same initial-boundary value problem described by partial differential equations (PDE)s for the conservation of mass, momentum and energy (along with many parametrizations for physical and chemical processes in the atmosphere as well as sub-grid scale models of unresolved flow features). The real distinction is that climate researchers and weather forecasters care about different functions or statistics of the solutions over different time-scales. A weather forecast depends on providing time-accurate predictive distributions of the state of the atmosphere in the near-future. A climate prediction is trying to provide a predictive distribution of a time-averaged atmospheric state which is (hopefully) independent of time far enough into the future (there’s an implicit ergodic hypothesis here, which, as reader Tom Vonk points out, still requires some theoretical developments to justify for the PDEs we’re actually interested in).

Another concept I’d like to introduce before I show example results from the Lorenz ’63 system is the entropy of a probability distribution. This can be viewed as a measure of informativeness of the distribution. [update: this paper presents the idea in context of climatic predictions] So a very informative distribution would have low entropy, but an uninformative distribution would have maximum entropy. For weather forecasting we would like low entropy in our predictive distributions, because that means we have significant information about what is going to happen. For climate, the distribution itself is what we want, and we are actually in some sense looking for the maximum entropy distribution that is consistent with our constraints. Measuring the “distance” between distributions with different constraints is what climate forecasting is all about.

Now for the toying with the Lorenz63 system. The “ensembles” I’m running are just varying the initial condition in the x-component (using the same approach shown here). I’m also using two slightly different b_x parameters in the forcing function. Figure 1 shows the x-component of the resulting trajectories.

---

(a)Forcing1 (b)Forcing2 Figure 1: X-component trajectory under different forcings

---

The analogy to the weather / climate difference is pretty well illustrated by these results. Up to t ~ 3 we could do some pretty decent “weather” prediction, the ensemble members stay close together. After that things diverge abruptly (exponential growth in initial differences), this illustrates the need to “spin-up” a model when you are interested in the climate distributions. After t ~ 10 we could probably start estimating the “climate” of the two different forcings (histograms for 10 < t < 12 shown in Figure 2).

---

Figure 2: Two different “climate” distributions

---

These trajectories illustrate another interesting aspect of deterministic chaos, our uncertainty in the future state does not increase monotonically, it will grow and shrink in time (for instance compare the spread in the ensemble members at t = 4 to that at t = 8 shown in Figure 1). The entropy in the distribution of the trajectories as a function of time is shown in Figure 3. A Python function to calculate this for my ensemble (stored in a 2d array) is shown below.

def ensemble_entropy(x): 
    # the ensemble runs accross the first dimension of x, the time 
    # runs accross the second, return entropy as a function of time 
    eps = 1e-16 
    ent = sp.zeros(x.shape[1], dtype=float) 
    bin_edges = sp.linspace( 
        min(x.ravel()), max(x.ravel()), int(sp.sqrt(x.shape[0]))) 
    for i in xrange(x.shape[1]): 
        # the histogram function returns a probability density with normed=True 
        p = sp.histogram(x[:,i], bins=bin_edges, normed=True) 
        # we would like a probability mass for each bin, so we need to 
        # multiply by the width of the bin: 
        dx = p[1][1] - p[1][0] 
        p = dx * p[0] 
        # normalize (it’s generally very close, this is probably 
        # unnecessary), and take care of zero p bins so we don’t get 
        # NaNs in the log: 
        p = p / sum(p) + eps 
        ent[i] = -sum(p * sp.log(p)) 
    return(ent)

---

Figure 3: Trajectory-Distribution Entropy as a Function of Time

---

The entropy gives us a nice measure of the informativeness of our ensemble. In the initial stages (t < 4) we’ve got small entropy (we could make “weather” predictions here). There’s a significant spike around t = 4, and then we see the magnitude of entropy drop off for a bit (or a nat, ha-ha) around t = 8, which matches the eye-balling of the ensemble spread we did earlier.

That’s it for toy model results, now for some conclusions and opinions.

There are two honest concerns with climate forecasting (if you know of more let me hear about them, if you don’t think my concerns are honest, let me hear that too). First, are the things we can predict with climate modeling useful for planing mitigation and adaption policies? So many of the alarming predictions of costs and catastrophes attributed to climate change in the press (and even in the IPCC’s reports) depend on particular regional (rather than global) climate changes. I think an open research question is how the time averaging (and large time-steps) involved in calculating the long trajectories for estimating equilibrium climate distributions (setting aside the theoretical underpinning of this ergodic assumption) affect the accuracy of the predicted spatial variations and regional changes (it is after-all a PDE rather than an ordinary differential equations (ODE)). This seems to be an area of research that is just beginning. Also, the papers I’ve been able to find so far don’t seem to report any grid convergence index results for the solutions (please link to papers that do have this in the comments if you know of any, thanks). This is an important part of what Roache calls calculation verification (as opposed to the code verification demonstrated here).

Second, what about the details of the ’averaging window’ (in both space and time)? How useful to policy are the results of long time-averages? Are the equilibrium distributions things we will ever actually reach? These two concerns about the usefulness of the climate modeling product for policy makers (and in a Republic like mine, the public), and the details about the averaging, seem to be the motivation for Pielke Sr.’s quibble over in this thread about the definition of climate and the implications of chaos. As Pielke points out, take your averaging volume small enough, and things start looking pretty chaotic.

My personal view is that deterministic chaos is neat in toy problems (and a fun challenge for applying the method of manufactured solutions), but in the real world our uncertainties about everything (and the stochastic nature of many of the forcings) swamp the infinitesimals. The thing that makes the dynamics of the climate-policy-science “system” interesting is the tension between giving useful insight to decision makers and ensuring that insight is not overly sensitive to our inescapable uncertainties. Right now, the state of that system is far from equilibrium.

[Update: These survey results (courtesy of Roger Pielke Sr) are interesting.  They seem to indicate that "climate scientists" view climate prediction as an IVP.
The question is, 16. How would you rate the ability of global climate models to:  (very poor 1  2  3  4  5  6 very good) 16c. model temperature values for the next 10 years 16d. model temperature values for the next 50 years.  The mean response for the longer term prediction is actually lower than for the short term prediction (3.7 vs. 4.2), though the difference isn't that big considering the standard deviation.

The comparing the "reproduce observations" questions (16a. and 16b.) with the "predict future values" (16c. and 16d.) goes to the validation question.  Would you bet your life on your models predictions?]

Easy Full-Factorial Ensemble Experiments

This is another installment of the Lorenz63 series. In this post I wanted to demonstrate an ’experiment’ on the Lorenz ’63 system. The basics of these sorts of computer experiments as they relate to climate modeling are covered at climateprediction.net.

The main factors underlying our experimental design are initial conditions, parameters and forcing. Since we followed Roache’s advice about making our code Verifiable with the method of manufactured solutions (MMS) the ability to specify forcings is already included. Python has some nice built-in functions in the itertools module that make doing a computational experiment really easy. First all of the arguments for the time integrator are defined as lists (the ones which aren’t part of the experiment only have one element, so they remain constant):

# factors in our experimental design (all of the inputs to the system 
# integrator, some only have one level): 
nt = 2000 
savevery = [10] 
nsaves = nt / savevery[0] + 1 
T = 2.0 
# the first element of the array is the initial condition, the rest 
# are over-written with the solution: 
x = [1.05*sp.ones(nsaves,dtype=float,order=’F’), 
     0.95*sp.ones(nsaves,dtype=float,order=’F’)] 
y = [-1.05*sp.ones(nsaves,dtype=float,order=’F’), 
      -0.95*sp.ones(nsaves,dtype=float,order=’F’)] 
z = [-10.5*sp.ones(nsaves,dtype=float,order=’F’), 
         -9.5*sp.ones(nsaves,dtype=float,order=’F’)] 
t = [sp.zeros(nsaves, dtype=float, order=’F’)] 
Ra = [27.0, 28.0, 29.0] 
Pr = [9.0, 10.0, 11.0] 
b = [8.0/3.0 - 0.1, 8.0/3.0 + 0.1] 
h = [T / float(nt)] 
tol = [1e-14] 
maxits = [20] 
forcing = [0, 1] 
bx = [0.0] 
by = [-1.0] 
bz = [-10.0]

Then a call to itertools.product() gives us the Cartesian Product of all of those lists, which is a full-factorial set of the inputs (2 ⋅ 2 ⋅ 2 ⋅ 3 ⋅ 3 ⋅ 2 ⋅ 2 = 288 treatments in this example).

# this gives us an iterator of every combination of our factor levels: 
full_factorial = itertools.product( 
        x, y, z, t, Ra, Pr, b, h, tol, maxits, forcing, bx, by, bz, savevery)

Once we have the iterator defined, we can calculate all of the trajectories associated with the input treatments in the experiment:

# iterate over all of the ’treatments’ in our design, the iterator 
# yields the arguments directly: 
for args in full_factorial: 
    L63.time_loop(*args) 
    xresult.append(args[0].copy()) 
    yresult.append(args[1].copy()) 
    zresult.append(args[2].copy()) 
    tresult.append(args[3].copy())

The resulting trajectories are shown in Figure 1. It’s pretty clear that realistic uncertainties in parameters and initial conditions lead quickly to pretty diffuse predictive distributions. Part of the challenge of attribution studies is determining the effect of each change in parameter or forcing. In realistic models, the number of ’factors’ to vary is quite large.

---

(a) X component (b) Y component (c) Z component Figure 1: Trajectories for the Full-Factorial Experiment

---

With such a low-dimension system, and few parameters it is easy enough to run the whole full-factorial, smarter experimental designs are required for experiments based on large simulation runs. A good reference on this is Design and Analysis of Computer Experiments.

Chaotic Time-Convergence and MMS

This previous post covered time-convergence of a single trajectory, and the growth of the spread in an ensemble. The next post covered some special considerations that chaos imposes on us when we try to apply the method of manufactured solutions (MMS). I’ll try to tie up some loose ends in this (perhaps final) installment of the Lorenz63 series. Two questions that I think I haven’t adequately addressed yet:

• What about the time-step convergence behavior of an ensemble?
• So can you successfully apply the MMS to verify implementations of numerical approximations to this system or not? (The literature indicates that you can, but I didn’t actually demonstrate the correctness of my implementation.)

To address the first question I wrote a little Python script that runs a range of trajectories with various time-steps and also perturbations to the initial conditions:

nt = [192001, 384001, 768001] # number of time-steps 
step = [1, 2, 4] 
T = 60.0 # integration period 
# parameters for the Lorentz 63 system to get chaotic response: 
Pr = 10.0 
Ra = 28.0 
b = 8.0 / 3.0 
# knobs for the Newton’s method: 
tol = 1e-14 
maxits = 20 
# number of different IC’s in our ensemble: 
nic = 5 
# number of pairs of differences between ensemble members: 
npairs = comb(nic, 2, exact=1) 
h = [] 
x = [] 
xdiffs = [] 
y = [] 
ydiffs = [] 
z = [] 
zdiffs = [] 
t = [] 
for n in nt: 
    h.append(T / float(n)) 
    x.append(sp.zeros((n,nic), dtype=float, order=’F’)) 
    xdiffs.append(sp.zeros((n,npairs), dtype=float, order=’F’)) 
    y.append(sp.zeros((n,nic), dtype=float, order=’F’)) 
    ydiffs.append(sp.zeros((n,npairs), dtype=float, order=’F’)) 
    z.append(sp.zeros((n,nic), dtype=float, order=’F’)) 
    zdiffs.append(sp.zeros((n,npairs), dtype=float, order=’F’)) 
    t.append(sp.linspace(0.0, T, n)) 
 
# some perturbations: 
eps = 1e-16 * sp.array([0, 1, -1, 2, -2]) 
for i in xrange(len(nt)): 
    x[i][0,:] = 1.0 + eps 
    y[i][0,:] = -1.0 + eps 
    z[i][0,:] = 10.0 + eps 
 
# calculate all of the trajectories: 
for i in xrange(len(nt)): # time-step variation loop 
    for j in xrange(nic): # perturbations loop 
        L63.time_loop( 
            x[i][:,j],y[i][:,j],z[i][:,j],t[i],Ra,Pr,b,h[i],tol,maxits)

Here the result for the x-component of the solution:

---

Figure 1: Convergence and Chaotic spread behavior of Lorenz ensemble

---

Figure 1 shows that we have a time-step converged solution out to t ≈ 30. After which, the trajectory between the different time-steps is significantly different, but the ensemble spread is still negligible. At t ≈ 40 the ensemble begins to spread noticeably for all three time-step sizes, and then shortly thereafter we get an abrupt divergence of the trajectories for all of the ensemble members.

The interesting question is whether the divergence is inevitable (that’s my impression from reading the ’chaos’ literature, let me know if you think/know something different), or if getting more and more accurate solutions will delay that abrupt divergence (the visible divergence seems to happen latter in the most resolved ensemble). Unfortunately, I’m already flirting with swap-thrash on my little laptop with these ensembles, so I need to do a little re-factoring before I’ll be able to demonstrate that one definitively (so this probably isn’t the final installment). What’s really required here is a time-step converged solution of the region t > 40 where the ensemble members diverge. These solutions aren’t there yet (smaller time-step or higher-order method is needed).

The second part of the post will deal with demonstrating an MMS-based verification of my (purportedly third-order accurate) implementation. As discussed in this post, since the system is chaotic, choosing a manufactured solution becomes a bit more difficult. Here’s a Python script which calculates trajectories with MMS forcing and calculates error norms:

nt = [501, 1001, 2001, 4001, 8001] # number of time-steps 
step = [1, 2, 4, 8, 16] 
T = 2.0 * sp.pi # integration period 
Pr = 10.0 
Ra = 28.0 
b = 8.0 / 3.0 
tol = 1e-14 
maxits = 20 
mms_switch = 1 
h = [] 
x = [] 
y = [] 
z = [] 
t = [] 
ms = [] 
for n in nt: 
    h.append(T / float(n)) # time-step size 
    x.append(sp.zeros(n, dtype=float, order=’F’)) 
    y.append(sp.zeros(n, dtype=float, order=’F’)) 
    z.append(sp.zeros(n, dtype=float, order=’F’)) 
    t.append(sp.linspace(0.0, T, n)) 
    ms.append(sp.zeros((n,3), dtype=float, order=’F’)) 
 
bx = 0.0 
by = -1.0 
bz = 10.0 
# manufactured solutions: 
for i in xrange(len(nt)): 
    ms[i][:,0] = sp.cos(t[i]) + bx 
    ms[i][:,1] = sp.sin(t[i]) + by 
    ms[i][:,2] = sp.sin(t[i]) + bz 
 
# trajectories: 
for i in xrange(len(nt)): 
    x[i][0] = bx + 1.0 
    y[i][0] = by 
    z[i][0] = bz 
    L63.time_loop( 
        x[i],y[i],z[i],t[i],Ra,Pr,b,h[i],tol,maxits,mms_switch,bx,by,bz) 
 
# error L2 norms: 
err_x = sp.zeros(len(nt), dtype=float) 
err_y = sp.zeros(len(nt), dtype=float) 
err_z = sp.zeros(len(nt), dtype=float) 
for i in xrange(len(nt)): 
    err_x[i] = sp.sqrt(sp.sum((x[i] - ms[i][:,0])**2)/nt[i]) 
    err_y[i] = sp.sqrt(sp.sum((y[i] - ms[i][:,1])**2)/nt[i]) 
    err_z[i] = sp.sqrt(sp.sum((z[i] - ms[i][:,2])**2)/nt[i]) 
 
# find the observed order of convergence: 
px = sp.polyfit(sp.log(h), sp.log(err_x), 1) 
py = sp.polyfit(sp.log(h), sp.log(err_y), 1) 
pz = sp.polyfit(sp.log(h), sp.log(err_z), 1)

Figure 2 shows the behavior of the solution with MMS forcing when the manufactured solution is not in the convergent region. The solution really is not converging.

---

(a) X-component with MMS forcing (b) Error convergence Figure 2: Convergence behavior with MMS forcing, bx = 0.0, by = -1.0, bz = 10.0

---

Contrast this with the behavior in Figure 3, which actually shows at least first order convergence (and the solution looks like the one we chose).

---

(a) X-component with MMS forcing (b) Error convergence Figure 3: Convergence behavior with MMS forcing, bx = 10.0, by = -10.0, bz = 40.0

---

The fact that the slope of the line in Figure 3 is 1 rather than 3 means I still have a lurking ordered-error in my implementation (likely in the error calculations themselves), or I’m not yet in the asymptotic range, though this seems unlikely given the small time-step size and very consistent slope. The good news illustrated by Figure 3 is that we can apply MMS to even chaotic systems. Now where’s that bug…