Testing Stochastic Motion via the Fluctuation-Dissipation Theorem

Theory

Stochastic motion on a particle level is often well described by Brownian3D or Langevin3D dynamics [1,2]. The Fluctuation-Dissipation theorem [3] gives estimates the for mean squared displacement (MSD) of a collection of such particles and can thus be used to test the numerical implementation.

Brownian Dynamics

In the case of Brownian3D dynamics

$$\begin{equation} \dot{X} = -\frac{D}{k_B T} \nabla V(X) + \sqrt{2D}R(t), \end{equation}$$

Einstein predicted that the diffusion constant is proportional to the particles mobility [4]. The MSD can be derived by applying the FDT to the probability density function (PDF) of a brownian particle

$$\begin{equation} P(x,t) = \frac{1}{\sqrt{4\pi D t}}\exp\left(-\frac{(x-x_0)^2}{4Dt}\right). \end{equation}$$

By calculating the fourier transform, we obtain the characteristic function

$$\begin{equation} G(k) = \int \text{e}^{ikx} P(x,t|x0)dx = \text{exp}(ikx_0 - k^2Dt). \end{equation}$$

We obtain the moments by differentiating the characteristic function $\kappa_n = (-i)^n\partial_k^n G(k)|_{k=0}$.

$$\begin{align} \kappa_1 &= x_0\\ \kappa_2 &= 2Dt \end{align}$$

For more than one spatial dimension, this approach can be generalized and we obtain

$$\begin{equation} \left<r^2(t)\right> = 2d D t \end{equation}$$

where $d$ is the number of spatial dimensions of our system.

Langevin Dynamics

The estimate of the MSD for Langevin3D dynamics

$$\begin{equation} M \ddot{X} = - \nabla V(X) - \lambda M \dot{X} + \sqrt{2M\lambda k_B T}R(t) \end{equation}$$

is given by [5]

$$\begin{equation} \left<r^2(t)\right> = v^2(0) \frac{1-\text{e}^{-\lambda t}}{\lambda^2} - \frac{d k_B T}{m\lambda^2} \left(1-\text{e}^{-\lambda t}\right) \left(3 - \text{e}^{-\lambda t}\right) + \frac{2 d k_B T}{m\lambda}t. \end{equation}$$

It is now a matter of writing code such that the described system can be solved, analyzed and tested.

Code

First, we define methods, to initialize the agents and set up the simulation parameters.

 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37


struct Parameters {
    n_particles: usize,
    domain_size: f64,

    dt: f64,
    n_steps: u64,
    save_interval: u64,

    diffusion_constant: f64,

    storage_name: std::path::PathBuf,
    random_seed: u64,
}

impl Default for Parameters {
    fn default() -> Self {
        Self {
            n_particles: 160,
            domain_size: 200.0,

            dt: 1e-3,
            n_steps: 200,
            save_interval: 1,

            diffusion_constant: 1.0,
            storage_name: "out/brownian".into(),
            random_seed: 0,
        }
    }
}

We specifically require that results are only stored on the disk (in json format) if #[cfg(not(debug_assertions))]. This is usually the case when running in release mode.

Define settings

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69


fn define_settings(
    parameters: &Parameters,
) -> Result<
    cellular_raza_core::backend::chili::Settings<FixedStepsize<f64>, true>,
    Box<dyn std::error::Error>,
> {
    let time = cellular_raza_core::time::FixedStepsize::from_partial_save_steps(
        0.0,
        parameters.dt,
        parameters.n_steps,
        parameters.save_interval,
    )?;

    // Only save results when debug_assertions are not active.
    let location = parameters.storage_name.to_owned();
    let storage = cellular_raza_core::storage::StorageBuilder::new()
        .priority([
            cellular_raza_core::storage::StorageOption::Memory,
            #[cfg(not(debug_assertions))]
            cellular_raza_core::storage::StorageOption::SerdeJson,
        ])
        .location(location)
        .init();

    Ok(cellular_raza_core::backend::chili::Settings {
        n_threads: 1.try_into().unwrap(),
        time,
        storage,
        show_progressbar: false,
    })
}

Calculate mean and standard deviation

 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107


fn calculate_mean_std_dev<const D: usize>(
    parameters: &Parameters,
    positions: impl IntoIterator<Item = (u64, Vec<SVector<f64, D>>)>,
) -> Result<std::collections::HashMap<u64, (f64, f64)>, Box<dyn std::error::Error>> {
    let initial_position = SVector::from([parameters.domain_size / 2.0; D]);
    let square_displacements = positions
        .into_iter()
        .map(|(iteration, positions)| {
            (
                iteration,
                positions
                    .into_iter()
                    .map(|pos| (pos - initial_position).norm_squared())
                    .collect::<Vec<_>>(),
            )
        })
        .collect::<std::collections::HashMap<u64, Vec<f64>>>();

    let means = square_displacements
        .iter()
        .map(|(iteration, displs)| (*iteration, displs.iter().sum::<f64>() / displs.len() as f64))
        .collect::<std::collections::HashMap<_, _>>();

    let means_std_dev = square_displacements
        .iter()
        .map(|(iteration, displs)| {
            let sigma = (displs
                .iter()
                .map(|displ| (displ - means[&iteration]).powf(2.0))
                .sum::<f64>()
                / displs.len() as f64)
                .sqrt();
            (*iteration, (means[&iteration], sigma))
        })
        .collect::<std::collections::HashMap<_, _>>();
    Ok(means_std_dev)
}

In order to test our numerical results, we compare the mean squared displacement at every iteration with the analytical estimate. Their difference should be within the margin of multiple $\sigma$ of each other. This method of testing will be valid for the chosen set of parameters. Should a future code refactor (which does not affect the tests) produce different results, this test will only fail when results are not in the specified confidence region. For our chosen threshold of $3.5\sigma$, the probability for any result to be outside this band is $p=0.0465\%$.

Brownian

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123


fn analyze_positions_brownian<const D: usize>(
    parameters: &Parameters,
    positions: impl IntoIterator<Item = (u64, Vec<SVector<f64, D>>)>,
) -> Result<(), Box<dyn std::error::Error>> {
    let means_std_dev = calculate_mean_std_dev(&parameters, positions)?;

    // Calculate probability for values to be identical
    for &iteration in means_std_dev.keys() {
        let expected =
            2.0 * D as f64 * parameters.diffusion_constant * iteration as f64 * parameters.dt;
        let (mean, std_dev) = means_std_dev[&iteration];
        assert!((expected - mean).abs() <= 3.5 * std_dev);
    }
    Ok(())
}

We now combine all functionalities in the test_brownian! macro. Every agent, is positioned in the middle of the simulation domain with no initial velocity (irrelevant for the Brownian case). We have chosen the domain size such that no agent reaches the outer boundary which change the MSD.

Combine previously defined functions with macro

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168


macro_rules! test_brownian {
    ($parameters:ident, $domain_name:ident, $particle_name:ident, $d:literal) => {{
        let domain_size = $parameters.domain_size;
        assert!(domain_size > 0.0);
        let mut domain =
            $domain_name::from_boundaries_and_n_voxels([0.0; $d], [domain_size; $d], [3; $d])?;
        domain.rng_seed = $parameters.random_seed;

        let initial_position = nalgebra::SVector::from([domain_size / 2.0; $d]);
        let particles = (0..$parameters.n_particles)
            .map(|_| $particle_name::new(
                initial_position.into(),
                $parameters.diffusion_constant,
                1.0
            ));

        let settings = define_settings(&$parameters)?;

        let storage_access = cellular_raza_core::backend::chili::run_simulation!(
            agents: particles,
            domain: domain,
            settings: settings,
            aspects: [Mechanics]
        )?;

        let positions = storage_access
            .cells
            .load_all_elements()?
            .into_iter()
            .map(|(iteration, particles)| {
                (
                    iteration,
                    particles
                        .into_iter()
                        .map(|(_, (p, _))| cellular_raza::concepts::Position::pos(&p.cell))
                        .collect::<Vec<_>>(),
                )
            })
            .collect::<std::collections::HashMap<u64, _>>();

        analyze_positions_brownian(&$parameters, positions)?;
        Result::<(), Box<dyn std::error::Error>>::Ok(())
    }}
}

This macro can afterwards be used in a collection of multiple functions to test various parameter configurations.

170
171
172
173
174
175
176
177
178


#[test]
fn diffusion_constant_brownian_3d_1() -> Result<(), Box<dyn std::error::Error>> {
    let mut parameters = Parameters::default();
    parameters.storage_name = "out/brownian_3d_1".into();
    parameters.diffusion_constant = 1.0;
    parameters.random_seed = 1;
    test_brownian!(parameters, CartesianCuboid3New, Brownian3D, 3)?;
    Ok(())
}

Langevin

Since the langevin equation introduces fluctuations for the velocity of the particle, the first integration step does not change the position. Thus we skip the initial step when comparing results.

260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288


fn analyze_positions_langevin<const D: usize>(
    parameters: &Parameters,
    positions: impl IntoIterator<Item = (u64, Vec<SVector<f64, D>>)>,
    damping: f64,
) -> Result<(), Box<dyn std::error::Error>> {
    let means_std_dev = calculate_mean_std_dev(&parameters, positions)?;

    // Calculate probability for values to be identical
    let mass = 1.0;
    let kb_temperature = kb_temp(mass, parameters.diffusion_constant, damping);
    for &iteration in means_std_dev.keys() {
        // Here we assume that the initial velocity v(0) = 0
        // https://en.wikipedia.org/wiki/Langevin_equation#Trajectories_of_free_Brownian_particles
        if iteration > 1 {
            let (mean, std_dev) = means_std_dev[&iteration];
            let time = iteration as f64 * parameters.dt;
            let expected = -(D as f64) * kb_temperature / mass / damping.powf(2.0)
                * (1.0 - (-damping * time).exp())
                * (3.0 - (-damping * time).exp())
                + 2.0 * D as f64 * kb_temperature / mass * time / damping;
            assert!((expected - mean).abs() <= 3.5 * std_dev);
        }
    }
    Ok(())
}

fn kb_temp(mass: f64, diffusion_constant: f64, damping: f64) -> f64 {
    diffusion_constant * damping * mass
}

Macro to compare Langevin Results

299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349


        assert!(domain_size > 0.0);
        let mut domain =
            $domain_name::from_boundaries_and_n_voxels([0.0; $d], [domain_size; $d], [3; $d])?;
        domain.rng_seed = $parameters.random_seed;
        let initial_position = nalgebra::SVector::from([domain_size / 2.0; $d]);
        let mass = 1.0;
        let kb_temperature = $parameters.diffusion_constant * $damping * mass;
        let particles = (0..$parameters.n_particles)
            .map(|_| $particle_name ::new(
                initial_position.into(),
                [0.0; $d].into(),
                mass,
                $damping,
                kb_temperature,
            ));
        let settings = define_settings(&$parameters)?;

        let storage_access = cellular_raza_core::backend::chili::run_simulation!(
            agents: particles,
            domain: domain,
            settings: settings,
            aspects: [Mechanics],
        )?;

        let positions = storage_access
            .cells
            .load_all_elements()?
            .into_iter()
            .map(|(iteration, particles)| {
                (
                    iteration,
                    particles
                        .into_iter()
                        .map(|(_, (p, _))| cellular_raza::concepts::Position::pos(&p.cell))
                        .collect::<Vec<_>>(),
                )
            })
            .collect::<std::collections::HashMap<u64, _>>();

        analyze_positions_langevin(&$parameters, positions, $damping)?;
        Result::<(), Box<dyn std::error::Error>>::Ok(())
    }}
}

#[test]
fn diffusion_constant_langevin_3d_1() -> Result<(), Box<dyn std::error::Error>> {
    let mut parameters = Parameters::default();
    parameters.storage_name = "out/langevin_3d_1".into();
    parameters.diffusion_constant = 80.0;
    parameters.random_seed = 1;
    test_langevin!(parameters, CartesianCuboid3New, Langevin3D, 3, 10.0)?;

Example Test For Langevin Results

351
352
353
354
355
356
357
358


}

#[test]
fn diffusion_constant_langevin_3d_2() -> Result<(), Box<dyn std::error::Error>> {
    let mut parameters = Parameters::default();
    parameters.storage_name = "out/langevin_3d_2".into();
    parameters.diffusion_constant = 40.0;
    parameters.random_seed = 2;

Running the Tests

To run the supplied test navigate to the cellular_raza folder containing the bundle of all dependencies and execute the tests in release mode.

git clone https://github.com/jonaspleyer/cellular_raza
cd cellular_raza/cellular_raza
cargo test -r

# For visualization
python tests/plot_brownian_langevin.py

This will generate results stored in json format in the out directory. To visualize the results, use the provided python script. It requires matplotlib,numpy and scipy as dependencies.

Script to plot results

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249


import numpy as np
import matplotlib.pyplot as plt
import json
from glob import glob
import scipy as sp

def get_last_save_dir(storage_name: str) -> str:
    return list(sorted(glob("out/{}/*".format(storage_name))))[-1]

def get_trajectories(storage_name: str) -> np.ndarray:
    last_save_dir = get_last_save_dir(storage_name)

    # Obtain all values for cells
    iterations_cells = []
    for iteration_dir in sorted(glob(last_save_dir + "/cells/json/*")):
        cells = []
        for batch in list(sorted(glob(iteration_dir + "/*"))):
            with open(batch) as f:
                cells.extend(json.load(f)["data"])
        iterations_cells.append(cells)

    # Calculate the trajectories
    trajectories = np.array(
        np.array([[
            values_at_iter[j]["element"][0]["cell"]["pos"]
            for values_at_iter in iterations_cells
        ] for j in range(len(iterations_cells[0]))]
    ))
    return trajectories

def get_domain_boundaries(storage_name: str) -> tuple[np.ndarray, np.ndarray]:
    last_save_dir = get_last_save_dir(storage_name)
    # Obtain all values for subdomains
    iteration_dir = glob(last_save_dir + "/subdomains/json/*")[0]
    single = glob(iteration_dir + "/*")[0]
    with open(single) as f:
        subdomain = json.load(f)["element"]
        dmin = np.array([subdomain["domain_min"][0]])
        dmax = np.array([subdomain["domain_max"][0]])
        # return np.ndarray([subdomain["domain_min"]]), np.ndarray([subdomain["domain_max"]])
        return dmin, dmax

def plot_2d_only(trajectories: np.ndarray, domain_middle: np.ndarray, last_save_dir: str):
    # Plot the obtained results for each iteration
    dh = np.max(np.abs(trajectories - domain_middle), axis=(0,1))
    s = dh[0] / dh[1]
    lim_lower = domain_middle - 1.1 * dh
    lim_upper = domain_middle + 1.1 * dh
    xlim = [lim_lower[0], lim_upper[0]]
    ylim = [lim_lower[1], lim_upper[1]]

    fig, ax = plt.subplots(figsize=(8, s*8))
    ax.set_title("Trajectories")
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    for traj in trajectories:
        ax.plot(traj[:,0], traj[:,1], color="k", linestyle="-")
    fig.tight_layout()
    fig.savefig("{}/trajectories.png".format(last_save_dir))

    # Plot a heatmap of the total explored space
    heatmap, _, _ = np.histogram2d(
        trajectories[:,:,0].reshape((-1,)),
        trajectories[:,:,1].reshape((-1,)),
        range=[xlim, ylim],
        bins=50
    )
    extent = [*lim_lower, *lim_upper]
    fig, ax = plt.subplots(figsize=(8, 8*s))
    ax.imshow(heatmap.T, extent=extent, origin='lower')
    ax.set_title("Heatmap of explored space")
    fig.tight_layout()
    fig.savefig("{}/heatmap.png".format(last_save_dir))
    plt.close(fig)

def plot_msd(trajectories: np.ndarray, domain_middle: np.ndarray):
    # Plot the mean squared displacement per iteration
    msd = np.mean(np.sum((trajectories - domain_middle)**2, axis=2), axis=0)
    msd_err = np.std(np.sum((trajectories - domain_middle)**2, axis=2), axis=0)\
        / trajectories.shape[0]**0.5

    x = np.arange(msd.shape[0])
    fig, ax = plt.subplots()
    ax.errorbar(x, msd, msd_err, color="gray", linestyle="-", label="Mean Displacements")
    return fig, ax, x, msd, msd_err

def plot_brownian(
        storage_name: str,
        diffusion_constant: float,
        dimension: int,
        dt: float,
    ):
    print(storage_name)

    # Get trajectories
    last_save_dir = get_last_save_dir(storage_name)
    trajectories = get_trajectories(storage_name)

    # Get Domain size
    domain_min, domain_max = get_domain_boundaries(storage_name)
    domain_middle = 0.5 * (domain_min + domain_max)

    fig, ax, x, msd, msd_err = plot_msd(trajectories, domain_middle)

    def prediction_brownian(t, dim, diffusion):
        return 2 * dim * diffusion * t

    y = prediction_brownian(dt * x, dimension, diffusion_constant)
    popt, pcov = sp.optimize.curve_fit(
        lambda t, D: prediction_brownian(t, D, dimension),
        dt * x,
        msd,
        sigma=msd_err,
    )

    ax.plot(
        x,
        y,
        label="Prediction $2nDt$ with $D={}$".format(diffusion_constant),
        color="k",
        linestyle=":",
    )
    ax.plot(
        x,
        prediction_brownian(dt * x, popt[0], dimension),
        label="Fit $D={:4.3} \\pm {:4.3}$".format(popt[0], pcov[0][0]**0.5),
        linestyle="--",
        color="orange",
    )

    ax.legend()

    ax.set_title("Mean Squared Displacement {}".format(storage_name))
    fig.tight_layout()
    fig.savefig("{}/mean-squared-displacement.png".format(last_save_dir))
    plt.close(fig)

    if trajectories.shape[2] == 2:
        plot_2d_only(trajectories, domain_middle, last_save_dir)

# Note that the values here do need to match the values of the
# tests as defined by the cellular_raza test suite.
# See cellular_raza/tests/brownian_diffusion_constant_approx.rs
BROWNIAN_VALUES = [
    {"storage_name":"brownian_1d_1", "diffusion_constant": 1.0, "dimension":1, "dt":1e-3},
    {"storage_name":"brownian_1d_1", "diffusion_constant": 1.0, "dimension":1, "dt":1e-3},
    {"storage_name":"brownian_1d_2", "diffusion_constant": 0.5, "dimension":1, "dt":1e-3},
    {"storage_name":"brownian_1d_3", "diffusion_constant":0.25, "dimension":1, "dt":1e-3},
    {"storage_name":"brownian_2d_1", "diffusion_constant": 1.0, "dimension":2, "dt":1e-3},
    {"storage_name":"brownian_2d_2", "diffusion_constant": 0.5, "dimension":2, "dt":1e-3},
    {"storage_name":"brownian_2d_3", "diffusion_constant":0.25, "dimension":2, "dt":1e-3},
    {"storage_name":"brownian_3d_1", "diffusion_constant": 1.0, "dimension":3, "dt":1e-3},
    {"storage_name":"brownian_3d_2", "diffusion_constant": 0.5, "dimension":3, "dt":1e-3},
    {"storage_name":"brownian_3d_3", "diffusion_constant":0.25, "dimension":3, "dt":1e-3},
]

for kwargs in BROWNIAN_VALUES:
    plot_brownian(**kwargs)

def plot_langevin(
        storage_name: str,
        damping: float,
        diffusion: float,
        dim: int,
        dt: float
    ):
    print(storage_name)
    kb_temperature_div_mass = diffusion * damping

    # Get trajectories
    last_save_dir = get_last_save_dir(storage_name)
    trajectories = get_trajectories(storage_name)

    # Get Domain size
    domain_min, domain_max = get_domain_boundaries(storage_name)
    domain_middle = 0.5 * (domain_min + domain_max)

    # Plot the mean squared displacement per iteration
    fig, ax, x, msd, msd_err = plot_msd(trajectories, domain_middle)

    def prediction_langevin(t, damping, kb_temperature_div_mass, dim):
        return - dim * kb_temperature_div_mass / damping**2\
            * (1.0 - np.exp(- damping * t))\
            * (3.0 - np.exp(- damping * t))\
            + 2.0 * dim * kb_temperature_div_mass * t / damping

    popt, pcov = sp.optimize.curve_fit(
        lambda t, damping, kb_temp_div_mass: prediction_langevin(
            t, damping, kb_temp_div_mass, dim
        ),
        dt * x[1:],
        msd[1:],
        # sigma=msd_err[1:],
        p0=(damping, kb_temperature_div_mass),
    )

    y = prediction_langevin(dt * x, damping, kb_temperature_div_mass, dim)
    ax.plot(
        x,
        y,
        label="Prediciton $\\left<r^2(t)\\right>$",
        color="k",
        linestyle=":",
    )
    ax.plot(
        x,
        prediction_langevin(dt * x, *popt, dim),
        label="Fit $\\lambda={:4.3} \\pm {:4.3}, D={:4.3}\\pm {:4.3}$".format(
            popt[0],
            pcov[0][0]**0.5,
            popt[1] / popt[0],
            ((pcov[1][1]/popt[0]**2) + (pcov[0][0]*popt[1]**2/popt[0]**4))**0.5
        ),
        linestyle="--",
        color="orange",
    )

    ax.legend()

    ax.set_title("Mean Squared Displacement")
    fig.tight_layout()
    fig.savefig("{}/mean-squared-displacement.png".format(last_save_dir))

    if trajectories.shape[2] == 2:
        plot_2d_only(trajectories, domain_middle, last_save_dir)

# Note that the values here do need to match the values of the
# tests as defined by the cellular_raza test suite.
# See cellular_raza/tests/brownian_diffusion_constant_approx.rs
LANGEVIN_VALUES = [
    {"storage_name": "langevin_3d_1", "diffusion": 80.0, "dim": 3, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_3d_2", "diffusion": 40.0, "dim": 3, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_3d_3", "diffusion": 20.0, "dim": 3, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_3d_4", "diffusion": 40.0, "dim": 3, "damping":  1.0, "dt": 1e-3},
    {"storage_name": "langevin_3d_5", "diffusion": 40.0, "dim": 3, "damping":  0.1, "dt": 1e-3},
    {"storage_name": "langevin_2d_1", "diffusion": 80.0, "dim": 2, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_2d_2", "diffusion": 40.0, "dim": 2, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_2d_3", "diffusion": 20.0, "dim": 2, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_2d_4", "diffusion": 20.0, "dim": 2, "damping":  1.0, "dt": 1e-3},
    {"storage_name": "langevin_2d_5", "diffusion": 20.0, "dim": 2, "damping":  0.1, "dt": 1e-3},
    {"storage_name": "langevin_1d_1", "diffusion": 80.0, "dim": 1, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_1d_2", "diffusion": 40.0, "dim": 1, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_1d_3", "diffusion": 20.0, "dim": 1, "damping": 10.0, "dt": 1e-3},
    {"storage_name": "langevin_1d_4", "diffusion": 20.0, "dim": 1, "damping":  1.0, "dt": 1e-3},
    {"storage_name": "langevin_1d_5", "diffusion": 20.0, "dim": 1, "damping":  0.1, "dt": 1e-3},
]

for kwargs in LANGEVIN_VALUES:
    plot_langevin(**kwargs)

Results

Brownian Results

Langevin Results

ℹ️

Note that the displayed errorbars are not identical to the ones used for automated testing. We plot the standard error $\sigma / \sqrt{N}$ while the standard deviation $\sigma$ is used to compare numerical results.

The figures above show results for the brownian and langevin cases. We first compare the MSD with its predicted values and obtain parameters by fitting the predictor to our simulated values. In the second plots, all trajectories of all particles are shown. We can clearly see that the motions in the Brownian case are much more fluctuating while the curves of the Langevin case are more smooth although this behaviour depends on the chosen parameters. In the last step, a heatmap is shown which counts the number of times each bin is visited by a particle.

Discussion

In practice, our testing scheme is most sensitive in the early stages of simulation since the deviation of the MSD is lowest at this point in time. This is a practical challenge since the initial steps of the solver are less accurate and thus the overall uncertainty of the numerically obtained results.

References

[1] D. S. Lemons and A. Gythiel, “Paul Langevin’s 1908 paper ‘On the Theory of Brownian Motion’ [‘Sur la théorie du mouvement brownien,’ C. R. Acad. Sci. (Paris) 146, 530–533 (1908)],” American Journal of Physics, vol. 65, no. 11. American Association of Physics Teachers (AAPT), pp. 1079–1081, Nov. 01, 1997. doi: 10.1119/1.18725.

[2] R. Brown, “XXVII. A brief account of microscopical observations made in the months of June, July and August 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies,” The Philosophical Magazine, vol. 4, no. 21. Informa UK Limited, pp. 161–173, Sep. 1828. doi: 10.1080/14786442808674769.

[3] H. B. Callen and T. A. Welton, “Irreversibility and Generalized Noise,” Physical Review, vol. 83, no. 1. American Physical Society (APS), pp. 34–40, Jul. 01, 1951. doi: 10.1103/physrev.83.34.

[4] A. Einstein, “Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen,” Annalen der Physik, vol. 322, no. 8. Wiley, pp. 549–560, Jan. 1905. doi: 10.1002/andp.19053220806.

[5] N. G. VAN KAMPEN, “THE LANGEVIN APPROACH,” Stochastic Processes in Physics and Chemistry. Elsevier, pp. 219–243, 2007. doi: 10.1016/b978-044452965-7/50012-x.