diff --git a/Jaslo/sheet8.py b/Jaslo/sheet8.py new file mode 100644 index 0000000000000000000000000000000000000000..1d69fffb8308a5a5063cdff0d1df499f4b767e5f --- /dev/null +++ b/Jaslo/sheet8.py @@ -0,0 +1,199 @@ +from timeit import timeit + +import matplotlib.pyplot as plt +import numpy as np + +from scipy.special import factorial + +from tqdm import tqdm + +# Problem 8.1 +# c) +N_0 = 0 +N_try = int(1e6) +zeta = 10 + +N = np.zeros(N_try + 1, dtype=int) +N[0] = N_0 +for i in range(N_try): + if np.random.random() < 0.5: + if np.random.random() < zeta / (N[i] + 1): + N[i+1] = N[i] + 1 + continue + else: + if np.random.random() < N[i] / zeta: + N[i+1] = N[i] - 1 + continue + + N[i+1] = N[i] + +print(f"<N> = {np.mean(N)} ± {np.var(N)}") + +hist, bins = np.histogram(N, bins=np.max(N)) + +x = bins[:-1] +plt.plot(x, hist / np.max(hist), label="MC", zorder=3) + +x = np.linspace(x[0], x[-1]) +y = zeta**x * np.exp(-zeta) / factorial(x) +plt.plot(x, y / np.max(y), label="Poisson") +plt.xlabel(r"$N$") +plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)") +plt.legend() +plt.savefig("problem8.1c.png") +plt.show() + +# d) +N = np.zeros(N_try + 1, dtype=int) +N[0] = N_0 +for i in range(N_try): + if np.random.random() < 0.5: + if np.random.random() < zeta**3 / ((N[i] + 1) * (N[i] + 2) * (N[i] + 3)): + N[i+1] = N[i] + 3 + continue + else: + if np.random.random() < N[i] * (N[i] + 1) * (N[i] + 2) / zeta**3 and N[i] > 2: + N[i+1] = N[i] - 3 + continue + N[i+1] = N[i] + +print(f"<N> = {np.mean(N)} ± {np.var(N)}") + +hist, bins = np.histogram(N, bins=np.max(N)) + +x = bins[:-1] +plt.plot(x, hist / np.max(hist), label="MC") + +x = np.linspace(x[0], x[-1]) +y = zeta**x * np.exp(-zeta) / factorial(x) +plt.plot(x, y / np.max(y), label="Poisson") +plt.xlabel(r"$N$") +plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)") +plt.legend() +plt.savefig("problem8.1d-3-10.png") +plt.show() + + +zeta = 100 +N = np.zeros(N_try + 1, dtype=int) +N[0] = N_0 +for i in range(N_try): + if np.random.random() < 0.5: + if np.random.random() < zeta**3 / ((N[i] + 1) * (N[i] + 2) * (N[i] + 3)): + N[i+1] = N[i] + 3 + continue + else: + if np.random.random() < N[i] * (N[i] + 1) * (N[i] + 2) / zeta**3 and N[i] > 2: + N[i+1] = N[i] - 3 + continue + N[i+1] = N[i] + +print(f"<N> = {np.mean(N)} ± {np.var(N)}") + +hist, bins = np.histogram(N, bins=np.max(N)) + +x = bins[:-1] +plt.plot(x, hist / np.max(hist), label="MC") + +x = np.linspace(x[0], x[-1]) +y = zeta**x * np.exp(-zeta) / factorial(x) +plt.plot(x, y / np.max(y), label="Poisson") +plt.xlabel(r"$N$") +plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)") +plt.legend() +plt.savefig("problem8.1d-3-100.png") +plt.show() + + +zeta = 10 +M = 3 +N = np.zeros(N_try + 1, dtype=int) +N[0] = N_0 +for i in range(N_try): + m = np.random.randint(1, M+1) + if np.random.random() < 0.5: + if np.random.random() < zeta**m / np.multiply.accumulate(np.arange(N[i] + 1, N[i] + m + 1))[-1]: + N[i+1] = N[i] + m + continue + else: + if np.random.random() < np.multiply.accumulate(np.arange(N[i], N[i] + m))[-1] / zeta**m and N[i] >= m: + N[i+1] = N[i] - m + continue + N[i+1] = N[i] + +print(f"<N> = {np.mean(N)} ± {np.var(N)}") + +hist, bins = np.histogram(N, bins=np.max(N)) + +x = bins[:-1] +plt.plot(x, hist / np.max(hist), label="MC") + +x = np.linspace(x[0], x[-1]) +y = zeta**x * np.exp(-zeta) / factorial(x) +plt.plot(x, y / np.max(y), label="Poisson") +plt.xlabel(r"$N$") +plt.ylabel(r"$p^\textrm{eq}(N)$ (normalised)") +plt.legend() +plt.savefig("problem8.1d-m-100.png") +plt.show() + +# e) +zeta = 10 +N = np.zeros(N_try + 1, dtype=int) +N[0] = N_0 +for i in range(N_try): + if np.random.random() < 0.5: + if np.random.random() < zeta / (N[i] + 1): + N[i+1] = N[i] + 1 + continue + else: + if np.random.random() < N[i] / zeta: + N[i+1] = N[i] - 1 + continue + + N[i+1] = N[i] + +k = np.arange(100, N_try - 20, 20, dtype=int) +C = np.zeros((len(k), 21)) +l = np.arange(21, dtype=int) +for i in tqdm(l): + for j in range(len(k)): + C[j,i] = np.mean(N[k[j] + i] * N[k[j]]) - np.mean(N)**2 + +plt.plot(l, np.mean(C, axis=0)) +plt.xlabel(r"$l$") +plt.ylabel(r"$C(l)$") +plt.savefig("problem8.1e.png") +plt.show() + +# f) +def mc(): + zeta = 1000 + M = 3 + N = np.zeros(N_try + 1, dtype=int) + N[0] = N_0 + for i in range(N_try): + m = np.random.randint(1, M+1) + if np.random.random() < 0.5: + if np.random.random() < zeta**m / np.multiply.accumulate(np.arange(N[i] + 1, N[i] + m + 1))[-1]: + N[i+1] = N[i] + m + continue + else: + if np.random.random() < np.multiply.accumulate(np.arange(N[i], N[i] + m))[-1] / zeta**m and N[i] >= m: + N[i+1] = N[i] - m + continue + N[i+1] = N[i] + +def poisson(): + zeta = 1000 + l = np.exp(-zeta) + for _ in range(N_try): + p = 1 + k = 0 + + while True: + p *= np.random.random() + k += 1 + + if p <= l: + break