Quantitative Analysis of Physical and Statistical Properties of Flocks.

(1)

Theoretical Physics

Quantitative Analysis of Physical and Statistical

Properties of Flocks

Joakim Bergdahl

jbergda@kth.se

Lars Almgren

larsalmg@kth.se

SA104X Degree Project in Engineering Physics, First Level

Department of Theoretical Physics

Royal Institute of Technology (KTH)

Supervisor: Jack Lidmar

(2)

Abstract

Flocking behavior is a common phenomenon in nature in the form of, for instance, flocks of birds or schools of fish. Making the assumption that the members of a flock can be considered a system of interacting particles it is possible to use methods from statistical physics to quantitatively analyze flocks and their properties. This way, a flock can exist in various thermodynamic phases and exhibit phase transitions depending on changes within the flock. In this report the flock analysis is performed with the help of a model originally created by Vicsek et al. The model is governed by certain parameters controlling the interaction between individual flock members. From the results it is possible to see that even small deviations in these parameter values can lead to great phase alterations as well as phase transitions, which strengthens to the assumption that a flock can be considered a system of particles.

(3)

Chapter 1 Introduction

This chapter will provide an introduction to the report. The chapter begins with a brief review of the history of flock modeling. Later in the introduction, the structure of the report will be introduced.

1.1 History of Flock Modeling

The phenomenon of flocks is something that the general public is considerably familiar with. It is certain to say that almost everyone has seen instances of flocking behavior in flocks of birds emigrating in large numbers to warmer latitudes or colonies of ants crossing forest pathways. Something a little less familiar is the evident existence of flocking behavior on other levels than macroscopic ones. For instance, cells and even microorganisms display group behavior that, in a sense, could be compared to that of birds [1]. As a result of this versatile nature of flocking behavior, the available areas of research are endless. Previous work has led to a variety of models that tries to explain how a large group of organisms behave [2]. How is it possible that these collections of organisms coherently move together? Given such an interesting field of research; how can the nature of such a collection be quantified using sound methodology and accurate means?

When observing a flock you can easily tell that there exists some kind of internal dynamics between its members. Using observable behavioral properties it is possible to construct a model for simulating the internal interactions and natural motion of a flock [1]. As early as 1952, John T. Emlen observed and performed an empirical study on the behavior of cliff swallows and noticed that they followed a distinct interaction pattern [3]. This report will focus on bird flocks in flight where the birds initially will be considered a collection of moving members in a general two-dimensional plane. This is a reasonable assumption for low flying birds whom are constrained by the ground, or more precisely by the wish to avoid colliding with the ground. Craig Reynolds, one of the pioneers in

(6)

the field of computational flocking simulation, proposed a simulation model using just a few key ingredients in 1987. His idea was that each moving flock member, which he referred to as a boid, followed three simple rules. Firstly, each boid had a desire to avoid its flock neighbors in order to limit cluster generation. Secondly, it tried to align itself with its neighbors. Lastly, it had a desire to move towards the center of mass of its flock neighbors [4].

1.2 Approach using Statistical Physics

In a later stage in the field of flock modeling, Vicsek et al. proposed an approximate model where each flock member was represented as a self-propelled particle with a cer-tain interaction potential and velocity. This approach made it possible to analyze the collective behavior using methods originating from statistical physics. A flock could, on a microscopic level, display patterns related to states and phases similar to those of actual solid-state matter, gases and liquids [1]. Controlling the interaction made it possible to force state and phase transitions only by adjusting the model parameters. Analogous to a ferromagnetic equilibrium state, Vicsek found that the velocities of the flocking members self-aligned like the magnetic spins of iron atoms in a ferromagnet without the need of an external influence. This meant that an aligned flock moving in a distinct direction could be achieved solely by the inter-bird interaction without the birds themselves following a preferred direction [2].

Toner et al. suggested a generalization of Vicsek’s particle model based on the results made by biologists Andreas Huth and Christian Wissel [5]. It was proposed that the original Vicsek model (VM) lacked the specific interaction needed for the flock members to keep a preferred distance between each other [2] making it slightly less optimal for simulation purposes.

1.3 Synopsis

As a point of reference, a further revision of this modified Vicsek model [2] will be used to simulate the flocking behavior in this study, which is further discussed under section 2.2. The revision originates from the work performed by Tu et. al. and will, for the sake of simplicity, be referred to as the Tu-Vicsek model (TVM) in this study [6]. The results by Tu et al. was based on very large flocking systems. One of the purposes in this study is to investigate if comparable results can be obtained for a smaller system. The produced simulation data is used for the purpose of quantification of the statistical properties of the occurring flock phases [1] with the analytical means presented under section 2.3. Concluding the study, the results are featured and further discussed in chapter 3.

(7)

further research in the area of flocking modeling. The main target group for this thesis is under graduates and academics with an interest in the field of statistical physics. Using the material presented in chapter 3, it is possible to get a better understanding of the flocking phenomenon and develop viable methods for further research regarding, for instance, computer simulation.

(8)

Chapter 2 Investigation

In this chapter the investigation of the study will be presented. It starts with the def-inition of the problem. In the model section, the model in the report and some other models will be introduced. The analysis section explains the tools used for analyzing a flock. The results of the study will then be presented in the results section. The results are then explained and discussed in the final section.

2.1 Problem

To be able to describe the motion of a bird within a flock there is a key set of parameters to consider in order to describe its motion as naturally as possible. To simulate a moving flock a number of rules that regulates the inter-flock dynamics has to be postulated. These rules can be used to, for instance, form the basis to a realistic simulation model. The main problem that is considered in this thesis is the one of the physical quantification of simulation data produced by such a model. Vicsek et al. suggested that using particles with corresponding interaction potentials to emulate flocking behavior gave way to the application of successful methods originating in statistical physics in order to classify collective motion [1]. This is only one approach to the modeling problem and there are others to consider before the most beneficial model can be selected. For instance, Hildenbrandt et al., Cucker-Smale and Hayakawa each suggested models describing the phenomenon [7–9].

When a model has been selected, it can be used to produce simulation data. A series of statistical measurements can then be applied to this data in order to quantify the flocking behavior resulting from the model. Among the statistical quantities that should be observed, the one of order is the most apparent. In this study, order will be mea-sured using the magnitude of the polarization of the flock, which describes how a flock transforms over time due to flock member alignment. Furthermore, the variance of this parameter will be observed for flocks in stationary states in order to measure the

(9)

long-lasting influence of statistical noise. The density and mean distance within a flock will change during the course of the simulation, making them important tools for tracking the flocks’ reactions from changes in the model parameters. Using all of these tools, the results of a change in model parameters can be investigated and critical values can be obtained.

One of the purposes of this report is to investigate if results from earlier work, where large flocks with about 320,000 birds have been used [2], can be reproduced for a smaller system with a lower number of birds.

2.2 Model

In this section, some of the existing models for simulating flocking behavior will be discussed. It is also explained why the TVM is chosen for this study, as well as how it is constructed, what set of rules it is based upon and how it works.

2.2.1 Different Models

The field of flocking behavior is vast and there exists a variety of models describing the phenomenon. Here, three models will be examined closer, constructed by Hildenbrandt et al., Cucker-Smale and Hayakawa [1, 7–9].

Hildenbrandt et al.

In the model Hildenbrandt et al. presented, the behavior of an individual was based on its cruise speed, its social environment, its attraction to the roost, and the simplified aerodynamics of flight that includes banking while turning. The model is three dimen-sional and built using SI units, where real parameters are used when available. This approach gives a more realistic environment for the problem. As in Hayakawa’s model, Hildenbrandt et al. also uses individuals with mass, which makes it possible to use forces that act on the individuals, such as a force that act as a speed control. If the individual deviate from its cruise speed v0, the force will compel it to accelerate or decelerate. Other

forces includes separation, cohesion, alignment, roost attraction and flight forces. The topological interaction in this model makes each individual adapt its interaction range in such a way that individuals attempt to interact only with a constant number of their closest neighbors [9].

(10)

Cucker-Smale

Cucker-Smale postulated a model, both continuous and discrete, for the evolution of a flock. As a motivating example, they considered a population of birds or fish, whose members are moving in R3_{. They considered an adjacency matrix that contained the}

interaction information within the population. The velocity of a member is described by

vi(t + 1) − vi(t) = k

X

j=1

aij(vj(t) − vi(t)), (1)

where the bird adjusts its velocity viby the addition of a weighted average of the difference

between the velocity of the bird and the velocities of the other flock members. aij is a

weight that quantifies how the birds influence each other and is defined by

aij = η ||xi− xj||2 ,

where η is a continuous function converging to zero for increasing distances between bird members. These weights constitute the elements of the previously mentioned adjacency matrix. As Hildenbrandt et al., Cucker and Smale also made the assumption that the influence is a function of the distance between the birds. For the change of position, a natural equation is used in the shape of

x(t + 1) = x(t) + ∆tv(t), (2)

where x(t) and v(t) are the birds position and velocity at a time t, and ∆t is the time step [8].

Hayakawa

Hayakawa’s model is, in comparison with the other two models, a simplified model that is based on field measurements. The model introduces a one dimensional coordinate system for flying geese. Hayakawa assumes that each individual of the flock has a linear trajectory with a constant lateral spacing a, comparable to the length of a wingspan. The equation of motion Hawakawa uses includes the mass m of an individual, a self-driving force α0, an interaction force Fij between individual i and j, and a parameter γ including

a drag coefficient, the density of air, and the frontal area S of an individual.

m¨v = α0+

X

j<i

Fij − γ( ˙yi)2 (3)

Furthermore, Hayakawa also assumes that an individual only interacts with its nearest neighbors [7].

(11)

The Benefits of the TVM

All the models mentioned above have their own special ways of addressing the flocking problem. Hayakawa’s model treats a one dimensional flock and is, because of this, less complex. The model from Hildenbrandt et al. is three dimensional and have a lot of calculations for different forces. Since the model is so complex, an actual implementation would be demanding. The Cucker-Smale model is very much implementable - preferably in matlab due to its heavy use of linear algebra regarding the model matrices. The reason for using the TVM in this study, is the implementation of the model using object-oriented programming paradigm, where the use of matrices and differential equations is not impossible, but more demanding. Using object-oriented programming, it is easier to create an effective, low-leveled simulation.

2.2.2 The Rules of the TVM

When inspecting a flock of birds, the complex behavior and movements can be astonish-ing. A way to simplify this behavior is to create a list of rules that each individual has to follow. It is often easier to explain these rules and why they are good assumptions rather than explaining the real behavior of the flock. These rules will also help when building or implementing a model. When the rules are set up, you will have to look at how the single units in the flock interact with the rest of the individuals. The TVM that is used in this study is built on a collection of simple rules explaining how individuals interact with each other, as well as some necessary constraints on the flock system [1, 2]

1. The individuals of a flock are identical

2. A conservation law applies for the flock, meaning units cannot be created nor vanish during flight.

3. Individuals do not have an internal compass, i.e. they do not know which direction they are supposed to fly in.

4. Each individual of a large flock moves in a dimension d (R2) and it tries at all times to follow its neighbors.

5. The interactions between the units in the flock are short-ranged, meaning that the individuals only interact with neighbors within a given interaction-distance R, independent of the flock size L.

6. An individual i will interact with another individual j within an interaction range 0 < rij < R. If the individual j is within the range l0 < rij < R of individual i they

will attract each other. If individual j is within the range 0 < rij < l0 of individual

(12)

7. There is a short-ranged stochastic noise in the following between the individuals that will cause the interactions to not be ideal.

8. Each individual moves with a constant absolute velocity v0 and it can only change

its direction.

9. The model must have complete rotational symmetry, with the same probability for all directions.

To assume that each individual only interacts with other individuals within a given radius is a good approximation because each individual only has a limited view. This will indirect make the individual interact with the whole flock through chain reactions. Using this assumption, the simulation becomes less complex. The complementing l0 term

symbolizes the preferred internal distances within the flock between each flock member. When interaction occurs, l0 is the boundary distance for when the interaction changes

between attraction and repulsion, which will be further described in the mathematical formulation in the following subsection. The noise exists because the movements of the individuals are not perfect. There is always a chance an individual will leave the flock and it is not rational to believe that an individual knows exactly how the other individuals in the flock will act.

2.2.3 The Mathematical Formulation

The rules mentioned above, can be formulated mathematically as a time-discrete model. Each computational iteration using the model will be performed for a time step ∆t. Beginning with spatial description, the position of a bird i is updated using

ri(t + ∆t) = ri(t) + v0∆t (cos(θi(t)), sin(θi(t))) , (4)

where the positional update of the bird is its current position ri(t) plus the distance

traveled due to its velocity over time ∆t in its current direction. The directional update, θi(t + ∆t), for the bird is calculated using

θi(t + ∆t) = Θ 1 M M X j=1 (vj(t) + gij(t)) + ηi(t) ! (5)

The new direction is based on Θ, M, vj, gij and ηi, all depending on the time t. Here Θ

is the polar angle, M is the number of birds within the interaction distance R, vj is the

velocity of bird j, gij is the attraction between bird i and bird j, and ηi is the stochastic

noise.

(13)

vj = v0(cos(θj(t)), sin(θj(t))) (6)

M is found by comparing the absolute distance rij from bird i to each other bird j in

the interaction distance R. If a bird j is positioned within R of bird i, the information of its velocity is used to update the new direction of bird i.

rij = |ri− rj| < R (7)

The interaction vector gij comes from the equation

gij = g0 lo rij 3 − l0 rij 2! (ri− rj), (8)

where the vector difference (ri − rj) is directed from bird j to bird i, and g0 is the

interaction magnitude parameter. If the birds gets within the interaction distance R defined above, two kinds of interactions will occur. If the distance is larger than the preferred inter-flock distance l0, the factor

lo rij 3 −l0 rij 2

becomes negative due to the dominating squared term, causing the bird i to experience an attractive "force" to bird j, in the opposite direction of the vector difference (ri − rj). If the distance rij

is smaller than l0, the factor will be positive due to the dominance of the cubed term

resulting in an interaction vector having the same direction as the vector difference. The bird will experience this as a repulsive force directed from bird j. The coefficient function that is scaling the vector difference (ri − rj) can be considered a potential

function gp(rij) = g0 l0 rij 3 −l0 rij 2

. A potential of this nature, as shown in Figure 2.1, resembles a Lennard-Jones potential, which describes the interaction between neutral atoms and molecules in gases [10].

(14)

0 1 2 3 4 5 −1 0 1 2 3 4 5 6 7 8 9 10

Distance Between Two Birds in Pair

Potential as Function of Distance

Figure 2.1: Illustration of the potential g_p(rij) for l0 = 0.5 and g0 = 5 as a function of the

distance rij between an arbitrary pair of birds i and j. In this case, the interaction distance R

is left unspecified. For a non-zero value of R, the potential function would yield zero potential for R < r_ij. Note the change in sign of the potential when passing the preferred distance l₀, indicating the attraction-repulsion-boundary.

The noise term η_i is calculated using a parameter ∆v combined with an error term ei,

were ei is a random number uniformly distributed between [−1, 1].

η_i(t) = ∆v(cos(πei), sin(πei)) (9)

The function Θ(u) finds the principal argument of a vector u corresponding to some complex number z = ux+ iuy = eiΘ. Simply using arctan(u_uy_x) would not be ideal, since

arctan is not defined for ux = 0, which occurs for birds moving perpendicularly to the

x-axis.

2.2.4 Brief Discussion of Implementation

A few constraints on the implementation of the mathematical model has been applied in this study in order to limit time-complexity for simulation purposes. The domain in which the flocking birds are moving is square and has been given dimensions of 64 × 64 length units, or pixels. The domain has been given periodical boundary conditions, which

(15)

emulates a flock of infinite size. During the initial phase of a simulation, each generated bird is given a uniformly, pseudo-randomized positional vector within the domain as well as a starting direction θInitial in the boundary of θInitial ∈ [0, 2π]. The number of

simulated birds are consistently 512 for each simulation. The specific domain dimensions are further explained under the analysis chapter, but they are not essential to the study itself. When simulating, each discrete time step ∆t will be equal to 1.

2.3 Analysis

2.3.1 Limitations in the Tu-Vicsek Model

The TVM contains three key components - the influential, cohesive force of internal interaction, the velocities of the birds as well as the stochastic noise as a result of the decisions the birds make when correcting their courses. The model does not take the individual properties of the birds into account. All birds might not have an azimuthal field of vision corresponding to its entire surrounding, thus making the birds unaware of other bird neighbors traveling behind each other. Another aspect to consider in the model would be the one of weight and inertia. If this was to be implemented, the overall course corrections would be slower due to the inability for the birds to instantaneously change their directions of movement. In a model with incorporated inertia, each bird would have to perform work in order to stray from its current course, which is closer to a real-life scenario. Each bird would have a limited power and would have only be able to perform a limited amount of work. These could be ideas of improvements of the TVM, though this is not taken into account in this study.

2.3.2 Flock Density

The average flock density might be an indicator of sorts, to whether a starting condition of a simulation will be interesting or not. The domain average flock density, ρDA, is

simply calculated using

ρDA =

N

A, (10)

Here N is the number of birds in the system and A is the area of the simulation domain. Note that this is only favorable in some cases. Tu et al. found that a certain relation be-tween the average density and the preferred inter-flock distances could determine whether cluster formations would occur or not.

To further classify flocking behavior, a more stringent method has to be used. Since certain phases yield clustering bird formations, it is important to be able to locate these

(16)

discrete clusters in order to compute accurate flock densities. A usable method of finding these densities is using a uniform mesh of squares. In this study, the flocking domain has dimensions of 64 × 64, which has side-lengths that are proportional to the power of two. During each time-step each bird is located via the two-dimensional mesh. The sum of areas of each occupied mesh square represents the macroscopic flock size. Dividing the N number of birds with this area results in the current flock density. The total flock area, AF lock(t), is composed of the sum of each occupied mesh square following the paradigm

of

AF lock(t) =

X

i,j

Aij ∀Aij containing birds. (11)

As shown in Figure 2.2, the density calculation mesh for a 64 × 64 domain is composed of square cells. For the purpose of illustration, each mesh cell has dimensions of 4 × 4 which also is allowed according to our grid specifications. These cells are referred to via their positional indices (i, j), where i = j = 1 is the first cell starting by the origin. During an actual simulation, the mesh cells will be of dimensions 2 × 2.

0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 x y

Figure 2.2: Density calculation mesh for mesh cells with side length 4.

Figure 2.3 shows a three-dimensional plot associated with Figure 2.2 demonstrating the bird distribution over the mesh, where each vertex of the surface-squares corresponds to given (i, j) indices.

(17)

1 23 4 5 67 8 9 101112 13141516 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 160 5 10 15 20 25 i j

Birds Per Mesh Square

Figure 2.3: Bird distribution associated with the calculation mesh in Figure 2.2. Each vertex corresponds to the indicated indices on the axis.

The flock area AF lock varies between simulation steps. The momentaneous flock density,

ρt, is therefore a time dependent function following

ρt =

N AF lock(t)

, (12)

which is a better approximation of the current flock density, since it can manage trans-forming flock formations.

2.3.3 Measurement of Order

One idea behind quantifying the order of a flock could be to use some standard measure-ment of overall flock orientation. What is considered an ordered flock in this study, would be a flock where the majority of the flock members are moving in the same direction. This would be referred to as a flock with high polarization [11]. The polarization Φ is represented by the absolute normalized average velocity of the flock

Φ = 1 N X v_i ₍₁₃₎

(18)

resulting in a scalar Φ ∈ [0, 1]. Since the absolute velocity for each flock member is consistently v0 in the TVM, the equation is simplified as

Φ = 1 N v0 N X i=1 vi . (14)

Using the polarization it is possible to illustrate the varying degree of order within the flock over the course of a simulation.

2.3.4 Statistical Analysis

To see how a flock behaves, it is necessary to look at the statistical behavior of the flock. A couple of interesting things to look closer at is the mean velocity and its variance, mean direction and the average inter-flock distances. Since the model uses a constant absolute velocity, the velocity only differs by change in direction. The mean velocity tells which direction the flock is moving in general, and the variance shows how well the birds follow each other over a long period of time. If the variance is small, the possibility of finding a bird moving in a direction other than the whole flock due to the incorporated statistical noise is small. The equation for the mean velocity is

v = 1 N N X i=1 vi, (15)

where N is the total number of birds in the flock, and vi is the velocity of the i :th bird.

Further, finding the variance of the mean velocity is done using

V ar(¯v) = 1 N N X i=1 |vi − ¯v|2, (16)

where vi is the velocity for the i :th bird, ¯v the mean velocity for the flock and N is the

total number of flock members.

For an homogeneous flock where the flock members are evenly spread over the simulation domain, one way to find the mean inter-flock distances could be using

RM ean = 1 N2 N X i=1 N X j=1 |ri− rj| ∀i 6= j, (17)

where ri is the positional vector of bird i, rj is the positional vector of bird j, and N is

the total number of birds. Since homogeneity only occurs for a limited amount of flock states, it might be better to approach the inter-flock distance calculation with a method that circumvents the problem of cluster formations. Finding the mean distance of a

(19)

system with cluster formations is generally difficult since the clusters has to be separated and treated respectively. Instead, the mean distance will be approximated in a similar fashion as the momentaneous density defined by equation XXXX. In this case, the surface is divided in squares with size 2 × 2, which fill the whole surface 64 × 64. If a square contains birds, these birds are counted. For each square containing birds, the square root of the area divided by the number of birds in the square is calculated. The square roots are then summed and divided by the number of squares containing birds in order to get an average. 1 B B X m=1 r ASquare Nm . (18)

Here B is the number of squares that contain birds, ASquare is the area of a square, and

Nm is the number of birds in square m.

2.4 Results

Featured results from this study are solely based on various sets of parameter values in the TVM. The chapter is divided into subsections regarding specific results, described by their section title. Each simulation result is presented with its associated parameter set table used to generate the corresponding simulation data.

2.4.1 Ordered and Disordered Phases

Parameters

As previously mentioned, Vicsek noticed apparent ferromagnetic properties during his modeling from 1995 [1]. An ordered ferromagnetic flock can be obtained with use of the further developed TVM and the following model parameters. Note that for this default simulation, the noise parameter ∆v is equal to 0. Using this as a reference, ∆v will be incrementally increased until the state undergoes a phase transition. For all consecutive simulations, the interaction radius will selected as R = 3.

Table 2.1: Parameter values for achieved ferromagnetic ordered phase

Parameter g0 l0 v0 ∆v R N L

(20)

Evolution of an Ordered Flock

Figure 2.4 illustrates the initial and final states of a simulation, with parameter values from Table 2.1. 0 10 20 30 40 50 60 0 10 20 30 40 50 60

(a) Time step t = 1

0 10 20 30 40 50 60 0 10 20 30 40 50 60 (b) Time step t = 1000

Figure 2.4: Initial and final state of a flock of N = 512 birds, with a domain side-length of L = 64, g0 = 1, R = 3, v0 = 1, l0 = 1 and ∆v = 0. Note the heavy increase in flock alignment

after finished simulation.

The flock evolution is measured using the polarization (Figure 2.5 (a)), density (Figure 2.5 (b)), variance of mean velocity (Figure 2.6 (a)), and mean distance (Figure 2.6 (b)).

(21)

0 200 400 600 800 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time Polarization

(a) Polarization during simulation

0 200 400 600 800 1000 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 Time Density

(b) Density during simulation

Figure 2.5: Polarization and density illustrated during 1000 time steps for the simulation corresponding to the parameter set in Table 2.1. Note the high polarization due to the flock member alignment. The density also indicates that the flock has reached a more compact state since it increases towards the end of the simulation.

(22)

0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

Mean Velocity Variance

(a) Variance of mean velocity

0 200 400 600 800 1000 0.8 1 1.2 1.4 1.6 1.8 2 Timesteps

Mean Distance in Flock

(b) Mean distance

Figure 2.6: (a) Variance of mean velocity for the ferromagnetic flock with ∆v = 0.

(b) Evolution of mean distance during simulation with ∆v = 0 and parameters from Table 2.1. Note that the resulting mean distance converges to about 1, which is comparable to the preferred inter-flock distance l₀ = 1.

Effect of Increase in Noise

Following the previously mentioned temperature analogy, it is possible to see the im-pact of an increase in the stochastic noise parameter ∆v in the form of disorder during the course of simulation. Initially, ∆v is increased to 0.1 with earlier parameter values unchanged, see Figure 2.7.

(23)

0 10 20 30 40 50 60 0 10 20 30 40 50 60

(a) Time step t = 1

0 10 20 30 40 50 60 0 10 20 30 40 50 60 (b) Time step t = 1000

Figure 2.7: An increase of ∆v from 0 to 0.1 results in a slightly less ordered ferromagnetic flock. The bird alignment is still evident, but the internal cohesion seems to be lower as the flock is divided into separate groups.

Similar to the previous simulation, the flock evolution is measured using the polarization (Figure 2.8 (a)), density (Figure 2.8 (b)), variance of mean velocity (Figure 2.9 (a)), and mean distance (Figure 2.9 (b)).

(24)

0 200 400 600 800 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time Polarization

(a) Polarization during simulation

0 200 400 600 800 1000 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Time Density

(b) Density during simulation

Figure 2.8: Similar to the previous illustration, the polarization and density increases during simulation. An important difference is that the maximum values for both has decreased due to the increase in ∆v. Previously used parameters have been left unchanged.

0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(a) Variance of mean velocity

0 200 400 600 800 1000 0.8 1 1.2 1.4 1.6 1.8 2 Timesteps

(b) Mean distance

Figure 2.9: (a) The mean velocity variance ferromagnetic flock with increased ∆v = 0.1. (b) The mean distance for a simulation with ∆v = 0.1 is slightly increased compared to the previous illustration with ∆v = 0.

(25)

Phase Transition

Further increasing ∆v with increments of 0.1, Figure 2.10 illustrates the effect on the polarization where maximum polarization magnitude for a given simulation has been plotted to its associated value of ∆v starting from 0.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Noise Maximum Polarization

Figure 2.10: Maximum polarization for simulation with corresponding ∆v ∈ [0, 1]

As shown in Figure 2.10, the maximum polarization magnitude in a given simulation decreases as ∆v increases. Between ∆v = 0.6 and ∆v = 0.8, the polarization starts to reach a minimum. For ∆v >> 0.6, the polarization converges to a constant value of about 0.18. Figure 2.11 illustrates initial and final states of a simulation at this limit ∆v = 0.6.

(26)

0 10 20 30 40 50 60 0 10 20 30 40 50 60

(a) Time step t = 1

0 10 20 30 40 50 60 0 10 20 30 40 50 60 (b) Time step t = 1000

Figure 2.11: Initial and final state of a simulation with the approximate limit of ∆v = 0.6 found using results displayed in Figure 2.9. The final frame indicates that cohesion occurs, but that the flock members has started to lose its ability to align with each other. In other words, the flock is starting to display a paramagnetic behavior.

A comparison is made between the mean velocity variance of the ordered ferromagnetic and the disordered paramagnetic phases. The comparison, seen in Figure 2.12, shows the increase of mean velocity fluctuations due to the higher ∆v value for the disordered paramagnetic phase.

(27)

0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(a) ∆v = 0 0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(b) ∆v = 0.6

Figure 2.12: Variance of mean velocity for (a) a ferromagnetic phase with ∆v = 0 and (b) a paramagnetic phase for ∆v = 0.6 .

For the ordered ferromagnetic phase, the mean velocity variance decreases to an average converging to zero, which means that the risk of a mean velocity deviation is minimal for the flock as it reaches a stationary flight after a long period of time. In the simulation, a period of 200 time steps is enough for the flock to reach this state. As the noise parameter is increased, the mean velocity variance becomes more erratic and the mean velocity is more likely to diverge. Figure 2.13 depicts the mean velocity variation for the transition over ∆v = 0.2, 0.4, 0.6 and 0.8.

(28)

0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(a) ∆v = 0.2 0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(b) ∆v = 0.4 0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(c) ∆v = 0.6 0 200 400 600 800 1000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time

(d) ∆v = 0.8

Figure 2.13: Increase in ∆v indicates a higher mean velocity variance. (a) For ∆v = 0.2, the variance converges to a value just above zero, meaning the system is more sensitive to mean velocity fluctuations than for ∆v = 0 in Figure 2.11.(a). (b) ∆v = 0.2 yields a higher variance with the exception of the pronounced peak occurring approximately between 310 ≤ t ≤ 550. (c) and (d) displays similar properties where an increase from ∆v = 0.6 to ∆v = 0.8 has no substantial impact.

Finding the mean distances corresponding to Figure 2.13, the noise impact on the flock cohesion can be visualized, as seen in Figure 2.14.

(29)

0 200 400 600 800 1000 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 Timesteps

(a) ∆v = 0.2 0 200 400 600 800 1000 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 Timesteps

(b) ∆v = 0.4 0 200 400 600 800 1000 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 Timesteps

(c) ∆v = 0.6 0 200 400 600 800 1000 1.72 1.74 1.76 1.78 1.8 1.82 1.84 1.86 Timesteps

(d) ∆v = 0.8

Figure 2.14: The figures show the variation of the mean distance within the flocks for each time step for four different values of ∆v. In (a) the mean distance goes to the preferred value l0 = 1. For all values of ∆v, the mean distances within the flock decline, but when ∆v is

increased, as in (b),(c) and (d) it is possible to see that mean distance fails to tend to l0 = 1.

The graphs also show that the mean distance varies more for a higher ∆v.

2.4.2 Other Flock Phases

A number of other interesting flock behaviors can be obtained with some modifications to the model parameters. In Figure 2.15, simulated flocks corresponding to parameters

(30)

in Table 2.2 are shown. Note that each still frame represents no specific time step during simulation. 0 10 20 30 40 50 60 0 10 20 30 40 50 60 (a) 0 10 20 30 40 50 60 0 10 20 30 40 50 60 (b) 0 10 20 30 40 50 60 0 10 20 30 40 50 60 (c) 0 10 20 30 40 50 60 0 10 20 30 40 50 60 (d)

Figure 2.15: Noteworthy states achieved with the associated parameters in Table 2.2. Table 2.2: Parameter values for achieved ferromagnetic ordered phase

Parameter g0 l0 v0 ∆v R N L

(a) 20 1 1 0.2 3 512 64

(b) 1 5 1 0.2 20 512 64

(c) 3 0.5 1 0.2 3 512 64

(31)

Two interesting properties apparent for simulation (a) and (b), in Figure 2.15, is that the generated clusters remain static during simulation, with only a slight vibration due to the non-zero value of ∆v. Case (a) resembles a solid crystal with amorphous lattice points.

2.5 Discussion

2.5.1 Interpreting the Statistical Noise

The TVM incorporates statistical noise in the form of η_i. Toner et al. suggested that the noise comes from the inability of the birds to make consistent decisions when choosing their movement direction [2], which has an impact on the evolution of a simulation. In the TVM, the magnitude of the stochastic noise is controlled with ∆v. A higher value of ∆v results in a bird behavior similar to that of vibrating molecules or atoms in matter undergoing temperature increase due to an influx of heat. Using this similarity, it is possible to increase disorder within the flock via perturbation due to an increase in ∆v -almost like a system of particles experiencing heating. A ferromagnet in such a situation will, under incremental heating, eventually lose its ferromagnetic properties due to the decrease in magnetic spin alignment, resulting in a paramagnetic phase. The magnet will transition into the paramagnetic phase at a certain critical temperature TC called the

Curie Temperature [12]. In this study, the noise parameter ∆v is the foremost variable to consider when phase transition from an ordered state into a disordered is desired.

Investigating the data from the results, it possible to see that changing the noise param-eter has a large effect in the interval [0, 0.6], as shown in Figure 2.10. The effect seems to have an approximate 1 − C∆v2 _{dependence up to about ∆v = 0.6, meaning the}

po-larization decreases when the noise increases. For ∆v >> 0.6, the popo-larization seems to be approximately invariant. This effect can be compared to the high-temperature limit in statistical physics, where, for instance, Cv, the heat capacity at constant volume, goes

to 3N Kb [12].

In figure 2.14 it is possible to see that the mean distance within the flock declines. This is a sign of the flocking behavior of the birds. When ∆v = 0.2 (as in (a)) the mean distance even converges to the preferred inter flock distance l0 = 1. Because of the increasing ∆v

in (b),(c),(d) the birds are subject to a larger noise. This is shown by the increase of the mean distance within the flock. The effect is that the birds have a difficulty forming flocks that stay together, i.e. the birds tend to try to fly away from the flock. Inspecting the mean velocity variance in Figure 2.14, it is indicated that the system becomes more sensitive to mean velocity deviation as ∆v increases for longer time periods of flights. For ∆v = 0, the variance converges to about 0, which indicates exceptionally low risk for sudden mean velocity deviation. For the solid crystal phase as shown in Figure 2.15 (a), it might be possible to qualitatively analyze the impact of ∆v by examining phonon energies for the crystal lattice. This, however, goes beyond the scope of this study.

(32)

With the parameters g0 = 1, v0 = 1, l0 = 1, R = 3, ∆v = 0 it is possible to simulate a

growing flock that is rapidly aligning. The birds go from uniformly randomized positions with random starting-directions Θ ∈ [0, 2π] to a cluster of a couple of hundred birds with almost identical movement directions [13]. These cluster formations are similar to formations that can be seen in the nature (Figure 2.16).

Figure 2.16: Real-life occurrences of flocking behavior. Left image depicts a flock of starlings in flight near Athens, Greece. Right image shows snow geese in Merced National Wildlife Refuge, California [14, 15].

By setting the interaction parameter g0 to 20 the behavior changes completely. Now the

simulation renders multiple static small flocks with constant flock distances (Figure 2.15 (a)).

2.5.2 Transition from Ferromagnetic to Paramagnetic Phase

With the right set of parameters, it is possible to get a ferromagnetic looking state of the flock. Shown in Figures 2.4-9 are the result of a parameter sets, where l0 is about

66 percent of R which makes the attraction area quite small. Because of the large area where the birds repel each other, the birds keep the distance between each other but the interaction causes the birds to travel in almost the same direction. The ferromagnetic state is the most obvious for low values of ∆v, where the magnitude of the polarization is large. If ∆v is increased, the magnitude of the polarization decreases and the order of the flock becomes more and more chaotic. This effect is apparent up to ∆v = 0.6. When ∆v is increased further, the polarization is approximately constant, see Figure 2.10. The constant polarization for high values of the noise implies that only ∆v of the interval [0, 0.6] is of interest. When the flock order decreases, we see that the ferromagnetic state devolves to a paramagnetic state. This effect has an evident relationship to the real ferromagnetic-paramagnetic transition at the critical Curie-temperature Tc [12].

(33)

Chapter 3 Summary and Conclusion

With the established rules, it is possible to set up a flocking model. The model not only describes a flock of birds, but also describes the different phases of the flock. Simulations have shown that there are several flock phases that can be compared with the phases of matter. By using methods from the statistical physics the flock phases have been identified. If these phases appear in real flocks of birds is left unsaid. Even though a small scale system was used, results from earlier works have successfully been reproduced. The results from the report confirms that the model works for a small scale system with only 512 birds in a 64 × 64 system. This can be beneficial for people who do not have the possibility to run large scale simulations, but still want to explore the phenomenon of flocking behavior. The analytical tools used in this study have successfully shown that variations in the model parameters provide an impact on the flocking behavior within the system, meaning they should be viable in proving the occurrence of flocking behavior in models other than the TVM used in this report as well.

Something that must be noted is that the model does not only describe behavior of flock-ing birds, but also describes the behavior of many different kinds of flockflock-ing, for example fish schools and microorganisms. If the conservation law is ignored and a stochastic increase of the individuals is implemented, even the flocking behavior of bacteria could be simulated. Further interesting work with the model would be to implement obstacle avoidance or a predator-prey behavior. This would give the model another dimension that could be of interest to describe for example a herd of buffalo. The set of analytical tools could be complemented with the pair correlation function or the radial distribution function in order to verify the correlation between intermediate flock members due to interaction.

(34)

Bibliography

[1] A. Czirók and T. Vicsek, “Collective motion,” in Statistical Mechanics of Biocom-plexity (D. Reguera, J. Vilar, and J. Rubí, eds.), vol. 527 of Lecture Notes in Physics, pp. 152–164, Springer Berlin Heidelberg, 1999.

[2] J. Toner, Y. Tu, and S. Ramaswamy, “Hydrodynamics and phases of flocks,” Annals of Physics, vol. 318, no. 1, pp. 170 – 244, 2005. Special Issue.

[3] J. T. Emlen, “Flocking behavior in birds,” The Auk, vol. 69, no. 2, pp. 160–170, 1952.

[4] C. W. Reynolds, “Flocks, herds and schools: A distributed behavioral model,” SIG-GRAPH Comput. Graph., vol. 21, pp. 25–34, Aug. 1987.

[5] A. Huth and C. Wissel, “The simulation of the movement of fish schools,” Journal of Theoretical Biology, vol. 156, no. 3, pp. 365 – 385, 1992.

[6] Y. Tu, J. Toner, and M. Ulm, “Sound waves and the absence of galilean invariance in flocks,” Phys. Rev. Lett., vol. 80, pp. 4819–4822, May 1998.

[7] Y. Hayakawa, “Spatiotemporal dynamics of skeins of wild geese,” EPL (Europhysics Letters), vol. 89, no. 4, p. 48004, 2010.

[8] F. Cucker and S. Smale, “Emergent behavior in flocks,” Automatic Control, IEEE Transactions on, vol. 52, pp. 852–862, May 2007.

[9] H. Hildenbrandt, C. Carere, and C. Hemelrijk, “Self-organized aerial displays of thousands of starlings: a model,” Behavioral Ecology, vol. 21, no. 6, pp. 1349–1359, 2010.

[10] J. E. Jones, “On the determination of molecular fields. ii. from the equation of state of a gas,” Proceedings of the Royal Society of London. Series A, vol. 106, no. 738, pp. 463–477, 1924.

[11] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale, “Scale-free correlations in starling flocks,” Proceedings of the National Academy of Science, vol. 107, pp. 11865–11870, June 2010.

(35)

[12] C. Kittel, Introduction to Solid State Physics, Eighth Edition. John Wiley & Sons, Inc., 2005.

[13] L. Almgren and J. Bergdahl, “Simulation of modified vicsek model.” https:// vimeo.com/95856614.

[14] Flickr, “Starlings near athens.” https://www.flickr.com/photos/mwf2005/ 3092761353/.

[15] Flickr, “Snow geese in merced national wildlife refuge.” https://www.flickr.com/ photos/stevecorey/12723335535/.

Quantitative Analysis of Physical and Statistical Properties of Flocks.