Проблема.

Я пытаюсь создать симуляцию мегаполиса по 2D-модели Изинга. Мой код воспроизводит некоторые из ожидаемых характеристик: т.е. есть критический переход при неопределенно правильной температуре, но:

1. Энергия имеет плавный переход. Намагничивание - это резкий шаг вниз от 1 до 0, в то время как энергия продолжает расти и не показывает никаких изменений градиента.
2. Теплоемкость показывает пик, но вместо увеличения, когда я увеличиваю размер решетки, пик становится меньше.
3. Теплоемкость шумная: пик часто удваивается, и приспособить его к лоренцеву - кошмар. Что еще хуже, при многократном запуске моделирования получаются одни и те же данные (о чем свидетельствует небольшая панель ошибок.

Мой проект - это программа на C++, имитирующая единую решетку, и скрипт на Python, который анализирует результат. На данный момент сценарий передает стандартный вывод программы и назначает его внутренним переменным. Это неэффективно, но для целей отладки.

полный исходный код можно найти здесь. Для удобства ознакомьтесь с основными процедурами программы на C++ и скриптом Python ниже.

Код:

решетка.h

...

/* Converts given pair of indices to those with periodic boundary conditions. */
int inline to_periodic(int row, int col, int size) {
  ...
}

class lattice {
private:
  unsigned int size_;
  std::vector<short> *spins_;
  float J_;
  float H_;

public:
  lattice() noexcept : size_(0), spins_(NULL), J_(1.0), H_(0.0) {}
  lattice(int new_size, double new_J, double new_H) noexcept
      : size_(new_size), spins_(new std::vector<short>(size_ * size_, 1)),
        J_(new_J), H_(new_H) {}
  lattice(const lattice &other) noexcept
      : lattice(other.size_, other.J_, other.H_) {
...
  float compute_point_energy(int row, int col);
};

#endif

simulation.h

#ifndef simulation_h
#define simulation_h

#include "lattice.h"
#include "rng.h"
#include <gsl/gsl_rng.h>

class simulation {
private:
  unsigned int time_ = 0;  // Current time of the simulation.
  rng r_ = rng();
  lattice spin_lattice_;
  double temperature_;
  double mean_magnetisation_ = 1;
  double mean_energy_;
  double total_magnetisation_;
  double total_energy_;
  unsigned int print_interval_ = 1;
  void advance();

public:
  void set_print_interval(unsigned int new_print_interval) { print_interval_ = new_print_interval; }

  simulation(int new_size, double new_temp, double new_J, double new_H)
      : time_(0), spin_lattice_(lattice(new_size, new_J, new_H)), temperature_(new_temp),
        mean_energy_(new_J * (-4)), total_magnetisation_(new_size * new_size),
        total_energy_(compute_energy(spin_lattice_)) {}

  void print_status(FILE *f) {
    f = f==NULL? stdout : f;
    fprintf(f, "%4d\t%e \t%e\t%e\n", time_, mean_magnetisation_,
            mean_energy_, temperature_);
  }
  void advance(unsigned int time_steps, FILE *output);
  double compute_energy(lattice &other);
  double compute_dE(int row, int col) {
    return -2*spin_lattice_.compute_point_energy(row, col);
  }
  void set_to_chequerboard(int step);
  void print_lattice(){
    spin_lattice_.print();
  };
  // void load_custom(const lattice& custom);
};

#endif

simulation.cpp

...

void simulation::advance(unsigned int time_steps, FILE *output) {
  unsigned int area = spin_lattice_.get_size() * spin_lattice_.get_size();
  for (unsigned int i = 0; i < time_steps; i++) {    
    if (time_ % print_interval_ == 0) {
      total_magnetisation_ = spin_lattice_.total_magnetisation();
      mean_magnetisation_ = total_magnetisation_ / area;
      total_energy_ = compute_energy(spin_lattice_);
      mean_energy_ = total_energy_ / area;
      print_status(output);
    }
    advance();
  }
}

void simulation::advance() {
  // #pragma omp parallel for collapse(2)
  for (unsigned int row = 0; row < spin_lattice_.get_size(); row++) {
    for (unsigned int col = 0; col < spin_lattice_.get_size(); col++) {
      double dE = compute_dE(row, col);
      double p = r_.random_uniform();
      if (dE <0 || exp(-dE / temperature_) > p) {
        spin_lattice_.flip(row, col);
      }
    }
  }
  time_++;
}



double simulation::compute_energy(lattice &other) {
  double energy_sum = 0;
  unsigned int max = other.get_size();
  #pragma omp parallel for reduction(+ : energy_sum)
  for (unsigned int i = 0; i < max; i++) {
    for (unsigned int j = 0; j < max; j++) {
      energy_sum += other.compute_point_energy(i, j);
    }
  }
  return energy_sum/2;
}

lattice.cpp

.....

void lattice::print() {
  ...
}

float lattice::compute_point_energy(int row, int col) {
  int accumulator = get(row + 1, col) + get(row - 1, col) + get(row, col - 1) +
                    get(row, col + 1);
  return -get(row, col) * (accumulator * J_ + H_);
}

int lattice::total_magnetisation() {
  ...
}

vestigator.py

.....
# -------------------------------------------------------------------------


def smart_duration(temperature, multiplier=1.):
    return int(((10**3)*base_duration*multiplier)/(
                (temperature - t_c)**2 + breadth))


def investigate_temperature_dependence(temps=None, lattice_sizes=None,
                                       **kwargs):
    if temps is None:
        temps = append(linspace(1.5, 2.2, 8), linspace(2.2, 2.3, 8))
        temps = append(temps, linspace(2.3, 4, 8))
    if lattice_sizes is None:
        lattice_sizes = [50, 250, 300]
    fig, ax = plt.subplots(nrows=2, sharex='col')
    ax[0].set_ylabel('mean magnetisation / arb. u.')
    ax[1].set_ylabel('mean energy / arb. u.')
    ax[0].axvline(x=t_c, color='k', ls='--')
    ax[1].axvline(x=t_c, color='k', ls='--')
    for l in lattice_sizes:
        print(l)
        simulations = [Simulation(lattice_size=l, temperature=t,
                                  duration=smart_duration(l, t)) for t in
                       temps]
        data = [multi_run(s, **kwargs) for s in simulations]
        terminal_magnetisations = array(data).T[0]
        sigma_magnetisations = array(data).T[1]
        terminal_energies = array(data).T[2]
        ax[0].plot(temps, terminal_magnetisations, '-',
                   label='N = ' + str(l))
        ax[1].plot(temps, terminal_energies, '-', label='N =' + str(l))
    plt.xlabel('temperature / arb. u.')
    plt.xticks(
        append(linspace(min(temps), t_c, 3), linspace(t_c, max(temps), 3)))
    save_plot('Critical temperature')


# -------------------------------------------------------------------------


def theoretical_capacity(x, a, b, c, d):
    return c + (b*x)/((x - a)**2 + d)


def investigate_heat_capacity(lattice_sizes=None, temps=None, **kwargs):
    if temps is None:
        temps = linspace(1.5, 3, 60)
    if lattice_sizes is None:
        lattice_sizes = [16,32, 34, 36]

    fig, axes = plt.subplots(nrows=len(lattice_sizes), sharex='col')
    crit_temps = []
    for l, ax in zip(lattice_sizes, axes):
        crit_temps.append(fit_and_plot_capacity(ax, l, temps, **kwargs))

    fig.set_size_inches(10.5, 10.5)
    plt.xlabel('temperature / arb. u.')
    save_plot('Heat capacity')
    return crit_temps


def fit_and_plot_capacity(ax, l, temps, **kwargs):
    """
    Plot the heat capacity of simulations at given temperature and lattice size.
    Afterwards fit a lorentzian and plot.

    Parameters:
    -----------
    ax : pyplot.axis
    what to plot to.

    l : int
    lattice size

    temps: numpy.array
    Array of temperatures where to evaluate heat capacity.

    Returns:
    -------
    popt[0]: float
    most likely critical temperature.
    """
    global use_disk
    use_disk = False
    simulations = [Simulation(lattice_size=l, temperature=t, duration=2)
                   for t in temps]
    # sigmas = [stdev(s.mean_energies[:]) for s in simulations]
    sigmas = array([multi_run(s, **kwargs) for s in simulations]).T[3]
    meta_sigmas = array([multi_run(s, **kwargs) for s in simulations]).T[4]
    # print(meta_sigmas)
    Cs = [sigma**2/(temp**2)*10**3 for temp, sigma in zip(temps, sigmas)]
    C_errs = Cs[:]*meta_sigmas[:]
    try:
        popt, pcov = curve_fit(theoretical_capacity, temps, Cs,
                               sigma=meta_sigmas, bounds=(
                [min(temps) - .2, 0, 0, 0],
                [max(temps) + .2, inf, inf, .7]))
    except RuntimeError:
        popt = [temps[argmax(Cs)], (max(Cs) - min(Cs))/4, min(Cs),
                (max(temps) - min(temps))/4]
        print('I\'m too dumb to fit')

    ax.axvline(x=popt[0], ls='--', color='g')
    ax.plot(temps, theoretical_capacity(temps, *popt), 'g-',
            label='N = ' + str(l) + 'd = ' + str(popt[3]) + ' fit')
    ax.errorbar(temps*10, Cs, fmt='b.', yerr=C_errs, label='N = ' + str(l))
    ax.set_ylabel(r'C $\cdot 10^3$/ arb. u.')
    ax.axvline(x=t_c, ls='-.', color='k')
    ax.set_xticks(
        append(linspace(min(temps), t_c, 6), linspace(t_c, max(temps), 6)))
    ax.legend(loc='best')
    return popt[0]


def multi_run(sim, re_runs:int=2, take_last:int=300):
    """
    Re run the Ising model simulation multiple times, and gather statistics.

    Parameters:
    ----------
    sim: simulation
    re_runs: int
    number of times to repeat simulation
    take_last: int
    How many of the final points to take statistics over.

    Returns:
    list:

    """
    global use_disk
    use_disk = False
    sim.duration = smart_duration(sim.temperature)
    print(sim.duration)
    magnetizations = []
    sigma_magnetizations = []
    energies = []
    sigma_energies = []
    for i in range(re_runs):
        sim.data = sim.run()
        # Make each run take 50% longer, so that we
        # can see if a system is still settling
        sim.duration *= 2

        last_magnetizations = sim.mean_magnetizations()[-take_last:]
        magnetizations.append(abs(mean(last_magnetizations)))
        sigma_magnetizations.append(std(last_magnetizations))

        last_energies = sim.mean_energies()[-take_last:]
        energies.append(mean(last_energies))
        sigma_energies.append(std(last_energies))

    return [mean(magnetizations), mean(sigma_magnetizations),
            mean(energies), mean(sigma_energies),
            std(sigma_energies)/mean(sigma_energies)]


# -------------------------------------------------------------------------

def finite_size_scale(N, t_inf, a, v):
    return t_inf + a*(N**(-1/v))


def investigate_finite_size_scaling(critical_temperatures, lattice_sizes,
                                    **kwargs):
    if critical_temperatures is None:
        critical_temperatures = investigate_heat_capacity(lattice_sizes,
                                                          **kwargs)
    args, cov = curve_fit(finite_size_scale, lattice_sizes,
                                  critical_temperatures)
    plt.plot(lattice_sizes, critical_temperatures, 'b+', label='data')
    plt.plot(lattice_sizes, finite_size_scale(lattice_sizes, *args), 'r-',
             label='fit')
    plt.ylabel('critical temperature / arb. u.')
    plt.xlabel('Lattice size')
    save_plot('Finite size scaling')
    return args[0], sqrt(cov[0, 0])


# -------------------------------------------------------------------------

t_c = 2/log(1 + sqrt(2))
use_disk = False
breadth = 1
base_duration = 50
sizes = [16, 32, 40, 50, 70]
# investigate_time_evolution()
# investigate_temperature_dependence(lattice_sizes=sizes, re_runs=4)

critical_temps = investigate_heat_capacity(lattice_sizes=sizes,
                                           take_last=300)

# temp_inf = investigate_finite_size_scaling(critical_temps, sizes)

#  print(temp_inf)
# print((t_c - temp_inf[0])/ temp_inf[1], ' Standard errors away')

Фотографий

Чтобы проиллюстрировать проблему, посмотрите здесь:

Как я уже сказал, намагничивание - это хорошо, а вот энергия - нет. Резкого перехода нет.

Теплоемкость еще хуже. Он начинает выходить на плато для

Мои мысли.

Я мог бы просто не запускать симуляцию в течение достаточно длительного времени, десятки тысяч MCS и усреднение по слишком небольшому количеству точек данных (последние 300), но тогда полосы ошибок будут больше или, по крайней мере, видны.