Complex Systems

comp90072 – The Art of Scientific Computing

1 Introduction

1  Introduction

One might ‘intuitively’ argue that in general, catastrophic events have equally catastrophic causes; or more formally, large perturbations to the state of a system must have proportionally large causes. For example, a planet orbiting its star will deviate from its orbit in proportion to the size of the object that hits it such that a small meteorite crashing into the Earth has a negligible effect whilst having a black hole fly through the Solar System might result in the Earth being flung out into interstellar space. Although there are countless examples in physics were the cause and effect are proportional, this is not always the case as we will investigate in this lab.

Although the scenario of planet orbiting a star is quite easy to describe, not all systems are quite so simple. Many systems (and often the more interesting ones) in fact contain an extremely large number of interconnected components giving rise to an enormous number of degrees of freedom. Naïvely, one might assume that if we can describe how each component behaves with its immediate neighbours, we might be able to draw conclusions as to the behaviour of the system as a whole; but this is in general not possible save for the very simplest of microscopic interactions. Even if we are able to write down the Lagrangian of the system, solving the equations of motion will be intractable.

In this lab, we will be investigating an (initially simplistic) example of a sand pile on a lattice. The general idea is that we will be adding grains of sand onto a lattice and we will have to ‘topple’ points of the lattice according to a set of rules.  Despite the rules being extremely simple to begin with, we will see that the complex behaviours arise.

1.1  Self-organized Criticality and Scale Invariance

In 1987, Bak, Tang and Wiesenfeld introduced the concept of self-organized criticality (soc) [1,2]. It is a property of certain systems (typically slowly driven complex systems) whereby it will evolve towards a particular configuration without the need to finely tune parameters of the system. This particular configuration the system is evolving towards is called the attractor as it ‘attracts’ the system to it.

In addition to having an attractor, systems exhibiting socgenerally exhibit some sort of scale invariance. As the name suggests, scale invariance is the property whereby a system looks the same over a broad range of scales.  Most likely you have seen this property in mathematics in fractals whereby the pattern repeats itself as you zoom in as seen in figs. 1 and 2. Scale invariance can also arise in real physical systems and manifests itself by forming some sort of power-law for the frequency of an event compared to its size. Examples of phenomena which exhibit power-law behaviour include:

•  earthquakes (which are described by the Gutenberg–Richter Law);

• solar flares;

Figure 1: Mandelbrot set which exhibits scale invariance as certain patterns repeat at different scales (above pictures showing different ’zoom’ into the Mandelbrot set).

• forest fires;

•  epidemics; and,

•  fluctuations in financial markets.

Note that all of the systems above are slowly driven. That is, they are only ever affected by small localized effects and thus are never greatly perturbed.

Question 1:  Why is it that a system exhibiting scale invariance will have phenom- ena which are described by a power law?

Self-organizing criticality as a framework has been used to study an extraordinarily diverse range of systems. It has been used to characterise, analyze and predict earthquakes, financial markets, evolution and extinction events, pulsar glitches, neuronal behaviour and much more.

Figure 2: The btw sand pile exhibiting scale-invariant patterns. In this sand pile, the sand was dropped continuously onto the centre.

1.2  Non-equilibrium Statistics

As discussed above, the systems we will be investigating have an enormous number of degrees of freedom making it difficult to work with analytically and impossible to solve the exact dynamics of the system. In order to still get an insight on the system, we must rely on statistical tools and in particular, tools which work with non-equilibrium systems.

Non-equilibrium statistics underlies how we think about systems that display self-organizing criticality and does so in a holistically. When a system is critical, microscopic phenomena    (such as a local increases in stress) are not responsible for the macroscopic observables of the    system (such as an earthquake’s size, duration, area etc). In particular, the proportions of    large and small events do not depend on the exact microscopic details of what is happen-    ing. Consequently, one cannot analyze the building blocks of a system in order to explain    the large scale phenomena. In other words, a detailed understanding of the microscopic    interactions of the system do not translate into system-level understanding.

1.3 Aim

The primary aim of this lab will be to build a basic sand pile model and study socwithin this system and the associated non-Gaussian statistics. The basic model can subsequently be extended in any of a number of ways to replicate more physically accurate systems or investigate how the system behaves under different conditions.

2 Bak–Tang–Wiesenfeld Model

2  Bak–Tang–Wiesenfeld Model

When modelling a complex system, it is often difficult (or impossible) to create mathematical formalism that is both sufficiently realistic and computationally/theoretically tractable. As a result, researchers will consider a simplified model that captures the important elements of the physical system in question. The btwmodel is an excellent example of a simplified model which exhibits soc.

The btwmodel describes a (somewhat abstract) sand pile. The grains of sand are located on an 푁 × 푀 lattice such that 푝 t(푖, 푗) is the number of grains located on site (푖, 푗) at time 푡, where 푖 ∈  [1, . . . , 푁], 푗  ∈  [1, . . . , 푀] and 푡 ∈ N0. At each time step 푡, a new grain of sand is added to a random site in the lattice causing the sand piles at each site to grow.

If an individual pile of sand reaches a height of 4 or more, then the pile topples and distributes four grains of sand to its four immediate neighbours (we are not considering diagonal neighbours):

푝 t(푖, 푗)  = 푝 t(푖, 푗) - 4,

푝 t(푖 ± 1, 푗) = 푝 t(푖 ± 1, 푗) + 1, 푝 t(푖, 푗 ± 1) = 푝 t(푖, 푗 ± 1) + 1.

Note that this process is instantaneous and does not increment 푡 . Furthermore, if the toppling happens on the edge of the lattice (i.e. 푖 ∈ {1, 푁 } and/or 푗 ∈ {1, 푀 }),the grains that would topple to a site off the table are deleted from the system.

Question 2: In what ways does the btwmodel satisfy the properties required for a system to exhibit soc?

Question 3:  Why is it important for the table to have finite boundaries?

What if we had an infinite lattice, but dropped the sand on only a (finite) subset of the lattice?

If the lattice begins with no sand at all, we might expect the sand to just accumulate to begin with as very few grains of sand will topple off the edge. Eventually though, the slow addition of sand should balance out with those grains of sand toppling off the edge of the grid.

Figure 3: An example of a simple avalanche taking place on a 3 × 3 grid.

Question 5:  Explain the difference between: a steady state, a statistically station- ary state, and an equilibrium state.

Do you expect the btw model reach an equilibrium state or statistically stationary state? In terms of system observables, how might we characterize this state?

The addition of a grain of sand to the lattice has a chance of creating an unstable site which topples; but the toppling in turn can create further unstable sites which will subsequently topple. As a result, the addition of a single grain of sand can create a chain of toppling which we call an avalanche as depicted in fig. 3. It will be these avalanches that we will be most interested in. In particular, we will consider the following properties:

Topples  The number of topples required to take the grid from its initial unstable config- uration to a stable one;

Area  The number of unique sites which have been toppled;

Loss  The number of grains of sand which are toppled off the lattice; and, Length  A ‘length’ which characterizes the avalanche.

Note that we are purposefully leaving the length definition quite vague. You’ll have to decide where from and to this length is defined, and which metric to use (such as the 퐿1- norm or 퐿2-norm, a.k.a. Manhattan distance and Euclidean distance respectively). Don’t forget to justify your choice.

Question 6:   Is there a difference between sampling a system over time and sampling an ensemble of systems? If so, under what circumstances?

Hint: Lookup the property of ergodicity and argue whether or not you believe the btw model is ergodic.

Question 7:  Record the frequency with which avalanches have a specific number of topples, area, loss and length. What kind of distribution do you get for each?     If possible, can you find a fit for the distribution? (E.g. for a Gaussian, find the mean and width of the distribution; for a power law, find the exponent.)

Question 8:  Do you expect any correlation between the four properties we are recording for the avalanche? Justify why/why not, then verify this with your simulation.

2.1  Abelian Property

Although the btwmodel is a model of‘sand piles’, it is extremely unrealistic—after all,sand piles aren’t formed on a lattice and even if they were, a site wouldn’t go from four grains of sand to none from toppling on its neighbours. Nevertheless, it exhibits an extremely useful mathematical property when it comes to analysing its behaviour: the topplings are Abelian. In other words, it doesn’t matter in which order we execute each toppling, the end result is the same.

To give a concrete example, consider the following start to an avalanche:

There are now two unstable sites, (1, 2) and (3, 2).  If the model were not Abelian, we would need to keep track of the order in which sites are toppled and make sure that we topple the correct next unstable site; however as the btwsand pile model is Abelian, we can do this in whatever order we wish and arrive at the same result thereby greatly simplifying the computational complexity of avalanches.

You might take it for granted that the btwmodel is Abelian, but this can be fairly easily shown rigorously. To begin with, we can consider the lattice of sand piles to be 퐿 which is a finite subset of Z푑 where 푑 is the dimensionality of the lattice. For any site 풙 ∈ 퐿, the number of grains of at each point of the lattice is given by the function

푝 : 퐿 → N0 .                                       (3)

Astable configuration corresponds to the scenario where 푝(풙 ) < 푘 for all 풙 ∈ 퐿 while an unstable one has at least one site 풙  ∈  퐿 where 푝(풙 )  ≥  푘 where 푘 defines is the stability threshold.

From an unstable configuration 푝, we must introduce a way to describe the toppling; for this, we introduce the toppling matrix Δ (풙 , 풚). This stabilizes an unstable site 풙 in the configuration 푝 by redistributing the sand to site 풚 . The toppling matrix must satisfy the following conditions:

∀풙 , 풚 ∈ 퐿 where 풙 ≠ 풚 ,    Δ (풙 , 풚) = Δ (풚, 풙 ) ≥ 0,                (4a)

∀푥 ∈ 퐿,    Δ (풙 , 풙 ) < 0,                            (4b)

∀푥 ∈ 퐿,                          (4c)

풙                              (4d)

Question 9:  For each condition in eq. (4), explain the significance and justify its necessity.

The actual definition of the toppling matrix in the btwmodel is:

2푑  re nearest neighbours .                   (5)

>:0        Otherwise

Question 10:  Explain and justify the definition of the toppling matrix in the btwmodel and verify that it satisfies all the conditions set out in eq. (4). What is the stability threshold 푘 ?

Given the definition of the toppling matrix, it is now possible to define the toppling operator 푇풙 which maps a configuration 푝 to a new configuration 푝 /:

Question 11:  Show that 푇풙 and 푇풙/ commute for unstable configurations.

In general, multiple topplings are required in order to stabilize a configuration 푝 and we can introduce a stabilization operator which maps an unstable configuration to a stable configuration

T : 푈L  → 푆 L                                                                     (7)

where 푈L is the space of all possible height configurations and 푆 L  C 푈L is the subset of stable height configurations. To give a concrete example, the stabilization operator for fig. 3 consists of four topplings:

푝 / = T푝 = 푇(2,3)푇(1,3)푇(1,2)푇(2,2)푝 .                          (8)

3 Extending the btw Model

Thus we see that a stabilization operator can be represented as

푁

where 푁 is the number of instabilities through an avalanche.

Question 12:  Show that Tis well defined. In other words, show that given some unstable configuration 푝, 푝/ = T푝 is unique (you can’tend with two different stable configurations). This is equivalent to the sequence of 푇풙 푖  required is unique up to re-ordering.

Hint: Read the original papers for hints as to how this can be done.

3  Extending the btwModel

The btwmodel is clearly unrealistic; however, you have just seen that it still exhibits very interesting statistics. You can now extend this model in any of a number of ways. Irrespective of which extension(s) you choose, remember to make predictions as to what you expect and then test these predictions with your simulation. You should also compare your extension to the baseline btwmodel: in what ways is it similar and in what ways does it differ?

A few extensions you may consider (and combine or extend as your please):

1.  Change the toppling so that it depends on the gradient: if a site has 푘 more grains of sand than its neighbour it is unstable and it drops 1 grain of sand onto this neighbour and repeat until it is stabilized.

2.  Change the toppling to be random: instead of redistributing the four grains of sand to all four neighbours, it selects one of the four neighbours at random for each grain of sand. What about adjusting the probability to drop to a neighbour based on the height difference?

3.  Change the boundary conditions:

a)  Make the toppling ‘wraparound’ in 푥 (like a cylinder) so that grains of sand can only drop off the top or bottom of the grid.

b)  What happens if grains of sand can’t topple off the edge, but instead there is a hole in the centre of the grid (imagine a kind of hourglass).

4.  Change how fast grains of sands are added to the grid. Instead of adding 1 grain of sand at a time, add 2, or 3, or 푛 . Do you add each of these grains to different sites or the same?

5.  Change how you are adding grains of sand. Instead of being uniformly random on the grid, what if grains of sand were only dropped on a small subset of the grid, or followed a different distribution (such as a Gaussian centred on the middle of the grid).

6.  Make the avalanche happen really slowly so that the toppling rule, eq. (1), becomes:

p 푡+1(i, j) = p 푡 (i, j) - 4,

p 푡+1(i ± 1, j) = p 푡 (i ± 1, j) + 1,                        (10)

p 푡+1(i, j ± 1) = p 푡 (i, j ± 1) + 1.

Perhaps you don’t want to increase the timestep every topple, but every nth topple.

7.  Instead of a rectangular lattice, change it to a hexagonal lattice (have a look at hexagonal coordinatesystems).

8.  Change the toppling rule so that instead of dropping to the 4 directly adjacent sites, it drops to the 8 surrounding sites. Or maybe change it so that it drops to all sites with a radius of 2 of the initial site.

9.  Add a force that will push grains of sand in a particular direction (such as wind, or if the grid was sloped).

References

References

[1]   P. Bak, C. Tang and K. Wiesenfeld, ‘Self-organized criticality: an explanation of the 1,푓 noise’, Phys. Rev. Lett. 59, 381–384 (1987)1 0 .  1 1 0 3 / P h y s R e v L e t t  . 5 9  . 3 8 1 .

[2]   P. Bak, C. Tang and K. Wiesenfeld, ‘Self-organized criticality’, Phys. Rev. A 38, 364– 374 (1988)1 0  .  1 1 0 3 / P h y s R e v A . 3 8  . 3 6 4 .