import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import random
import math

### Problem 4

4a

def randwalk(n, d):
polymer = np.zeros((n, d))
for i in range(n):
polymer[i] = polymer[i-1]
blah = random.randint(0, d-1)
if random.random() < 0.5:
polymer[i, blah] = polymer[i, blah] - 1
else:
polymer[i,blah] = polymer[i,blah]+1
return polymer
random.seed(100)
polymer = randwalk(100, 3)
fig = plt.figure()
ax = plt.axes(projection = '3d')
ax.plot(polymer[:, 0], polymer[:, 1], polymer[:, 2], color = 'mediumseagreen')
plt.show

This polymer looks relatively elongated

4b

def mean_displ(n, dim, counts):
#Mean should be zero
R_exp = np.zeros((counts, dim))
for i in range(counts):
blah = randwalk(n = n, d = dim)
for j in range(dim):
R_exp[i, j] = blah[-1, j]
return R_exp
#plotting and stuff
random.seed(1000)
values = mean_displ(500, 2, 1000)
#plotting histogram
fig, (ax1) = plt.subplots(nrows=1,ncols=1)
fig.suptitle('Values of <R>')
ax1.set_title('Histogram representation')
ax1.hist(values, density = True)
num_hit = len([x for x in values if x[0] == 0 and x[1] == 0])
print(f'There were {num_hit} instances in our 1000 simulations where we observed R = 0.')
print('However, the average of the distribution appears to be close to zero')
#plotting scatter plot
#ax2.set_title('Scatter plot representation')
#ax2.scatter(values[:,0],values[:,1], color = 'green')
#ax2.scatter(values[:,0].mean(),values[:,1].mean(),color='red')

```
There were 0 instances in our 1000 simulations where we observed R = 0.
However, the average of the distribution appears to be close to zero
```

4c

While the average value of R of the simulations seems to be pretty close to zero, the simulations themselves appear to rarely achieve that value. Using the RMS value will give us a non-zero value of R, and will likely give us a value that is closer to what is actually observed in polymer behavior.

4d

We would expect the RMS values to scale by a factor of sqrt(n).

#Ask Siyuan how to print things neatly!!
sim_vals = {}
for i in [100, 500, 1000]:
R_vals = mean_displ(n = i, dim = 3, counts = 1000)
MS_vals = np.sum(np.square(R_vals), axis = 1)
sim_vals[i] = np.sqrt(np.mean(MS_vals))
for a, b in sim_vals.items():
predictions = np.sqrt(a)
calculated = b
print(f'n= {a}, Predicted RMS: {predictions}, Calculated RMS: {calculated}')

```
n= 100, Predicted RMS: 10.0, Calculated RMS: 9.96212828666646
n= 500, Predicted RMS: 22.360679774997898, Calculated RMS: 23.039748262513633
n= 1000, Predicted RMS: 31.622776601683793, Calculated RMS: 31.721380802228644
```

From the above results, it seems that the simulation does roughly scale by the predicted factor of sqrt(n)