Monte Carlo analysis of the Lorenz equations – Part II

Cont'd from here.

The non-linearity in the model is clear from this scatter plot. For small perturbations in the inputs from the mean values (depicted as the square dots below), the solutions can end up being anywhere in the space depicted by the scatter plots.

Stated in Lorenz’s words, a flap of the butterfly’s wings for the input conditions can create a drastic change in outputs, possibly resulting in a tornado. In Meteorological terms, this is one of the main reasons why longer term weather forecasts remain a difficulty. Even a statistical analysis such as the one performed here is worthless in predicting the longer term effects of weather, since conditions such as temperatures can lie anywhere in the contour of the scatter plots and thus mean nothing in terms of accurate predictions. This is the crux of Lorenz’s work; that accurate longer term predictions cannot be made when it comes to weather. Further, he leaves the question of the butterfly creating a tornado in a different hemisphere of the world unanswered.

Predictability of phenomena is more reliable for linear systems. While short term predictions can be made for non-linear systems, exact longer term predictions are both unreasonable and impossible. All that one can predict precisely about the location of the electronic butterfly in response to a tap on the jar is that it can be found some place inside the bottle after a few seconds. In a similar (and a rather cynical) vein, all that the Meteorologist can predict with certainty for longer term forecasts is that it will be summer in July in the Northern hemisphere. It only takes a mental extrapolation to extend this analogy to even longer term predictions with climate change, and the reasons for the discussions that ensue!  

Monte Carlo analysis of the Lorenz equations – Part I

Cont'd from here.

Solving a differential equation on a computer in the early 70’s was by no means an easy task, let alone solving a set of 12 aperiodic coupled non-linear ODE’s. Lorenz states that he got a break when he met with Dr. Barry Saltzman, whose work on thermal convection resulted in a set of 7-equations, 4 of which approached zero. Lorenz was able to modify the other three equations and solve for the aperiodicity of the outputs. These three equations are commonly referred to as the Lorenz model.

How does one study the aperiodicity of the Lorenz model? By following Lorenz’s lead and introducing minor perturbations to the initial values. And that is exactly where the Monte Carlo technique comes in useful. As stated earlier, the MC technique works on the premise of allocating randomized probabilities for perturbed inputs and studying the behavior of the system. The inputs used for the process were as follows:

Randomized initial values were thrown as inputs to the model within the bounds shown above. After 2000 such throws, the outputs were obtained for a specific instant of time, 15 time units in this case. The results of the scatter phase-plots are shown below.To be cont'd.

 

Edward Lorenz and the Butterfly Theory

Predictability; Does the flap of a butterfly’s wings in Brazil set off a tornado in Texas? This was the question proposed by Dr. Edward Lorenz at the 135^th Meeting of the American Association for the Advancement of Science in 1972. Dr. Lorenz can be stated as one of the pioneers of modern day Meteorology, and laid the foundations for equations in weather prediction. He went on to research on whether the behavior of atmosphere is “unstable” with respect to small perturbations in amplitudes when it comes to long term weather predictions.

In his other seminal work, “On the Prevalence ofAperiodicity in Simple Systems”, Lorenz states that weather predictions can be likened to a set of non-linear differential equations, whose solutions are heavily governed by the initial values. He states that periodic solutions would be naïve in weather forecasting, since if there were indeed a periodicity, the statistical prediction would be a trivial one.

It turns out that Lorenz started out with solving a set of 12 differential equations and stumbled upon the aperiodicity in the solutions with minor changes to input conditions. He noticed that re-solving the equations with small differences in the initial values resulted in vastly different outputs over a period of time. To quote Lorenz, “These perturbations were amplified quasi-exponentially, doubling in about four simulated days, so that after two months the solutions were going their separate ways”. Meaning, the instability in the solution was “due to the lack of its periodicity”. Lorenz’s work on non-periodic equations for weather predictions has been the founding blocks for Chaos Theory and the “Butterfly Effect” . 

To be cont'd. 

Implicit ODE’s, Rosenbrock Formulation & the Josephson’s Junction:

ODE’s of the form

are referred to as implicit ODE’s. While a matrix inversion can be used to convert these to an explicit ODE, the form of the Rosenbrock integration scheme can be modified to directly solve for an implicit 1^st order set of ODE’s.

The complete integration can be performed using the same procedures described earlier.

Getting back to numerically solving the Josephson’s Junction equations, which is an implicit ODE, the above formulation can also be used as an alternative and avoid the matrix inversion mentioned earlier.

2nd order ODE of Josephson’s Junction

The system of two coupled non-linear 2^nd order differential equations that represent the Josephson’s Junction phenomenon on superconductors is presented here, as adapted from Hairer et al . 

 

The equations describe the electron-pair interactions between two superconducting elements that are separated by a thin insulator, under the influence of an applied current source. The obvious first step in attempting to numerically integrate these equations is to decouple it.

Coupled equations can be decoupled in a variety of ways. Perhaps the easiest one is the brute-force method, where in a few algebraic sleights are performed to eliminate one the variables in either equation. In the subject case, multiplying the second equation by "alpha" and adding it to the first to eliminate one of the 2^nd derivative variables. Subsequently, this equation can be substituted into the first equation to obtain the decoupled versions. 

The algebraic sleight works best with a limited set of equations (usually up to two at the most). As the number of equations increases, the algebraic decoupling process becomes rather cumbersome. An alternative is to use the matrix formulation instead. The Josephson equations then become:

 A matrix inversion can then be performed to obtain the decoupled equations: 

Since the matrix that requires the inversion is a constant, only a single inversion is required in the integration process.

And the two 2^nd order ODE’s can be considered as being successfully decoupled. Numerical integration is just a step away, once the 2^nd order ODE’s are converted to 1^st order ODE’s.

Plotted below is the numerical solution of the equations, performed using a Rosenbrock technique. The initial conditions and values for the constants from Hairer et al were used for the analysis. 

Some observations on the Monte Carlo simulations

Setting aside the triviality of the premise for “calculating” pi (i.e., one needs to be cognizant of the fact that the area of the circle is represented using a pi term in it, implying that pi is already known to the user), this simple example illustrates the basics of the Monte Carlo technique. The following general observations can be made from the pi-experiment:

•          A mathematical model needs to be in place in order to perform a sensitivity analysis of the system.
•          Given a model, a prediction can be made for the steady state output of the system, for a set of input variables.
•         Should a question arise on the validity of the outputs, the MC technique can be consulted. Specifically, the question can arise in the form of what variation can be statistically forecast in the outputs in the event of a variation in the inputs.
•          Having posed the question thusly, the parametric statistical variation can be performed to the inputs over a bounded interval, and the behavior of the mathematical model can be obtained.

If the system being analyzed is linear, a bell-curve distribution can be expected for large sample sizes randomized multiple times. One needs to be cognizant of the fact that MC simulations are highly statistical, and depend heavily on the validity of both the input variables and the mathematical model. In the case of the pi-estimation experiment, a sample size of 1000 can only have a discrete number of marbles fall within the confines of the circle. Thus, the accuracy of the estimate for pi is limited to the 1000^ths place; i.e., only three decimal places can be predicted for pi with a 1000-sample bin. For a 10,000 size bin, only a 4-digit accuracy can be predicted. However, repeated throws and subsequent averaging using the MC simulation can statistically result in values of pi that include far more number of digits. This is clear from the histograms shown earlier.

That greater accuracies can be seemingly predicted with a MC technique that is initially not available in the model is a rather disturbing thought. Adequate understanding of the Physics of the system being modeled, combined with proper bin sizes for the data analysis in a real-world system is thus an absolute necessity before one can launch off on an MC simulation to predict meaningful results.