Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from dataclasses import dataclass
- import numpy as np
- import matplotlib.pyplot as plt
- @dataclass
- class TimeGrid:
- t_min: float = 1e-3
- t_max: float = 1e4
- dt: float = 0.002
- def values(self) -> np.ndarray:
- span = self.t_max - self.t_min
- num = int(np.floor(span / self.dt)) + 1
- return np.logspace(np.log10(self.t_min), np.log10(self.t_max), num)
- def msd_free_times(D: float, N: int, t: np.ndarray, seed: int | None = 0):
- rng = np.random.default_rng(seed)
- x = np.zeros(N)
- msd = np.zeros_like(t)
- t_prev = 0.0
- for k, tk in enumerate(t):
- dt = tk - t_prev
- x += np.sqrt(D * dt) * rng.standard_normal(N)
- msd[k] = np.mean(x * x) - np.mean(x) ** 2
- t_prev = tk
- return t, msd
- @dataclass
- class CompositeParams:
- D: float = 1.0
- gamma0: float = 1.0
- gamma: float = 1.0
- k: float = 1.0
- M: int = 0
- N: int = 1000
- seed: int | None = 0
- def msd_composite_times(params: CompositeParams, t: np.ndarray):
- rng = np.random.default_rng(params.seed)
- N, M = params.N, params.M
- x = np.zeros(N)
- sat = np.zeros((N, M)) if M > 0 else np.zeros((N, 0))
- msd = np.zeros_like(t)
- t_prev = 0.0
- for k, tk in enumerate(t):
- dt = tk - t_prev
- if M > 0:
- diff = x[:, None] - sat
- force_on_x = -params.k * np.sum(diff, axis=1)
- x += (dt / params.gamma0) * force_on_x + (np.sqrt(params.D * dt) / params.gamma0) * rng.standard_normal(N)
- sat += (dt / params.gamma) * (params.k * diff) + (np.sqrt(params.D * dt) / params.gamma) * rng.standard_normal((N, M))
- else:
- x += (np.sqrt(params.D * dt) / params.gamma0) * rng.standard_normal(N)
- msd[k] = np.mean(x * x) - np.mean(x) ** 2
- t_prev = tk
- return t, msd
- def style_log_axes(ax: plt.Axes, xlabel: str = "time", ylabel: str = r"⟨(δx)²⟩") -> None:
- ax.set_xscale("log")
- ax.set_yscale("log")
- ax.set_xlabel(xlabel)
- ax.set_ylabel(ylabel)
- ax.grid(True, which="both", ls=":", lw=0.7, alpha=0.6)
- ax.set_xlim(1e-3, 1e4)
- ax.set_ylim(1e-4, 1e4)
- def show_fig37() -> None:
- grid = TimeGrid(dt=0.002)
- t_line = grid.values()
- t1, msd1 = msd_free_times(D=0.1, N=200, t=t_line, seed=1)
- _, msd2 = msd_free_times(D=1.0, N=200, t=t_line, seed=2)
- idx = np.linspace(0, len(t_line) - 1, 60).astype(int)
- t_mark = t_line[idx]
- msd1_mark = msd1[idx]
- msd2_mark = msd2[idx]
- fig, ax = plt.subplots(figsize=(8, 6), dpi=140)
- style_log_axes(ax)
- ax.plot(t_line, msd1, '-', lw=1.5, color='royalblue', label='D = 0.1')
- ax.plot(t_line, msd2, '-', lw=1.5, color='red', label='D = 1.0')
- ax.plot(t_mark, msd1_mark, 'o', ms=3, color='royalblue')
- ax.plot(t_mark, msd2_mark, 'o', ms=3, color='red')
- tt = np.logspace(-3, 4, 50)
- ax.plot(tt, tt, color='0.3', lw=1.0)
- ax.plot(tt, 0.1 * tt, color='0.3', lw=1.0)
- ax.text(3, 2, 't', fontsize=9)
- ax.text(3, 0.2, '0.1t', fontsize=9)
- ax.legend(frameon=False, loc='upper left')
- ax.set_title('Середньоквадратичне зміщення: вільна частинка')
- plt.tight_layout()
- plt.show()
- def show_fig() -> None:
- grid = TimeGrid(dt=0.002)
- base = CompositeParams(D=1.0, gamma0=1.0, gamma=1.0, k=1.0, N=1000)
- p0 = CompositeParams(**{**base.__dict__, 'M': 0, 'seed': 3})
- p10 = CompositeParams(**{**base.__dict__, 'M': 2, 'seed': 4})
- t_line = grid.values()
- t0, msd0 = msd_composite_times(p0, t_line)
- _, msd10 = msd_composite_times(p10, t_line)
- idx = np.linspace(0, len(t_line) - 1, 60).astype(int)
- t_mark = t_line[idx]
- msd0_mark = msd0[idx]
- msd10_mark = msd10[idx]
- fig, ax = plt.subplots(figsize=(8, 6), dpi=140)
- style_log_axes(ax)
- ax.plot(t_line, msd0, '-', lw=1.5, color='crimson', label='M = 0, N = 1000')
- ax.plot(t_line, msd10, '-', lw=1.5, color='darkgreen', label='M = 2, N = 1000')
- ax.scatter(t_mark, msd0_mark, s=5, color='crimson')
- ax.scatter(t_mark, msd10_mark, s=5,color='darkgreen')
- tt = np.logspace(-3, 4, 50)
- ax.plot(tt, tt, color='0.3', lw=1.0)
- ax.set_title('MSD композитної частинки')
- ax.legend(frameon=False, loc='upper left', title='$\\gamma_0 = 1$\n$\\gamma = 1$\n$k = 1$')
- plt.tight_layout()
- plt.show()
- def show_fig41() -> None:
- grid = TimeGrid(dt=0.002)
- base = CompositeParams(D=25.0, gamma0=5.0, gamma=1.0, k=1.0, N=1000)
- p0 = CompositeParams(**{**base.__dict__, 'M': 0, 'seed': 3})
- p10 = CompositeParams(**{**base.__dict__, 'M': 10, 'seed': 4})
- t_line = grid.values()
- t0, msd0 = msd_composite_times(p0, t_line)
- _, msd10 = msd_composite_times(p10, t_line)
- idx = np.linspace(0, len(t_line) - 1, 60).astype(int)
- t_mark = t_line[idx]
- msd0_mark = msd0[idx]
- msd10_mark = msd10[idx]
- fig, ax = plt.subplots(figsize=(8, 6), dpi=140)
- style_log_axes(ax)
- ax.plot(t_line, msd0, '-', lw=1.5, color='crimson', label='M = 0, N = 1000')
- ax.plot(t_line, msd10, '-', lw=1.5, color='darkgreen', label='M = 10, N = 1000')
- ax.scatter(t_mark, msd0_mark, s=5, color='crimson')
- ax.scatter(t_mark, msd10_mark, s=5,color='darkgreen')
- tt = np.logspace(-3, 4, 50)
- ax.plot(tt, tt, color='0.3', lw=1.0)
- ax.set_title('MSD композитної частинки')
- ax.legend(frameon=False, loc='upper left', title='$\\gamma_0 = 5$\n$\\gamma = 1$\n$k = 1$')
- plt.tight_layout()
- plt.show()
- def main() -> None:
- show_fig37()
- show_fig()
- show_fig41()
- if __name__ == "__main__":
- main()
Advertisement
Add Comment
Please, Sign In to add comment
