Source code for flocking

import sys
import copy
from math import *
import numpy as np
import numpy.random as random
import matplotlib.pyplot as plt
import matplotlib.animation as ani

#
# input file keywords
#
# note: global variables like these are
# often written in ALL CAPS
#
COHESION_TAG = "cohesion"
SEPARATION_TAG = "separation"
ALIGNMENT_TAG = "alignment"
AVOIDANCE_TAG = "avoidance"
CHASE_TAG = "chase"
MAX_TAG = "max"
RADIUS_TAG = "radius"
STRENGTH_TAG = "strength"

PREY_TAG = "prey"
PREDATOR_TAG = "predator"
SKY_TAG = "sky"
TIME_TAG = "time"
N_TAG = "count"
SPEED_TAG = "speed"
D_TAG = "dimension"
SIZE_TAG = "size"
TOTAL_TAG = "total"
DT_TAG = "dt"
RDT_TAG = "recording-dt"

STRENGTH_INDEX = 0
MAX_INDEX = 1
RADIUS_INDEX = 2
N_INDEX = 0
SPEED_INDEX = 1
D_INDEX = 0
SIZE_INDEX = 1
DT_INDEX = 0
TOTAL_INDEX = 1
RDT_INDEX = 2


[docs]def find_info(lines, tag): """ Searches for the information wrapped in the given tag among the given lines of text. If tag is, e.g., "foo", the function searches for the start tag <foo> and the end tag </foo> and returns the lines of information between them. The function only finds the first instance of the given tag. However, in order to catch multiple instances of the tag, the function also returns all the information in lines following the end tag. For instance, if lines contains the strings: .. code-block :: aa <foo> bb cc </foo> dd ee ff the function will return two lists: ["bb", "cc"], ["", "dd", "ee", "ff"]. Args: lines (list): the information as a list of strings tag (str): the tag to search Returns: list, list: the lines between start and end tags, the lines following the end tag """ info = [] is_relevant = False line_number = 0 # go through the data for i in range(len(lines)): line = lines[i] if is_relevant: # if we have found the starting tag, record information info.append(line) contents = line.strip() # remove whitespace at the start and end of the line if len(contents) > 0: # skip empty lines if contents[0] == "<" and contents[-1] == ">": # is this a tag? if contents[1:-1] == tag: # found the starting tag if not is_relevant: # we had not yet found the tag is_relevant = True # the following lines are relevant line_number = i else: # we had already started this tag print("Found tag <"+tag+"> while already reading <"+tag+">") raise Exception("parsing error") if contents[1:-1] == "/"+tag: # found the end tag return info, lines[i+1:] # we end up here, if we reach the end of the file if is_relevant: # the file ends while reading info (start tag was found, but no end tag) print("Reached the end of file while parsing <"+tag+"> from line "+str(line_number+1)) raise Exception("parsing error") elif info == []: # the tag was not found print("Tag <"+tag+"> was not found") return [], lines
[docs]def parse_line(line): """ Separates tag and info on a line of text. The function also removes extra whitespace and comments separated with #. For instance if line is " x : 1.23 # the x coordinate", the function returns ("x", "1.23"). Args: line (str): a string of information Returns: str, str: tag, info """ parts = line.split(":") tag = "" info = "" if len(parts) > 1: tag = parts[0].strip() info = parts[1].split("#")[0].strip() return tag, info
[docs]def read_force_info(lines, default=1.0): """ Reads parameters from given lines. Args: lines (list): information as a list of strings default (float): the default lattice parameter in all directions Returns: list: parameters """ parameters = [default]*3 for line in lines: tag, info = parse_line(line) if tag == STRENGTH_TAG: parameters[STRENGTH_INDEX] = float(info) elif tag == RADIUS_TAG: parameters[RADIUS_INDEX] = float(info) elif tag == MAX_TAG: parameters[MAX_INDEX] = float(info) return parameters
[docs]def read_animal_info(lines, default=1): """ Reads parameters from given lines. Args: lines (list): information as a list of strings default (float): the default lattice parameter in all directions Returns: list: parameters """ parameters = [default]*2 for line in lines: tag, info = parse_line(line) if tag == N_TAG: parameters[N_INDEX] = float(info) elif tag == SPEED_TAG: parameters[SPEED_INDEX] = float(info) return parameters
[docs]def read_box_info(lines, default=2): """ Reads parameters from given lines. Args: lines (list): information as a list of strings default (float): the default lattice parameter in all directions Returns: list: parameters """ parameters = [default]*2 for line in lines: tag, info = parse_line(line) if tag == D_TAG: parameters[D_INDEX] = float(info) elif tag == SIZE_TAG: parameters[SIZE_INDEX] = float(info) return parameters
[docs]def read_time_info(lines, default=0.1): """ Reads parameters from given lines. Args: lines (list): information as a list of strings default (float): the default lattice parameter in all directions Returns: list: parameters """ parameters = [default]*3 for line in lines: tag, info = parse_line(line) if tag == DT_TAG: parameters[DT_INDEX] = float(info) elif tag == TOTAL_TAG: parameters[TOTAL_INDEX] = float(info) elif tag == RDT_TAG: parameters[RDT_INDEX] = float(info) return parameters
[docs]def read_dynamics_file(filename): """ Reads the given file for force parameters. Args: filename (str): the file to read Returns: lists: cohesion, alignment, separation, avoidance, chase parameters """ f = open(filename) lines = f.readlines() f.close() info, dummy = find_info( lines, COHESION_TAG ) cohesion_params = read_force_info( info ) info, dummy = find_info( lines, SEPARATION_TAG ) separation_params = read_force_info( info ) info, dummy = find_info( lines, ALIGNMENT_TAG ) alignment_params = read_force_info( info ) info, dummy = find_info( lines, AVOIDANCE_TAG ) avoidance_params = read_force_info( info ) info, dummy = find_info( lines, CHASE_TAG ) chase_params = read_force_info( info ) return cohesion_params, alignment_params, separation_params, avoidance_params, chase_params
[docs]def read_system_file(filename): """ Reads the given file for system parameters. Args: filename (str): the file to read Returns: lists: prey, predator, box, time parameters """ f = open(filename) lines = f.readlines() f.close() info, dummy = find_info( lines, PREY_TAG ) prey_params = read_animal_info( info ) info, dummy = find_info( lines, PREDATOR_TAG ) predator_params = read_animal_info( info ) info, dummy = find_info( lines, SKY_TAG ) box_params = read_box_info( info ) info, dummy = find_info( lines, TIME_TAG ) time_params = read_time_info( info ) return prey_params, predator_params, box_params, time_params
[docs]def draw(frame, xtraj, ytraj, ztraj, bounds): """ Draws a representation of the system as a scatter plot. Used for animation. Args: frame (int): index of the frame to be drawn xtraj (array): x-coordinates of all animals at different animation frames ytraj (array): y-coordinates at all animals at different animation frames ztraj (array): z-coordinates at all animals at different animation frames bounds (array): list of lower and upper bounds for the plot as [[xmin, xmax], [ymin, ymax]] """ plt.clf() ax = plt.axes() ax.set_xlim(bounds[0]) ax.set_ylim(bounds[1]) ax.set_aspect('equal') size = (ztraj[frame]+1)*2 plt.scatter(xtraj[frame], ytraj[frame], marker='o', s=size )
[docs]def animate( prey, predators, box, multiply = [3,3] ): """ Animates the simulation. Args: prey (list): list of :class:`Prey` objects predators (list): list of :class:`Predator` objects box (flocking.PeriodicBox): supercell multiply (array): number of periodic images to draw in x and y directions """ nframes = len(prey[0].trajectory) print("animating "+str(nframes)+" frames") xtraj = [] ytraj = [] ztraj = [] Lx = box.lattice[0] Ly = box.lattice[1] # number of periodic images in x and y directions multix = multiply[0] multiy = multiply[1] margin = 0.1 bounds = np.zeros([3,2]) bounds[0,1] = Lx*multix bounds[1,1] = Ly*multiy bounds[:,0] -= margin bounds[:,1] += margin for i in range(nframes): xtraj.append([]) ytraj.append([]) ztraj.append([]) for p in prey: box.shift_inside_box(p.trajectory[i]) for ix in range(multix): for iy in range(multiy): xtraj[-1].append(p.trajectory[i][0] + ix*Lx) ytraj[-1].append(p.trajectory[i][1] + iy*Ly) ztraj[-1].append(1) for p in predators: box.shift_inside_box(p.trajectory[i]) for ix in range(multix): for iy in range(multiy): xtraj[-1].append(p.trajectory[i][0] + ix*Lx) ytraj[-1].append(p.trajectory[i][1] + iy*Ly) ztraj[-1].append(10) xtraj=np.array(xtraj) ytraj=np.array(ytraj) ztraj=np.array(ztraj) fig = plt.figure() motion = ani.FuncAnimation(fig, draw, nframes, interval=10, fargs=(xtraj, ytraj, ztraj, bounds) ) plt.show()
[docs]class PeriodicBox: """ Class representing a simulation box with periodic boundaries. The box is orthogonal, i.e., a rectangular volume. As such, it is specified by the lengths of its edges (lattice constants). Args: lattice (array): lattice constants """ def __init__(self, lattice): self.lattice = lattice
[docs] def shift_inside_box(self, position): """ If the given position (3-vector) is outside the box, it is shifted by multiple of lattice vectors until the new position is inside the box. That is, the function transforms the position vector to an equivalen position inside the box. Args: position (array): the position to be shifted """ # go over x, y and z coordinates for i in range(3): while position[i] < 0: position[i] += self.lattice[i] while position[i] > self.lattice[i]: position[i] -= self.lattice[i]
[docs] def distance_squared(self, position1, position2): """ Calculates and returns the square of the distance between two points, .. math :: r^2_{ij} = (x_i-x_j)^2 + (y_i-y_j)^2 + (z_i-z_j)^2. In a periodic system, each boid has an infinite number of periodic copies. Therefore the distance between two points is not unique. The function returns the shortest such distance, that is, the distance between the the periodic copies which are closest ot each other. Args: position1 (array): the first point position2 (array): the second point Returns: float: the squared distance :math:`r^2_{ij}` """ pass v = self.vector(position1,position2) return v @ v
[docs] def vector(self, position1, position2): """ Returns the vector pointing from position1 to position2, .. math :: \\vec{r}_{i \\to j} = \\vec{r}_j - \\vec{r}_i In a periodic system, each boid has an infinite number of periodic copies. Therefore the displacement between two points is not unique. The function returns the shortest such displacement vector. Args: position1 (array): the first point position2 (array): the second point Returns: array: components of :math:`\\vec{r}_{i \\to j}`, :math:`[x_{i \\to j}, y_{i \\to j}, z_{i \\to j}]` """ vector_1to2 = position2 - position1 # loop over x, y, and z directions for i in range(3): # If the absolute value of the separation # in this direction is larger than half the length # of the simulation box, there must be another # periodic image of boid2 closer to boid 1. # # We can find it by translating the coordinates by # an integer multiple of lattice vectors. # # Note: there is a more efficient way to calculate this. while vector_1to2[i] < -self.lattice[i]/2: vector_1to2[i] += self.lattice[i] while vector_1to2[i] > self.lattice[i]/2: vector_1to2[i] -= self.lattice[i] return vector_1to2
[docs]class Animal: """ A class representing an animal. All animals are treated as point-like agents in this model. Args: position (array): coordinates :math:`[x, y, z]` velocity (array): velocity components :math:`[v_x, v_y, v_z]` top_speed (fload): the maximum speed allowed for this animal mass (float): mass :math:`m` """ def __init__(self, position, velocity, top_speed, mass = 1.0): self.position = np.array(position) self.velocity = np.array(velocity) self.mass = mass self.force = np.zeros(3) self.trajectory = [] self.neighbors = [] self.top_speed = top_speed
[docs] def move(self, dt): """ Move the animal. Args: dt (float): time step :math:`\\Delta t` """ self.position += self.velocity * dt
[docs] def accelerate(self, dt): """ Set a new velocity as .. math:: \\vec{v}(t+\\Delta t) = \\vec{v}(t) + \\frac{1}{2m}\\vec{F} \Delta t If the speed would exceed the maximum speed of the animal, it is scaled down to the max value. Similarly, the animal must move at least with a speed that is half of the maximum speed. Args: dt (float): time step :math:`\\Delta t` """ self.velocity += self.force * dt/self.mass v2 = self.velocity @ self.velocity if v2 > self.top_speed**2: self.velocity *= self.top_speed/np.sqrt(v2) elif v2 < self.top_speed**2/4: self.velocity *= self.top_speed/np.sqrt(v2)/2
[docs] def save_position(self): """ Save the current position. Note: in a real large-scale simulation one would never save trajectories in memory. Instead, these would be written to a file for later analysis. """ self.trajectory.append( [ self.position[0], self.position[1], self.position[2] ] )
[docs]class Predator(Animal): """ A point-like model for a predator :class:`Animal`. Likes to harrass :class:`Prey` animals. Args: position (array): coordinates :math:`[x, y, z]` velocity (array): velocity components :math:`[v_x, v_y, v_z]` top_speed (fload): the maximum speed allowed for this animal mass (float): mass :math:`m` """ def __init__(self, position, velocity, prey, top_speed, mass = 1.0): Animal.__init__(self, position, velocity, top_speed, mass) # pick a random prey as the initial target self.target = random.choice( prey )
[docs] def pick_target(self, prey, box, radius): """ Randomly picks one :class:`Prey` within the given radius as the new target to chase. If there are no prey within range, nothing happens. Args: prey (list): list of :class:`Prey` objects box (flocking.PeriodicBox): supercell radius (float): maximum range for choosing a target """ targets = [] for boid in prey: dist = box.distance_squared( self.position, boid.position ) if dist < radius*radius: targets.append( boid ) if len(targets) > 0: self.target = random.choice( targets )
[docs] def add_chase_force(self, box, parameters): """ Adds a force driving the :class:`Predator` towards its target :class:`Prey`. Args: box (flocking.PeriodicBox): supercell parameters (list): force parameters """ chase_strength = parameters[STRENGTH_INDEX] max_force = parameters[MAX_INDEX] to_target = box.vector(self.position, self.target.position) force = chase_strength*to_target f2 = force @ force + 0.001 if f2 > max_force*max_force: force *= max_force/np.sqrt(f2) self.force += force
[docs] def add_separation_force(self, box, parameters): """ Adds a force to make predators avoid collisions. This is a short ranged repulsive force applied only if the animals are very close to each other. Args: box (flocking.PeriodicBox): supercell parameters (list): force parameters """ separation_strength = parameters[STRENGTH_INDEX] min_distance = parameters[RADIUS_INDEX] max_force = parameters[MAX_INDEX] force = np.zeros(3) for boid in self.neighbors: to_boid = box.vector(self.position, boid.position) d2 = to_boid @ to_boid + 0.001 if d2 < min_distance*min_distance: force -= separation_strength*to_boid/d2 f2 = force @ force + 0.001 if f2 > max_force*max_force: force *= max_force/np.sqrt(f2) self.force += force
[docs] def calculate_forces(self, box, force_parameters): """ Calculates all forces acting on this animal. The force sum is saved in Predator.force. Args: box (flocking.PeriodicBox): supercell force_parameters (list): force parameters """ separation_parameters = force_parameters[2] chase_parameters = force_parameters[4] self.force = np.zeros(3) self.add_chase_force(box, chase_parameters) self.add_separation_force(box, separation_parameters)
[docs]class Prey(Animal): """ A point-like model for a flock-forming :class:`Animal`. Can be harrassed by a :class:`Predator`. Args: position (array): coordinates :math:`[x, y, z]` velocity (array): velocity components :math:`[v_x, v_y, v_z]` top_speed (fload): the maximum speed allowed for this animal mass (float): mass :math:`m` """ def __init__(self, position, velocity, top_speed, mass = 1.0): Animal.__init__(self, position, velocity, top_speed, mass)
[docs] def add_cohesion_force(self, box, parameters): """ Adds a force to make :class:`Prey` seek each other. Cohesion force is a force pointing to the center (average coordinate) of all the neighbors of the animal. .. note :: This function is incomplete! Args: box (flocking.PeriodicBox): supercell parameters (list): force parameters """ cohesion_strength = parameters[STRENGTH_INDEX] max_force = parameters[MAX_INDEX] to_center = np.zeros(3) n_neighbors = len(self.neighbors) if n_neighbors == 0: return # calculate the vector from the animal to the center of its neighbors # normalize to_center and multiply by the strength parameter # todo force = 0.0 # check that the force is not too large f2 = force @ force + 0.001 if f2 > max_force*max_force: force *= max_force/np.sqrt(f2) self.force += force
[docs] def add_alignment_force(self, box, parameters): """ Adds a force to make :class:`Prey` move in the same direction Alignment force is calculated by first determining the average velocity vector of the neighbors. The force is then proportional to the difference between the current velocity of the animal and the calculated average. Args: box (flocking.PeriodicBox): supercell parameters (list): force parameters """ alignment_strength = parameters[STRENGTH_INDEX] max_distance = parameters[RADIUS_INDEX] max_force = parameters[MAX_INDEX] velocity = np.zeros(3) n_neighbors = len(self.neighbors) if n_neighbors == 0: return # calculate the average of the neighbor velocities for boid in self.neighbors: d2 = box.distance_squared(self.position, boid.position) if d2 < max_distance*max_distance: velocity += boid.velocity velocity /= n_neighbors # calculate difference between the velocity of the animal and the average # also multiply by the strength parameter to_velo = alignment_strength*(velocity-self.velocity) # normalize to_velo to_velo = to_velo / np.sqrt(to_velo @ to_velo + 0.001) # check that the force is not too large force = alignment_strength*to_velo f2 = force @ force + 0.001 if f2 > max_force*max_force: force *= max_force/np.sqrt(f2) self.force += force
[docs] def add_separation_force(self, box, parameters): """ Adds a force to make prey avoid collisions. This is a short ranged repulsive force applied only if the animals are very close to each other. Args: box (flocking.PeriodicBox): supercell parameters (list): force parameters """ separation_strength = parameters[STRENGTH_INDEX] min_distance = parameters[RADIUS_INDEX] max_force = parameters[MAX_INDEX] force = np.zeros(3) for boid in self.neighbors: to_boid = box.vector(self.position, boid.position) d2 = to_boid @ to_boid + 0.001 # apply repulsive force for too close neighbors # this implementation uses 1/r force if d2 < min_distance*min_distance: force -= separation_strength*to_boid/d2 # check that the force is not too large f2 = force @ force + 0.001 if f2 > max_force*max_force: force *= max_force/np.sqrt(f2) self.force += force
[docs] def add_avoidance_force(self, predators, box, parameters): """ Adds a force to make :class:`Prey` flee from :class:`Predator` animals. This force is calculated similarly to :meth:`Prey.add_separation_force`, but it is triggered by predators instead of other prey. Typically, this force should also be stronger and have a longer range than the separation force. Args: box (flocking.PeriodicBox): supercell parameters (list): force parameters """ avoid_strength = parameters[STRENGTH_INDEX] min_distance = parameters[RADIUS_INDEX] max_force = parameters[MAX_INDEX] force = np.zeros(3) for predator in predators: to_predator = box.vector(self.position, predator.position) d2 = to_predator @ to_predator # apply repulsive force for too close predators # this implementation uses 1/r force if d2 < min_distance*min_distance: force -= avoid_strength*to_predator/d2 f2 = force @ force + 0.001 # check that the force is not too large if f2 > max_force*max_force: force *= max_force/np.sqrt(f2) self.force += force
[docs] def calculate_forces(self, predators, box, force_parameters): """ Calculates all forces acting on this animal. The force sum is saved in Prey.force. Args: predators (list): list of :class:`Predator` objects box (flocking.PeriodicBox): supercell force_parameters (list): force parameters """ cohesion_parameters = force_parameters[0] alignment_parameters = force_parameters[1] separation_parameters = force_parameters[2] avoid_parameters = force_parameters[3] self.force = np.zeros(3) self.add_cohesion_force(box, cohesion_parameters) self.add_alignment_force(box, alignment_parameters) self.add_separation_force(box, separation_parameters) self.add_avoidance_force(predators, box, avoid_parameters)
[docs]def update_neighbors(animals, box, radius): """ Finds for each animal the other animals of the same type within the given radius. The animals can be both :class:`Prey` or :class:`Predator`. The neighbors are saved in the animal object as a list of animal objects. Args: animals (list): list of animals box (flocking.PeriodicBox): supercell radius (float): maximum distance to neighbors """ neighbors = [] r2 = radius*radius n_animals = len(animals) # empty the neighbor lists for i in range(n_animals): boid_i = animals[i] boid_i.neighbors = [] # loop over all pairs of animals and add new neighbors # # note: This is the most computationally expensive step in the whole # simulation. There are ways to make it faster, but this simple implementation # is sufficient for us. for i in range(n_animals): for j in range(i+1,n_animals): boid_i = animals[i] boid_j = animals[j] # if the animals are close to each other, add them as neighbors to both if box.distance_squared(boid_i.position, boid_j.position) <= r2: boid_i.neighbors.append(boid_j) boid_j.neighbors.append(boid_i)
[docs]def update_targets(prey, predators, box, radius): """ Possible updates the target for every predator. For each :class:`Predator`, there is a 10 % chance it seeks a new target :class:`Prey` via :meth:`Predator.pick_target`. Args: prey (list): list of :class:`Prey` objects predators (list): list of :class:`Predator` objects box (flocking.PeriodicBox): supercell radius (float): maximum distance to new targets """ for p in predators: r = random.random() if r < 0.1 or p.target == None: p.pick_target(prey, box, radius)
[docs]def calculate_forces(prey, predators, box, force_parameters): """ Calculates forces on all animals. Args: prey (list): list of :class:`Prey` objects predators (list): list of :class:`Predator` objects box (flocking.PeriodicBox): supercell force_parameters (list): force parameters """ radius = force_parameters[0][RADIUS_INDEX] update_neighbors(prey, box, radius) update_neighbors(predators, box, radius) update_targets(prey, predators, box, radius) for p in prey: p.calculate_forces(predators, box, force_parameters) for p in predators: p.calculate_forces(box, force_parameters)
[docs]def create_random_flock(prey_parameters, predator_parameters, box): """ Creates a flock of :class:`Prey` and possible also :class:`Predator` animals. The animals are placed randomly and they are given random velocities. Args: prey_parameters (list): amount and top speed of prey predator_parameters (list): amount and top speed of predators box (flocking.PeriodicBox): supercell """ n_prey = prey_parameters[N_INDEX] n_predator = predator_parameters[N_INDEX] s_prey = prey_parameters[SPEED_INDEX] s_predator = predator_parameters[SPEED_INDEX] # check if this is a 2D or 3D simulation if box.lattice[2] < 1: dimension = 2 else: dimension = 3 flock = [] predators = [] # create prey while len(flock) < n_prey: r = np.zeros(3) v = np.zeros(3) for i in range(dimension): r[i] = random.random()*box.lattice[i] v[i] = (2.5*random.random()-1)*s_prey/2 # don't put the prey too close to each other ok = True for boid in flock: if box.distance_squared(r, boid.position) < 0.5: ok = False break if ok: flock.append( Prey( r, v, s_prey ) ) # create predators while len(predators) < n_predator: r = np.zeros(3) v = np.zeros(3) for i in range(dimension): r[i] = random.random()*box.lattice[i] v[i] = (2.5*random.random()-1)*s_predator/2 ok = True for p in predators: if box.distance_squared(r, p.position) < 0.5: ok = False break if ok: predators.append( Predator( r, v, flock, s_predator ) ) return flock, predators
[docs]def update_positions_no_force(prey, predators, dt): """ Update the positions of all animals. .. math:: \\vec{r}(t+\\Delta t) = \\vec{r}(t) + \\vec{v} \Delta t Args: prey (list): a list of :class:`Prey` objects predators (list): a list of :class:`Predator` objects dt (float): time step :math:`\\Delta t` """ for p in prey: p.move(dt) for p in predators: p.move(dt)
[docs]def update_velocities(prey, predators, dt): """ Update the positions of all animals according to .. math:: \\vec{v}(t+\\Delta t) = \\vec{v}(t) + \\frac{1}{m}\\vec{F} \Delta t Args: prey (list): a list of :class:`Prey` objects predators (list): a list of :class:`Predator` objects force (array): array of forces on all bodies dt (float): time step :math:`\\Delta t` """ for p in prey: p.accelerate(dt) for p in predators: p.accelerate(dt)
[docs]def velocity_verlet(prey, predators, box, force_parameters, dt, time, trajectory_dt = 1.0): """ Verlet algorithm for integrating the equations of motion, i.e., advancing time. There are a few ways to implement Verlet. The leapfrog version works as follows: First, forces are calculated for all animals and velocities are updated by half a time step, :math:`\\vec{v}(t+\\frac{1}{2}\\Delta t) = \\vec{v}(t) + \\frac{1}{2m}\\vec{F} \Delta t`. Then, these steps are repeated: * Positions are updated by a full time step using velocities but not forces, .. math :: \\vec{r}(t+\\Delta t) = \\vec{r}(t) + \\vec{v}(t+\\frac{1}{2}\\Delta t) \Delta t. * Forces are calculated at the new positions, :math:`\\vec{F}(t + \\Delta t)`. * Velocities are updated by a full time step using the forces .. math :: \\vec{v}(t+\\frac{3}{2}\\Delta t) = \\vec{v}(t+\\frac{1}{2}\\Delta t) + \\frac{1}{m}\\vec{F}(t+\\Delta t) \Delta t These operations are done using the methods :meth:`calculate_forces`, :meth:`update_velocities` and :meth:`update_positions_no_force`. Because velocities were updated by half a time step in the beginning of the simulation, positions and velocities are always offset by half a timestep. You always use the one that has advanced further to update the other and this results in a stable algorithm. Args: prey (list): a list of :class:`Prey` objects predators (list): a list of :class:`Predator` objects box (flocking.PeriodicBox): supercell force_parameters (list): force parameters dt (float): time step :math:`\\Delta t` time (float): the total system time to be simulated trajectory_dt (float): the positions of prey are saved at these these time intervals - does not affect the dynamics in any way """ steps = int(time/dt) trajectory_wait = int(trajectory_dt / dt) calculate_forces(prey, predators, box, force_parameters) # get velocities at half time-step update_velocities(prey, predators, 0.5*dt) for i in range(steps): update_positions_no_force(prey, predators, dt) calculate_forces(prey, predators, box, force_parameters) update_velocities(prey, predators, dt) print_progress(i+1,steps) if i%trajectory_wait == 0: for part in prey: part.save_position() for part in predators: part.save_position() # to calculate proper kinetic energy, return to proper timestep update_velocities(prey, predators, -0.5*dt)
[docs]def main(dynamic_file = "flock_dynamics.txt", system_file = "flock_system.txt"): """ The main program. The program reads info from files, runs the simulation, and plots the trajectory. Args: dynamic_file (str): file with force parameters system_file (str): file with simulation parameters """ # read force parameters cohesion_p, alignment_p, separation_p, avoid_p, chase_p = read_dynamics_file(dynamic_file) force_p = (cohesion_p, alignment_p, separation_p, avoid_p, chase_p) # read simulation parameters prey_p, predator_p, sky_p, time_p = read_system_file(system_file) # create simulation box dimension = sky_p[D_INDEX] L = sky_p[SIZE_INDEX] if dimension == 2: box = PeriodicBox( np.array([L, L, 0.1]) ) else: box = PeriodicBox( np.array([L, L, L]) ) # create a random set of animals birds, predators = create_random_flock(prey_p, predator_p, box) dt = time_p[DT_INDEX] simulation_time = time_p[TOTAL_INDEX] recording_dt = time_p[RDT_INDEX] # run simulaiton velocity_verlet(birds, predators, box, force_p, dt, simulation_time, recording_dt) # animate simulation animate( birds, predators, box, [1,1] )
if __name__ == "__main__": random = random.default_rng() main()