Interesting Developments in Commercial Space

These interesting articles were collected in my daily AIAA news brief.

White House Reportedly Will Fund Commercial Spaceflight.

The Wall Street Journal (1/25, Pasztor, subscription required) website reports that the White House has chosen to fund private companies to take astronauts into space with a program likely to cost $3.5 billion over five years. According to the article, Congress is expected to challenge this because of safety concerns. Even though the article predicts conflicts because it will shift money from existing programs to commercial companies, it also contends there will not be any major program cancellations.

        Meanwhile, in continuing coverage, Space News (1/23, Klamper, subscription required) reported NASA, according to unnamed sources "will not be getting the $1 billion budget boost civil space advocates had hoped to see," but the fate of the Ares I is still unknown. However, the Orion capsule would not be cancelled, according to the sources. "While Obama's funding proposal deviates from the Augustine panel's push for a spending increase, sources said NASA's 2011 budget request is expected to align with the panel's so-called Flexible Path plan." According to the article, "In hindsight, NASA officials say the agency set Ares 1 and Orion on an unsustainable spending trajectory."

        Space Frontier Foundation Chairman Sees Killing Ares As "Clear" Choice. In an op-ed for the Orlando Sentinel (1/23), Bob Werb, chairman of the board of the Space Frontier Foundation, called for the cancellation of the "boondoggle" Ares I rocket program. It is "just pork dressed up as cost-effective human space transportation; it's not just wasteful, but destructive to future space exploration beyond Earth orbit." Werb believes there is an "overwhelming" case against the rocket, which he saw as a "bailout" for shuttle companies. "The commercial alternatives are based on well-tested, mature systems currently used to launch U.S. military, scientific and commercial satellites. Adapting these rockets to carry people is cheaper, faster and better." Werb saw this as a simple and "clear" decision.

Food fight's on!

I'm a little unclear on how you can call spaceflight whose only customer is the government "commercial" though...

In related news,

Bolden To Reveal NASA Budget In Press Conference Next Week.

Space News (1/26, Klamper, subscription required) reports, "NASA Administrator Charles Bolden will unveil the U.S. space agency's spending priorities for 2011 during a Feb. 1 press conference at NASA headquarters here, according to administration officials." The article notes Bolden is also "expected to discuss long-awaited details of the president's funding proposal in the morning, followed by a press conference hosted by the White House Office of Science and Technology Policy (OSTP) to rollout Obama's research and development priorities - including those that affect NASA goals and funding- for the coming budget year, these sources said." The budget is "expected to realign NASA's human spaceflight activities and investments to foster development of commercial systems." Bolden "tentatively" scheduled to conduct another press conference at the National Press Club the following day.

        According to Spaceflight Now (1/25, Clark), the budget "is expected to include new direction for NASA." A White House spokesperson told Spaceflight Now, "NASA is vital not only to spaceflight, but also for critical scientific and technological advancements. The expertise at NASA is essential to developing innovative new opportunities, industries and jobs. The President's budget will take steps in that direction." However, the "fate of the agency's vexed exploration program is still unclear." The upcoming budget "may provide no direct guidance on the Constellation program."

        Former NASA Official Cautions Over Commercial Human Spaceflight. In an op-ed for Space News (1/25, subscription required), Scott Horowitz, former associate administrator for the exploration systems directorate, wrote about the role he saw for commercial space companies in manned exploration. "To date, commercial participants have made limited progress and experienced many challenges. Contracts have been pulled, failures have occurred and schedules have slipped up to two years. This has been a stark reminder that rocket development programs are challenging." Horowitz believes that only after a commercial company delivers cargo to the International Space Station "will be time to carefully consider commercial human spaceflight." Horowitz, listing commercial ventures that have failed, felt that this is "not a time to discard these proven systems" of Orion and Ares "in favor of immature, high-risk commercial options."

Ares and Orion are proven? When did those flights happen?

Word Games: Rhetoric and Climate Policy

This post is about climate rhetoric (sorry, I can't do math all the time, I'm not Erdos). I was googling around (don't ask) and found an interesting article (if you can bear to wade through the academic jargon) which inspired me to play a word game. Here's my re-jiggering of a few of their opening paragraphs:

Configuring democracy science and dissent into a political incongruity, a contradiction of terms, is rhetorically strategic to dividing the world inextricably between good and evil, us versus them, in a deadly dual for global domination. Not since the Cold War has an American administration articulated an apocalyptic vision backed by such a massive commitment of military technocratic might, huge expenditure of economic resources, and wanton sacrifice of human life academic integrity. By presidential IPCC decree, everyone must decide whether they are allies or enemies of the United States Nations in a global war to eradicate terrorism carbon-based energy. No shades of grey, no differences of perspective, no room for dissent can be abided if freedom human civilization is to endure and democracy world government by committee is to prevail. The boundary must be drawn fast and firm between righteous truth and wicked persuasion. Thus, the domestic dissenter skeptical outsider symbolizes democracy's climate-alarmist policy's foreign threat, its enemy Other, a traitor to the people and their cause. Or so an empowered elite would have the public believe rather than suffer even a modicum of democratic rational self-rule policy making.

Accordingly, one might conclude that unmaking the oxymoron of democratic dissent climate-policy skepticism would be tantamount to striking at the rhetorical Achilles heal of a discourse that suppresses the actual practice of politics science in the very realm of the political. The political is the realm of antagonism endemic to human relations, as Chantal Mouffe and Ernesto Laclau emphasize, a realm marked by a basic condition of struggle, contested opinion, and "undecidability" (Laclau and Mouffe 2001, xi). Politics is the process of articulating and taking contingent decisions in a context of irreducible difference, conflict, and division through "persuasive redescriptions of the world" (Torfing 1999, 302). That is, by "the elaboration of a language providing us with metaphoric redescriptions of our social relations" we might achieve a revised, expanded, and provisional hegemony of interpretation and political motivation short of insisting on consensus "in a context crisscrossed by antagonistic forces" and contrary to enforcing an ideologically constructed reality of "fully constituted essences." (Mouffe 1993, 57; Torfing 1999, 116). Indeed, social division "without any possibility of a final reconciliation" is inherent to "the very possibility of a [pluralist] democratic politics," which in turn requires a lively dynamic between consensus and dissent (Laclau and Mouffe 2001, xvii, xiv; Mouffe 1996, 8). Politics, in short, is an "ensemble of practices, discourses and institutions which seek to establish a certain order and organize human coexistence in conditions that are always potentially conflictual" (Mouffe 2000, 101). The central question of democratic politics is how to tame and diffuse antagonism in human relations, not eliminate it, how to articulate strategically a practical but partial unity in a pluralistic context of conflict and diversity, by transforming sheer enemies into legitimate adversaries, i.e., by achieving what might be called a fluid condition of consubstantial rivalry. Thus, by this account, the "aim of democratic politics is to transform antagonism into agonism" so as to establish an "us/them" relation "compatible with pluralistic democracy" (Mouffe 2000, 103, 101).

Without open debate, governments tend to exaggerate the danger to the nation, target unpopular groups for vilification and repression, enact preexisting political agendas under the cover of national security scientific certainty, and generally spawn a culture of secrecy and suppression that fosters poor decision making with regrettable consequences.

Isn't it fun to imagine conspiracy and oppression? I think it is America's second favourite pastime. Silly word games aside, the basic fallacy of much of the rhetoric on both sides of climate policy debates is that scientific consensus should imply political consensus.

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…

Chaotic Method of Manufactured Solutions: Lorenz 63

I quickly brushed over the method of manufactured solutions (MMS) aspects in the previous post on Lorenz 63, you might wonder why. Well, it turns out that applying the MMS to chaotic systems includes an additional consideration that is not normally of concern. Here is the ’standard’ sort of checklist for choosing manufactured solutions (taken from Verification of Computer Codes in Computational Science and Engineering, emphasis mine):

1. Manufactured solutions should be sufficiently smooth on the problem domain so that the theoretical order-of-accuracy can be matched by the observed order-of-accuracy obtained from the test. A manufactured solution that is not sufficiently smooth may decrease the observed order-of-accuracy (however, see # 7 in this list).
2. The solution should be general enough that it exercises every term in the governing equation. For example, one should not choose temperature T in the unsteady heat equation to be independent of time. If the governing equation contains spatial cross-derivatives, make sure the manufactured solution has a non-zero cross derivative.
3. The solution should have a sufficient number of nontrivial derivatives. For example, if the code that solves the heat conduction equation is second order in space, picking T in the heat equation to be a linear function of time and space will not provide a sufficient test because the discretization error would be zero (to within round-off) even on coarse grids.
4. Solution derivatives should be bounded by a small constant. This ensures that the solution is not a strongly varying function of space or time or both. If this guideline is not met, then one may not be able to demonstrate the required asymptotic order-of-accuracy using practical grid sizes. Usually, the free constants that are part of the manufactured solution can be selected to meet this guideline.
5. The manufactured solution should not prevent the code from running successfully to completion during testing. Robustness issues are not a part of code order verification. For example, if the code (explicitly or implicitly) assumes the solution is positive, make sure the manufactured solution is positive; or if the heat conduction code expects time units of seconds, do not give it a solution whose units are nanoseconds.
6. Manufactured solutions should be composed of simple analytical functions like polynomials, trigonometric, or exponential functions so that the exact solution can be conveniently and accurately computed. An exact solution composed of infinite series or integrals of functions with singularities is not convenient.
7. The solution should be constructed in such a manner that the differential operators in the governing equations make sense. For example, in the heat conduction equation, the flux is required to be differentiable. Therefore, if one desires to test the code for the case of discontinuous thermal conductivity, then the manufactured solution for temperature must be constructed in such a way that the flux is differentiable. The resultant manufactured solution for temperature is nondifferentiable.
8. Avoid manufactured solutions that grow exponentially in time to avoid confusion with numerical instability.

As emphasized in # 4 we often have to use the free constants in our manufactured solution to bound the derivatives in the domain. For chaotic systems we’ll need to use not only the free constants, but the initial conditions as well, to achieve another end: ensure that the chosen solution is in the convergent region at all times. This is a sufficient but not necessary condition for the source term to be in the basin of entrainment. When the source term is in the basin of entrainment then we can expect to be able to verify the ordered convergence of our implementation as we would in the non-chaotic case. Otherwises we can expect the truncation error to eventually lead to an arbitrarily large departure from the manufactured solution. This terminology comes from the nonlinear dynamics and control literature [1] [3]. A worked example for the Lorenz ’63 system showing the boundary of the convergent region is in Chapter 6 of [2].

The convergent region is defined by the region of phase space for which the eigenvalues of the system Jacobian has negative real parts. As in any other application of the MMS, we start by choosing some convenient solution consistent with the desiderata enumerated above. Trigonometric functions are convenient in this case:

(1)

(2)

(3)

Now we need the Jacobian of the system (and then its eigenvalues), since our governing equations are definied in Maxima, that is easy enough:

/* Need the Jacobian of the system to find the ’convergent region’ so 
   we can select a useful IC/amplitude for our manufactured solutions */ 
sys_eqns : [eqn_1, eqn_2, eqn_3] $ 
J_cr : zeromatrix(3,3) $ 
for j : 1 thru 3 do 
( for i : 1 thru 3 do 
  ( 
    J_cr[i,j] : diff(rhs(sys_eqns[i]), unknowns[j]) 
    ) 
  ) $ 
[J_cr_eigvals, J_cr_eigmult] : eigenvalues(J_cr) $ 
J_cr_eigvals : fullratsimp(J_cr_eigvals) $ $

As before these expressions are output in Fortran90 format for compilation. Here are the eigenvalues in the complex plane for b_x = 1.0, b_y = -1.0, b_z = 10.0.

---

Figure 1: Jacobian eigenvalues for manufactured solution, bx = 1.0, by = -1.0, bz = 10.0

---

We can immediately see the problem from Figure 1: our manufactured solution is in the region of phase space where one of the Jacobian’s eigenvalues has a positive real part. We can adjust the initial conditions of our solution to get all of our eigenvalues on the left half-plane.

---

Figure 2: Jacobian eigenvalues for manufactured solution, bx = 10.0, by = -10.0, bz = 40.0

---

Being limited to only certain portions of the phase space for conducting verification is fine, because a single well-designed manufactured solution verifies the correctness of the entire implementation (as long as all of the terms are activated). This is the same reason that it is often emphasized in the V&V literature that there is no need for the manufactured solution to be physically meaningful.

References

[1]   Jackson, E.A., “Controls of dynamic flows with attractors,” Physics Review A, Vol. 44, Is. 8, 1991.

[2]   Wu, W., “Analytical and Numerical Methods Applied to Nonlinear Vessel Dynamics and Code Verification for Chaotic Systems,” Disertation, Virginia Polytechnic Institute and State University, Dec, 2009.

[3]   W. Wu, L.S. McCue, and C.J. Roy. “The method of manufactured solutions applied to chaotic systems,” Nonlinear Dynamics, 2009, under review.

Lorenz 63 Ensemble

Dan Hughes has started an ambitious series of posts about formally verifying (using the method of manufactured solutions (MMS)) a range of chaotic ordinary differential equations (ODE) systems.

One of the classic ones he’s looking at is Lorenz’s three-dimensional toy system from 1963:

(1)

(2)

(3)

where the system exhibits chaotic behaviour depending on the choice of parameter values (Ra = 28, Pr = 10.0, b = 8∕3 will be used throughout).

The first thing to do is define our governing equations in a computer algebra system (I like Maxima) because that makes writing bug free numerical software and complicated MMS forcing functions much easier.

  unknowns : [x,y,z] $ 
  depends(unknowns, t) $ 
 
  /* define the governing equations: */ 
  expr_1 : -Pr * x + Pr * y - diff(x,t) $ 
  eqn_1 : solve(expr_1, diff(x,t))[1] $ 
  expr_2 : -y + Ra * x - x * z - diff(y,t) $ 
  eqn_2 : solve(expr_2, diff(y,t))[1] $ 
  expr_3 : -b * z + x * y - diff(z,t) $ 
  eqn_3 : solve(expr_3, diff(z,t))[1] $

Then we need to define our discrete approximation to the governing equations (I’ll skip the MMS stuff for now)

/* elements of the Butcher tableau: */ 
A : matrix([1/4, -1/4], 
           [1/4, 5/12]) $ 
bT : matrix([1/4, 3/4]) $ 
c : matrix([0], [2/3]) $ 
 
/* identity matrix for doing Kronecker products, dimension depends on 
   the size of the system you are trying to solve */ 
sys_I : diagmatrix(3, 1) $ 
 
K : transpose([k[1], k[2], k[3], k[4], k[5], k[6]]) $ 
Xn : transpose([x[n],y[n],z[n],x[n],y[n],z[n]]) $ 
tn : transpose([t[n], t[n], t[n], t[n], t[n], t[n]]) $ 
F : transpose([f[1], f[2], f[3], f[4], f[5], f[6]]) $ 
 
A_rk : kronecker_product(A, sys_I) $ 
bT_rk : kronecker_product(bT, sys_I) $ 
hc_rk : kronecker_product(t[n] + h*c, sys_I) $ 
 
/* argument to our system right-hand-side operator: */ 
nonlin_arg : Xn + h * A_rk . K $ 
/* set up the two-stage system: */ 
rk_sys : zeromatrix(6,1) $ 
q_sys : zeromatrix(6,1) $ 
rk_sys[1,1] : k[1] - subst(x=nonlin_arg[1,1], subst(y=nonlin_arg[2,1], 
    subst(z=nonlin_arg[3,1], rhs(eqn_1)))) $ 
rk_sys[2,1] : k[2] - subst(x=nonlin_arg[1,1], subst(y=nonlin_arg[2,1], 
    subst(z=nonlin_arg[3,1], rhs(eqn_2)))) $ 
rk_sys[3,1] : k[3] - subst(x=nonlin_arg[1,1], subst(y=nonlin_arg[2,1], 
    subst(z=nonlin_arg[3,1], rhs(eqn_3)))) $ 
rk_sys[4,1] : k[4] - subst(x=nonlin_arg[4,1], subst(y=nonlin_arg[5,1], 
    subst(z=nonlin_arg[6,1], rhs(eqn_1)))) $ 
rk_sys[5,1] : k[5] - subst(x=nonlin_arg[4,1], subst(y=nonlin_arg[5,1], 
    subst(z=nonlin_arg[6,1], rhs(eqn_2)))) $ 
rk_sys[6,1] : k[6] - subst(x=nonlin_arg[4,1], subst(y=nonlin_arg[5,1], 
    subst(z=nonlin_arg[6,1], rhs(eqn_3)))) $

This is basically the same approached that I applied to the simple harmonic oscillator, but this time our system is non-linear, so we can’t invert the time-integration operator directly. We need to use an iterative Newton’s method to solve for each time-step, for which we need to define a Jacobian:

/* calculate the system Jacobian for use in a Newton’s method: */ 
J : zeromatrix(6,6) $ 
for j : 1 thru 6 do 
( for i : 1 thru 6 do 
  ( 
    J[i,j] : diff(rk_sys[i,1], k[j]) 
    ) 
  ) $ 
 
time_update : factor(h* bT_rk . K) $ $

The really nice thing about putting in a bit of work up-front with the symbol manipulation software is you can just dump the resulting expressions to compile-able code:

load(f90) $ 
file_output_append : false $ 
fname : ”rhs.f90” $ 
with_stdout(fname, print(”! generated by lorentz63_mms.mac”)) $ 
file_output_append : true $ 
with_stdout(fname, f90(f[1] = rk_sys[1,1])) $ 
/* ...much more boring file IO follows... */

With all of our useful expressions safely (and correctly) in Fortran, we can compile with F2py and play with them in Python. Integrating the system with a third-order implicit Runge-Kutta method gives the following time-step convergence behaviour:

---

(a) X component (b) Y component (c) Z component Figure 1: Time-step convergence for Lorenz ’63 system with 2-stage Radau-IA RK scheme

---

Figure 1 shows that as the time-step is decreased the trajectory stays with the ’pack’ for longer and longer times, but when the answers do diverge significantly the change is quite abrupt.

The most practically interesting thing with chaotic systems is the high sensitivity to initial conditions. This sensitivity competes with our uncertain knowledge of the true state of things in determining the limit on our ability to predict the trajectory of (some) dynamic systems (even if our knowledge were perfect, the finite precision of our machines still causes us significant trouble). Part of the way we deal with uncertainty is to do averaging over ensembles. To illustrate this idea I’ll start 7 realizations of the previous calculation, each adjusted by a small amount (on the order of the floating point round-off error for my machine).

nt = 400000 
T = 40.0 # integration period 
Pr = 10.0 
Ra = 28.0 
b = 8.0 / 3.0 
# knobs for the Newton’s method: 
tol = 1e-12 
maxits = 20 
nic = 7 # number of different IC’s in our ensemble 
h = T / float(nt) 
x = sp.zeros((nt,nic), dtype=float, order=’F’) 
y = sp.zeros((nt,nic), dtype=float, order=’F’) 
z = sp.zeros((nt,nic), dtype=float, order=’F’) 
t = sp.linspace(0.0, T, nt) 
# some perturbations: 
eps = 1e-16 * sp.array([0, 1, -1, 2, -2, 3, -3]) 
for i in xrange(nic): 
    x[0][i] = 1.0 + eps[i] 
    y[0][i] = -1.0 + eps[i] 
    z[0][i] = 10.0 + eps[i]

I’ll then calculate the pair-wise difference between each trajectory (thats 7 choose 2 or 21 differences). The combination function from the itertools module comes in handy here

from itertools import combinations 
from scipy.misc import comb 
 
# calculate our ensemble of trajectories: 
for i in xrange(nic): 
    L63.time_loop(x[:,i],y[:,i],z[:,i],t,Ra,Pr,b,h,tol,maxits) 
 
# calculate all the pair-wise differences between trajectories (nic 
# choose 2 of em): 
npairs = comb(nic, 2, exact=1) 
xdiffs = sp.zeros((nt,npairs), dtype=float, order=’F’) 
ydiffs = sp.zeros((nt,npairs), dtype=float, order=’F’) 
zdiffs = sp.zeros((nt,npairs), dtype=float, order=’F’) 
for i, pair in enumerate(combinations(xrange(nic), 2)): 
    xdiffs[:,i] = x[:,pair[0]] - x[:,pair[1]] 
    ydiffs[:,i] = y[:,pair[0]] - y[:,pair[1]] 
    zdiffs[:,i] = z[:,pair[0]] - z[:,pair[1]]

Here’s the resulting differences:

---

(a) X component (b) Y component (c) Z component Figure 2: Differences between members of the ensemble

---

You can see that things go along well in the early time, but eventually the differences between ensemble members grows suddenly to be on the same order of magnitude as the things we care about predicting.