# Workshop 7 optional material: Simulating dynamical systems


Physical systems that evolve in time can be easily (but not always quickly) simulated on a computer.   All you need to know is the function that takes the system from its current state to a future state:

$S_2 = f(S_1, dt)$

In the above equation, $S_1$ and $S_2$ represent the current and future states, respectively.  $S_2$ is the state of the system after a short amount of time, $dt$.  The function $f$, which encodes the dynamics of the system, often depends on the current state and the choice of $dt$.  Choosing a smaller time step, $dt$, generally improves the accuracy of any simulation.

For example, the position of a particle, $x(t)$, in one dimension can be written as a function of its velocity $v$:

$\frac{dx}{dt} = v$

You could find $x(t)$ using differential equation solving methods, or you could just simulate the motion of the particle.  To do this, we can express the dynamics as:

\begin{align}
\ x_2 & = f(x_1, dt) \\
\ & = x_1 + v*dt \\
\end{align}


```python
import matplotlib.pyplot as plt

velocity = 1.0
dt = 0.01

def update(x, t):
    '''This is like the function f above; you should
       think of the "state" of the system as its position
       at some time.  We update its state according to the
       dynamics of the system, and I prefer to also update
       the value of time in this same function--but this could
       also be done in the loop below.'''
    
    # update the state and time
    x = x + velocity*dt
    t = t + dt
    
    # return the updated values
    return x, t
```


```python
x = 5 # initial x position
t = 0 # initial time

# Some lists to store the positions and times
x_values = [x]  # I like putting the initial values in the list
t_values = [t]

# This is where the actual "simulation" occurs; the update
# function is used in a loop with however many iterations we want
for i in range(100):
    
    x, t = update(x, t)
    x_values.append(x)
    t_values.append(t)
    
# Finally, some plotting to visualize the results
plt.plot(x_values, t_values)
plt.ylabel('position')
plt.xlabel('time')
plt.show()
```

Wow! It's a straight line.  You probably could've guessed that or used some fancy, built-in ODE solving libraries.  But the point is, if you can write down an `update( )` function for some system, then you can simulate it--whether or not there's an analytical way to find a solution.  

There's a chance that your capstone project is just a variant of this problem, so I don't want to steal away too many project ideas in this notebook.  Once you start simulating a physical system, it's good to consider the accuracy of your simulated "solution".  If you don't have an exact solution to compare against, you can adjust the chosen value of $dt$ and see how it changes the simulation.

## Exercise 1 

Change the code above to simulate a particle under constant acceleration in one dimension.  The `update` function should now also update and return the velocity of the particle.  Plot the position of the particle as a function of time using at least two different choices of time step $dt$.  It may be helpful to include $dt$ as an input parameter to the `update` function so you can easily adjust $dt$ in your simulations.


```python

```

Stochastic systems (systems with some randomness) are another great use-case for this method.  [Brownian motion](https://en.wikipedia.org/wiki/Brownian_motion) describes the random motion of a particle (like a piece of dust) in a fluid.  The thermal motion of the fluid molecule bumps the dust particle around in a seemingly random, but well-described, process.  (In theory you could also model this by keeping track of the motion of many tiny fluid particles colliding with each other, the walls of some container, and the large dust particle.)  

Below, I've written down an update function for the position of a particle undergoing Brownian motion.


```python
import numpy as np
import matplotlib.pyplot as plt

def update(x, t, dt):
    
    x = x + np.sqrt(dt)*np.random.randn()
    t = t + dt
    
    return x, t
```


```python
x = 0
t = 0
x_values = [x]
t_values = [t]

for i in range(100):
    
    x, t = update(x, t, 0.1)
    
    x_values.append(x)
    t_values.append(t)
    
plt.plot(t_values, x_values)
plt.show()
```

## Exercise 2

Repeat the Brownian motion simulation many (~1000) times and plot the distribution of final positions.  The mean of the distribution should be zero, but there is some variance to this distribution.  How does the variance and/or standard deviation of this distribution scale with time? 


```python

```

## Viral spreading

(Very relevant in these times of covid...)

For a more complicated example, here's a simulation of viral spreading following the mathematics of [this paper](https://www.pnas.org/content/pnas/111/46/E4911.full.pdf). You can play with the parameter $\mu$ (mu) to see how it controls the nature of the spreading.


```python
import numpy as np

prob = 0.4  # Probability of infection
L = 2000  # Size of grid

# The distribution from which "jumps" are drawn
def jump(y, mu, L, C):
    return np.power((y*(np.power(L, -mu) - np.power(C, -mu)) + np.power(C, -mu)), -1/mu)

# function that updates the state of the population and list of infected individuals
def update(population, infected,  mu, L, C):
    for inf in infected:           
        if np.random.random() > prob:
            start_x = inf[0] 
            start_y = inf[1] 
            angle = np.random.uniform(0, 2*np.pi)
            dist = jump(np.random.random(), mu, L, C)
            end_x = int(start_x + np.cos(angle)*dist) % L
            end_y = int(start_y + np.sin(angle)*dist) % L
            if population[end_x, end_y] == 0:
                population[end_x, end_y] = 1
                infected.append([end_x, end_y])
    
    return population, infected

```


```python
population = np.zeros((int(L),int(L)))
population[1000,1000] = 1

infected = [[1000,1000]]
  
mu = 1.8
C = 1.5
pending = []
    
ts = []
rs = []
populations = np.zeros((30, int(L), int(L)))

# The simulation -- each iteration is more computationally expensive than the last,
# so be careful changing the number of iterations
for i in range(26):
    populations[i] = population
    population, infected = update(population, infected, mu, L, C)
```


```python
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm,rc
from matplotlib.colors import ListedColormap, LinearSegmentedColormap

dark = cm.get_cmap('Dark2_r', 256)
newcolors = dark(np.linspace(0, 1, 256))
white = np.array([1,1,1,1])
newcolors[:25, :] = white
newcmp = ListedColormap(newcolors)



fig = plt.figure(figsize=(8,8))

a = populations[25]  # You can change this number to get an earlier/later state of the population

im = plt.imshow(a, interpolation='none', aspect='auto', vmin=0, vmax=1, cmap = newcmp)
plt.show()
```


```python

```

096306665433d24b8bcda30006b1a8a632137abe
ipynb
Jupyter Notebook
Week09/WS07/Workshop07_optional.ipynb
ds-connectors/Physics-88-SP22
5f3a74cbd50063fd09b91c28da534176968bb3f2
BSD-3-Clause
Week09/WS07/Workshop07_optional.ipynb
ds-connectors/Physics-88-SP22
5f3a74cbd50063fd09b91c28da534176968bb3f2
BSD-3-Clause
Week09/WS07/Workshop07_optional.ipynb
ds-connectors/Physics-88-SP22
5f3a74cbd50063fd09b91c28da534176968bb3f2
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES