코드로 실험하는 일상 과학 · 1편 | 음펨바 효과와 평형에 도착하는 기준
특정 조건에서 뜨거운 물이 먼저 언다는 보고가 있다. 다만 조건과 판정 기준을 빼고 “뜨거우면 더 빨리 언다”고 일반화할 수는 없다. 이 글의 파이썬 코드는 물의 결빙을 계산하지 않는다. 입자가 어디에 있을지를 나타내는 확률분포가 평형으로 돌아가는 과정을 계산한다. 이 구분이 음펨바 효과를 이해하는 출발점이다. 물의 결빙을 둘러싼 논의와 더 넓은 완화 현상의 관계는 Kumar와 Bechhoefer의 2020년 Nature 논문에서도 확인할 수 있다.
계산에서는 더 뜨겁게 준비한 분포가 같은 평형 허용오차에 약 2,873배 빨리 들어갔다. 꽤 큰 차이다. 하지만 판정 기준을 바꾸자 덜 뜨거운 쪽이 먼저 도착하는 경우도 나왔다. 두 결과를 함께 보면 “더 빨리 식는다”는 말에 무엇이 빠져 있는지 드러난다.
먼저 결과부터 보면
환경 온도를 1로 고정하고, 초기 온도 2와 약 19.479545에서 준비한 두 분포를 비교했다. 온도와 시간은 모두 무차원 값이다. 19.48을 섭씨로 읽거나 아래 시간을 초로 읽으면 안 된다. 실제 냉동실이나 특정 입자에 맞춘 환산은 하지 않았다.
도착선은 총변동거리 TV가 0.001 이하가 되는 순간으로 정했다. TV는 현재 분포와 최종 평형 분포의 전체 차이를 재는 값이다. 0.001 이하라면 격자 위에서 어떤 영역을 고르더라도, 그곳에 입자가 있을 확률은 평형값과 최대 0.1퍼센트포인트만큼 다르다.
- 초기 온도 2: 도달 시간 1426.876
- 초기 온도 약 19.479545: 도달 시간 0.496636
- 두 도달 시간의 비: 약 2,873배
뜨거운 분포가 처음부터 유리한 위치에 있었던 것도 아니다. 초기 TV는 뜨거운 쪽이 약 0.498, 덜 뜨거운 쪽이 약 0.164였다. 전체 분포로 보면 약 세 배 더 멀리서 출발했다. 그런데도 나중에는 앞선다.
뜨거운 물은 찬물의 온도를 먼저 지나야 하지 않을까
단순한 냉각 모형에서는 맞는 직관이다. 물체의 상태를 온도 하나로 표현하고, 두 물체에 같은 환경 온도와 고정된 냉각계수를 적용해 보자. 뉴턴의 냉각 법칙을 풀면 두 온도의 차이는 다음처럼 줄어든다.
온도 차이(t) = 처음 온도 차이 × exp(−k × t), k > 0
오른쪽 값은 유한한 시간 동안 양수다. 처음 더 뜨거웠던 쪽이 계속 더 뜨겁고, 온도 곡선은 교차하지 않는다. 단지 “뜨거우면 열을 더 빨리 잃는다”는 말만으로 역전을 설명할 수 없는 이유다. 실제 물의 결빙에는 이 단순식에 넣지 않은 과정과 조건이 있으므로, 이 계산만으로 모든 컵의 결과를 판정할 수도 없다.
여기서 살펴볼 모형은 상태를 온도 하나로 압축하지 않는다. 공간 곳곳에 입자가 있을 확률을 모두 추적한다. 처음의 분포는 온도로 정하지만, 변화하는 도중의 분포에는 일반적으로 하나의 온도를 붙일 수 없다. 따라서 뜨거운 초기 분포가 덜 뜨거운 초기 분포와 똑같은 상태를 순서대로 거쳐야 할 이유도 없다.
두 우물 사이에서 오래 걸리는 일
입자가 움직이는 지형을 두 개의 골짜기와 가운데 장벽으로 만들었다. 물리에서는 이런 에너지 지형을 이중 우물 퍼텐셜이라고 부른다. 왼쪽 우물은 조금 더 깊고, 입자가 움직일 수 있는 범위도 좌우가 다르다. 입자는 양쪽 끝에서 밖으로 나가지 않는다. 역전의 구조를 보기 위해 의도적으로 구성한 모형이다.
각 우물 안에서 분포가 자리를 잡는 과정은 빠르다. 반면 왼쪽과 오른쪽에 있을 확률의 비율이 평형과 어긋나 있으면, 장벽을 넘는 이동을 통해 천천히 맞춰야 한다. 이 느린 재배분이 그래프에 긴 꼬리를 남긴다.
분포의 변화를 여러 감쇠 성분으로 나누면 가장 느리게 사라지는 성분을 찾을 수 있다. 이번 모형에서는 초기 온도를 약 19.48로 맞췄을 때 그 성분의 계수가 0을 지난다. 수치적으로 조율한 분포에는 오래 남을 성분이 거의 없어서, 더 빠른 성분들이 평형 접근을 결정한다. 이것이 2019년 PRX 논문에서 다룬 강한 음펨바 효과의 핵심이다.
개별 입자가 장벽을 순간이동으로 통과한다는 뜻은 아니다. 전체 분포에 남아 있는 느린 성분이 달라진 것이다. 더 뜨겁게 만들수록 무조건 좋아지는 것도 아니다. 초기 온도를 80으로 높이면 TV 0.001 도달 시간은 약 1043.726으로 다시 길어진다.
2,873배라는 숫자에 붙는 조건
먼저, 도착선을 바꾸면 순서가 달라질 수 있다. 분포 차이를 로그로 평가하는 KL 발산을 사용해 KL 0.01 이하를 기준으로 삼으면, 초기 온도 2의 시간은 0.07927이고 뜨거운 쪽은 0.19697이다. 이번에는 덜 뜨거운 쪽이 먼저다. TV와 KL에 같은 숫자를 넣어도 같은 엄격도의 검사가 되지는 않는다.
온도를 얼마나 정확히 맞췄는지도 중요하다. TV 0.001에서는 조율값에서 초기 온도를 ±1% 바꿔도 도달 시간이 약 0.50으로 비슷하다. 그러나 TV 0.0001까지 요구하면 조율한 경우는 약 0.679, ±1% 어긋난 경우는 약 177~184가 된다. 조금 남은 느린 성분이 더 엄격한 도착선에서는 오래 발목을 잡는다.
본문의 “약 19.48”은 읽기 위한 반올림이다. 재현 코드의 auto 옵션은 해당 격자에서 영점을 다시 찾는다. 표시된 숫자를 반올림한 채 직접 넣은 결과나, 더 촘촘한 격자의 결과가 아주 작은 허용오차까지 같으리라고 기대해서는 안 된다.
계산은 공간을 940칸으로 나눠 수행했다. 같은 초기 온도를 유지하며 1880칸으로 늘렸을 때, TV 0.001 도달 시간의 변화는 약 0.00062%였다. 별도의 시간 적분과 연속 방정식 계산도 해석을 뒷받침했다. 다만 이 검증은 구성한 모형의 수치 결과에 관한 것이다. 물이 어는 시간을 검증한 실험 데이터는 없다.
실제 연구에서는 어디까지 왔을까
Kumar와 Bechhoefer는 2020년 물속 콜로이드 입자를 이용한 제어된 계에서 빠른 평형 접근을 실험으로 보였다. 이 글의 코드는 그 연구와 강한 음펨바 이론에서 출발했지만, 논문의 실험 데이터를 재현하거나 같은 매개변수를 복제한 것은 아니다. Nature 논문과 공개 원고를 함께 읽을 수 있다.
2026년 3월 25일에는 고전계와 양자계의 음펨바 효과를 자원 이론으로 통합해 설명하는 연구도 PRX에 실렸다. 연구의 범위가 일상적인 물의 결빙을 넘어 확장되고 있다는 사례다. 이를 가정용 냉동실에서 뜨거운 물이 언제 먼저 어는지에 대한 새 해답으로 읽을 필요는 없다. Resource-Theoretical Unification of Mpemba Effects: Classical and Quantum.
파이썬으로 도착선을 바꿔 보기
아래의 run_experiment.py 전체 코드와 compute.py 전체 코드를 펼쳐 복사한 뒤, 각각의 파일명으로 같은 폴더에 UTF-8 텍스트로 저장한다. Python 3.10 이상과 NumPy, SciPy가 필요하다. Matplotlib이 있으면 결과 그래프도 저장된다.
터미널에서 두 파일이 있는 폴더로 이동한 뒤 기본 비교를 실행한다.
python run_experiment.py --warm 2 --hot auto --metric TV --epsilon 0.001 --grid 940 --output runs/basicruns/basic 안의 summary.txt에서 두 도달 시간을 읽을 수 있다. result.json에는 입력과 점검 결과가, arrivals.csv와 timeseries.csv에는 비교 수치가 저장된다. 다음 명령은 앞서 본 KL 기준의 순서 역전을 확인한다.
python run_experiment.py --warm 2 --hot auto --metric KL --epsilon 0.01 --grid 940 --output runs/kl_001결과 폴더가 이미 있으면 덮어쓰지 않고 멈춘다. 다시 실행할 때는 --output 뒤에 새 폴더명을 쓰면 된다. 초기 온도를 고정하려면 --hot auto를 --hot 19.2847496403717로 바꾸고, TV 기준을 0.0001로 낮춰 온도 오차의 영향을 살펴볼 수 있다.
이번 계산에서 확인한 것은 특정한 초기 분포가 느린 완화 성분을 줄여, 더 멀리서도 같은 평형 기준에 먼저 들어갈 수 있다는 사실이다. 약 2,873배라는 결과 옆에 모형, 거리, 허용오차를 함께 적어야 하는 이유도 여기에 있다. 코드를 다시 돌릴 때는 속도비 하나보다, 조건을 바꾸었을 때 그 숫자가 어떻게 달라지는지부터 비교해 보면 좋겠다.
직접 실행할 전체 코드
각 항목을 펼치면 생략 없는 원본 코드가 나온다. 코드 블록 안의 내용만 복사해 아래 파일명으로 저장한다. 두 파일을 같은 폴더에 둔 뒤 위의 실행 명령을 사용한다. 이 코드는 네트워크 접속이나 API 키 없이 로컬에서 계산한다.
run_experiment.py
run_experiment.py 전체 코드 펼치기
#!/usr/bin/env python3
"""A parameterized, deterministic experiment using the frozen compute.py solver.
Python >=3.10, NumPy, SciPy; Matplotlib is optional. All quantities are
nondimensional. This is a constructed numerical model, not a water experiment.
"""
from __future__ import annotations
import argparse
import csv
import hashlib
import importlib.util
import json
import math
import os
from pathlib import Path
import platform
import sys
# Set defaults before loading numerical libraries; an explicit user choice wins.
for _name in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ.setdefault(_name, "1")
GRIDS = (235, 470, 940, 1880)
TEMPERATURE_RANGE = (0.1, 200.0)
MIN_EPSILON = 1e-8
MAX_TIME = 1e7
def _parser() -> argparse.ArgumentParser:
p = argparse.ArgumentParser(
description="시작 온도와 도달 기준을 바꾸는 강한 음펨바 수치실험. 모든 값은 무차원입니다.",
epilog="예: python run_experiment.py --warm 2 --hot auto --metric KL --epsilon .01 --output runs/kl",
)
p.add_argument("--warm", type=float, default=2.0, help="warm 시작 온도: 0.1~200 (기본 2)")
p.add_argument("--hot", default="auto", help="hot 시작 온도: 0.1~200 또는 auto (해당 격자의 비자명근)")
p.add_argument("--metric", type=str.upper, choices=("TV", "KL"), default="TV")
p.add_argument("--epsilon", type=float, default=1e-3, help="도달 허용오차. TV: 1e-8~1; KL: 1e-8~10")
p.add_argument("--grid", type=int, choices=GRIDS, default=940, help="격자 셀 수 (기본 940; 모든 모드 사용)")
p.add_argument("--output", type=Path, default=Path("runs/default"), help="새 결과 폴더 (기본 runs/default)")
p.add_argument("--no-plot", action="store_true", help="선택적 PNG 그림 생성을 건너뜁니다")
return p
def _validate(p: argparse.ArgumentParser, args: argparse.Namespace) -> None:
low, high = TEMPERATURE_RANGE
if not math.isfinite(args.warm) or not low <= args.warm <= high:
p.error(f"--warm은 양의 유한한 온도여야 하며 지원 범위는 {low}~{high}입니다.")
if str(args.hot).lower() == "auto":
args.hot = "auto"
else:
try:
args.hot = float(args.hot)
except (ValueError, TypeError):
p.error("--hot은 auto 또는 숫자여야 합니다.")
if not math.isfinite(args.hot) or not low <= args.hot <= high:
p.error(f"--hot은 양의 유한한 온도여야 하며 지원 범위는 {low}~{high}입니다.")
maximum = 1.0 if args.metric == "TV" else 10.0
if not math.isfinite(args.epsilon) or not MIN_EPSILON <= args.epsilon <= maximum:
p.error(f"{args.metric}의 --epsilon 지원 범위는 {MIN_EPSILON:g}~{maximum:g}입니다.")
def _load_solver():
here = Path(__file__).resolve().parent
# Flat release layout first; parent fallback is only for source development.
candidates = (here / "compute.py", here.parent / "compute.py")
path = next((f for f in candidates if f.is_file()), None)
if path is None:
raise RuntimeError("compute.py를 run_experiment.py와 같은 폴더에 놓아 주세요.")
spec = importlib.util.spec_from_file_location("mpemba_frozen_compute", path)
module = importlib.util.module_from_spec(spec)
spec.loader.exec_module(module)
return module, path
def _arrival_time(model, temperature: float, epsilon: float, metric: str) -> float:
"""Same continuous-time criterion and Brent tolerances as frozen solver.
The explicit bracket cap makes even an unexpected numerical failure bounded.
No arrival time is inferred from the exported sampled plotting grid.
"""
from scipy.optimize import brentq
index = 0 if metric == "TV" else 1
def f(t):
value = float(model.metrics(temperature, t)[index]) - epsilon
if not math.isfinite(value):
raise RuntimeError("거리 계산이 유한하지 않습니다. 입력/수치 환경을 확인하세요.")
return value
if f(0.0) <= 0.0:
return 0.0
hi = 1.0
for _ in range(25):
if f(hi) <= 0.0:
return float(brentq(f, 0.0, hi, xtol=1e-10, rtol=1e-12))
hi = min(2.0 * hi, MAX_TIME)
raise RuntimeError(f"t<={MAX_TIME:g}에서 도달 구간을 찾지 못했습니다. 도달했다고 간주하지 않습니다.")
def _json_write(path: Path, data) -> None:
path.write_text(json.dumps(data, ensure_ascii=False, indent=2, allow_nan=False) + "\n", encoding="utf-8")
def _csv_write(path: Path, header, rows) -> None:
with path.open("w", encoding="utf-8", newline="") as f:
writer = csv.writer(f)
writer.writerow(header)
writer.writerows(rows)
def _plot(path, samples, states, metric, epsilon):
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return {"status": "unavailable", "reason": "Matplotlib is not installed; numerical results are complete."}
fig, ax = plt.subplots(figsize=(9.6, 5.6), constrained_layout=True)
colors = {"warm": "#087e8b", "hot": "#d95d39"}
for state in states:
name = state["label"]
ax.plot(samples[name]["time"], samples[name][metric], color=colors[name], linewidth=2,
label=f'{name}: T={state["temperature"]:.8g}; arrival={state["arrival_time"]:.6g}')
ax.scatter([state["arrival_time"]], [state["distance_at_arrival"]],
color=colors[name], s=38, zorder=4)
ax.axhline(epsilon, linestyle="--", color="#555555", label=f"criterion: {metric} <= {epsilon:g}")
# Symmetric-log axes have a linear neighborhood of zero, preserving t=0
# and exact-equilibrium distances without replacing zeros with fake values.
ax.set_xscale("symlog", linthresh=0.01)
ax.set_yscale("symlog", linthresh=1e-8)
ax.set_xlabel("Time (dimensionless; symlog axis, linear near 0)")
ax.set_ylabel(f"{metric} distance to the same equilibrium (symlog axis)")
ax.set_title("Change the initial state. Recompute the arrival time.", loc="left", fontweight="bold", pad=30)
ax.text(0.0, 1.02, "Constructed numerical experiment | bath T=1 | no seconds or kelvin", transform=ax.transAxes, fontsize=9)
ax.set_xlim(0.0, max(float(samples[s["label"]]["time"][-1]) for s in states))
ax.set_ylim(0.0, 1.2 * max(epsilon, *(float(samples[s["label"]][metric].max()) for s in states)))
ax.grid(alpha=0.18)
ax.legend(fontsize=9)
fig.savefig(path, dpi=160)
plt.close(fig)
return {"status": "created", "file": path.name, "axes": "symlog; exact zeros are retained"}
def _summary(result) -> str:
states = result["states"]
warm, hot = states
params = result["parameters"]
comparison = result["comparison"]
lines = [
"음펨바 수치실험 결과",
f'기준: {params["metric"]} <= {params["epsilon"]:g} | 격자 {params["grid"]}개 | 열탕 T=1',
f'자동 탐색 비자명근: T*={result["nontrivial_root_temperature"]:.12g}',
]
for s in states:
note = " (시작부터 기준 충족)" if s["initially_within_threshold"] else ""
lines.append(f'{s["label"]}: T={s["temperature"]:.12g}, 처음 거리={s["initial_distance"]:.9g}, 도달시간={s["arrival_time"]:.12g}{note}')
if comparison["winner"] == "both_already_within_threshold":
lines.append("두 입력 모두 시작부터 기준을 충족합니다. 더 빨리 도달한 입력은 없습니다.")
elif comparison["winner"] == "tie":
lines.append("수치 시간 허용오차 안에서 두 입력의 도달시간이 같습니다.")
else:
lines.append(f'{comparison["winner"]} 입력이 먼저 기준에 도달합니다.')
ratio = comparison["warm_time_over_hot_time"]
if ratio is None:
lines.append("warm 시간 / hot 시간: 정의하지 않음 (분모인 hot 시간이 0)")
else:
lines.append(f"warm 시간 / hot 시간 = {ratio:.12g} (시간의 비율이며 보편적 가속 배수가 아님)")
for warning in result["notes"]:
lines.append("참고: " + warning)
lines.append("온도·시간은 모두 무차원입니다. 물의 동결 실험이나 새로운 발견을 뜻하지 않습니다.")
lines.append("결과: result.json, arrivals.csv, timeseries.csv, summary.txt" + (", distances.png" if result["plot"]["status"] == "created" else ""))
return "\n".join(lines) + "\n"
def run_experiment(args: argparse.Namespace) -> dict:
import numpy as np
import scipy
solver, solver_path = _load_solver()
model = solver.Model(args.grid, full=True)
if not np.all(np.isfinite(model.rate)) or not np.all(model.rate > 0):
raise RuntimeError("유효하지 않은 감쇠율입니다. 결과를 저장하지 않습니다.")
root = float(model.root())
temperatures = (float(args.warm), root if args.hot == "auto" else float(args.hot))
index = 0 if args.metric == "TV" else 1
states = []
for name, temperature in zip(("warm", "hot"), temperatures):
initial_tv, initial_kl = (float(x) for x in model.metrics(temperature, 0.0))
initial = (initial_tv, initial_kl)[index]
arrival = _arrival_time(model, temperature, args.epsilon, args.metric)
states.append({"label": name, "temperature": temperature,
"initial_TV": initial_tv, "initial_KL": initial_kl,
"initial_distance": initial, "initially_within_threshold": initial <= args.epsilon,
"exact_bath_temperature": temperature == solver.TB,
"slowest_mode_coefficient_a2": float(model.slow(temperature)),
"arrival_time": arrival,
"distance_at_arrival": float(model.metrics(temperature, arrival)[index])})
tw, th = (state["arrival_time"] for state in states)
time_tolerance = max(2e-10, 1e-11 * max(tw, th))
winner = ("both_already_within_threshold" if tw == th == 0 else
"tie" if abs(tw - th) <= time_tolerance else "warm" if tw < th else "hot")
comparison = {"winner": winner, "time_comparison_tolerance": time_tolerance,
"warm_time_over_hot_time": tw / th if th > 0 else None,
"ratio_undefined_reason": "hot_arrival_time_is_zero" if th == 0 else None,
"hot_temperature_is_higher": temperatures[1] > temperatures[0],
"hot_initial_distance_is_larger": states[1]["initial_distance"] > states[0]["initial_distance"]}
notes = []
if temperatures[1] <= temperatures[0]:
notes.append("hot/warm은 입력 이름입니다. 이 설정에서 hot 온도는 warm보다 높지 않습니다.")
if min(temperatures) < solver.TB:
notes.append("열탕보다 낮은 시작 온도를 포함하므로 두 상태가 모두 냉각되는 비교는 아닙니다.")
if winner == "warm":
notes.append("이 조건에서는 warm이 먼저 도달합니다. 지표와 허용오차에 따라 순위가 달라집니다.")
if args.hot == "auto":
notes.append("auto는 선택한 격자에서 a2=0이 되는 온도를 다시 찾습니다. 실험 장치의 정확한 온도를 뜻하지 않습니다.")
notes.append("아주 엄격한 기준에서는 온도 오차와 격자 오차가 느린 꼬리를 되살릴 수 있습니다.")
if args.grid != 940:
notes.append("기본 검증 격자는 940개입니다. 변경한 격자의 결과를 연속계의 정확한 값으로 간주하지 마세요.")
horizon = max(1.5, 1.15 * max(tw, th))
times = np.unique(np.r_[0.0, np.geomspace(1e-6, horizon, 321), np.linspace(0, min(1.5, horizon), 181), tw, th])
samples, rows = {}, []
minimum_probability, mass_error = math.inf, 0.0
for state in states:
probability = model.distribution(state["temperature"], times)
tv = 0.5 * np.abs(probability - model.pi).sum(axis=1)
safe = np.maximum(probability, np.finfo(float).tiny)
kl = (safe * np.log(safe / model.pi)).sum(axis=1)
energy = probability @ model.U
samples[state["label"]] = {"time": times, "TV": tv, "KL": kl}
rows.extend((state["label"], state["temperature"], float(t), float(a), float(b), float(e))
for t, a, b, e in zip(times, tv, kl, energy))
minimum_probability = min(minimum_probability, float(probability.min()))
mass_error = max(mass_error, float(np.abs(probability.sum(axis=1) - 1).max()))
if mass_error > 1e-10 or minimum_probability < -1e-10:
raise RuntimeError("확률 보존/비음성 수치 점검 실패. 이 결과를 신뢰할 수 없어 저장하지 않습니다.")
output = args.output.expanduser().resolve()
# Each experiment is a new folder: never overwrite a previous result.
output.mkdir(parents=True, exist_ok=False)
result = {
"schema_version": 1,
"classification": "Constructed deterministic numerical experiment, not measured data or a new discovery",
"parameters": {"warm": temperatures[0], "hot_requested": args.hot, "hot": temperatures[1],
"metric": args.metric, "epsilon": args.epsilon, "grid": args.grid, "bath_temperature": solver.TB},
"units": "All temperatures, positions, potentials and times are nondimensional; no kelvin/seconds conversion",
"arrival_definition": "First continuous time D(p(t), pi)<=epsilon, bracketed root solve; exactly zero if already inside at t=0",
"nontrivial_root_temperature": root,
"nontrivial_root_a2_residual": float(model.slow(root)),
"spectrum": {"lambda2": float(model.rate[0]), "lambda3": float(model.rate[1]), "modes_used": len(model.rate)},
"states": states, "comparison": comparison, "notes": notes,
"sampled_checks": {"sample_count_per_state": len(times), "last_sample_time": horizon,
"max_mass_error": mass_error, "minimum_probability": minimum_probability,
"note": "Checks cover exported samples; the solver's separate audit supplies broader validation."},
"reproducibility": {"python": platform.python_version(), "numpy": np.__version__, "scipy": scipy.__version__,
"solver_sha256": hashlib.sha256(solver_path.read_bytes()).hexdigest(),
"driver_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
"randomness": "none", "all_nonstationary_modes": True},
"plot": {"status": "skipped"},
}
_csv_write(output / "arrivals.csv", ["label", "temperature", "metric", "epsilon", "initial_distance", "initially_within_threshold", "arrival_time", "a2"],
((s["label"], s["temperature"], args.metric, args.epsilon, s["initial_distance"], s["initially_within_threshold"], s["arrival_time"], s["slowest_mode_coefficient_a2"]) for s in states))
_csv_write(output / "timeseries.csv", ["label", "initial_temperature", "time", "TV", "KL", "mean_potential_energy"], rows)
if not args.no_plot:
result["plot"] = _plot(output / "distances.png", samples, states, args.metric, args.epsilon)
_json_write(output / "result.json", result)
(output / "summary.txt").write_text(_summary(result), encoding="utf-8")
return result
def main(argv=None) -> int:
parser = _parser()
args = parser.parse_args(argv)
_validate(parser, args)
output = args.output.expanduser()
if output.exists():
parser.error(f"결과 폴더가 이미 있습니다: {output}. 덮어쓰지 않습니다. 새 --output 경로를 지정하세요.")
try:
result = run_experiment(args)
except ImportError as exc:
print(f"필수 과학 계산 패키지가 없습니다: {exc}. Python 3.10+, NumPy, SciPy 환경이 필요합니다.", file=sys.stderr)
return 2
except (RuntimeError, ValueError, OSError) as exc:
print(f"계산을 완료하지 못했습니다: {exc}", file=sys.stderr)
return 2
print(_summary(result), end="")
print(f"저장 위치: {args.output.expanduser().resolve()}")
if result["plot"]["status"] == "unavailable":
print("Matplotlib이 없어 그림만 생략했습니다. JSON/CSV 수치 결과는 저장되었습니다.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
compute.py
compute.py 전체 코드 펼치기
#!/usr/bin/env python3
"""Deterministic strong-Mpemba computation. Python >=3.10; numpy, scipy only.
Run: OPENBLAS_NUM_THREADS=1 python compute.py
All x, U, T, t are nondimensional. Nothing here is experimental observation.
"""
from pathlib import Path
import json,csv,platform,time
import numpy as np
import scipy
from scipy.linalg import eigh_tridiagonal
from scipy.optimize import brentq
from scipy.sparse import diags
from scipy.sparse.linalg import expm_multiply
OUT=Path(__file__).resolve().parent
BARRIER=8.; TILT=.2; OUTER=4.; XLO=-3.; XHI=1.7; TB=1.; MOBILITY=1.
N_MAIN=940
THRESHOLDS=[1e-2,1e-3,1e-4,1e-5,1e-6]
def potential(x,tilt=TILT):
x=np.asarray(x)
return np.where(abs(x)<=1,BARRIER*(x*x-1)**2,OUTER*(abs(x)-1)**2)+tilt*x
def boltzmann(U,T):
p=np.exp(-(U-U.min())/T);return p/p.sum()
class Model:
def __init__(self,n=N_MAIN,xlo=XLO,xhi=XHI,tilt=TILT,full=True):
self.n=n; self.h=(xhi-xlo)/n; self.x=xlo+(np.arange(n)+.5)*self.h
self.U=potential(self.x,tilt); self.pi=boltzmann(self.U,TB); self.sq=np.sqrt(self.pi)
du=np.diff(self.U)/TB
self.rp=MOBILITY*TB/self.h**2*np.exp(-du/2)
self.rm=MOBILITY*TB/self.h**2*np.exp(du/2)
self.diag=-np.r_[self.rp,0]-np.r_[0,self.rm]
self.off=np.sqrt(self.rp*self.rm)
# stebz + inverse iteration; tight absolute eigenvalue tolerance.
kwargs={} if full else dict(select='i',select_range=(n-12,n-1))
val,V=eigh_tridiagonal(self.diag,self.off,lapack_driver='stebz',tol=1e-13,**kwargs)
val=val[::-1];V=V[:,::-1]
self.raw_zero=float(val[0]); self.rate=-val[1:];self.V=V[:,1:]
# Fix modal sign to be positive in the right-hand basin.
for k in range(self.V.shape[1]):
if np.dot(self.V[self.x>0,k],self.sq[self.x>0])<0:self.V[:,k]*=-1
self.R=self.sq[:,None]*self.V
# Exact mass conservation despite eigensolver nullspace roundoff.
self.R-=self.pi[:,None]*self.R.sum(axis=0)
self.L=self.V/self.sq[:,None]
self.L-=np.dot(self.pi,self.L)[None,:]
def p0(self,T):return boltzmann(self.U,T)
def coefficients(self,T):return (self.p0(T)-self.pi)@self.L
def slow(self,T):return float(self.coefficients(T)[0])
def root(self):return brentq(self.slow,5.,60.,xtol=5e-12,rtol=1e-14)
def distribution(self,T,t):
ts=np.atleast_1d(t);c=self.coefficients(T)
p=self.pi[None,:]+(np.exp(-ts[:,None]*self.rate[None,:])*c[None,:])@self.R.T
p[ts==0]=self.p0(T)
return p[0] if np.ndim(t)==0 else p
def metrics(self,T,t):
p=self.distribution(T,t);q=np.maximum(p,np.finfo(float).tiny)
return (.5*np.abs(p-self.pi).sum(axis=-1), (q*np.log(q/self.pi)).sum(axis=-1))
def threshold_time(self,T,epsilon,metric='TV'):
i=0 if metric=='TV' else 1
f=lambda t:float(self.metrics(T,t)[i]-epsilon)
if f(0)<=0:return 0.
hi=1.
while f(hi)>0:hi*=2
return float(brentq(f,0.,hi,xtol=1e-10,rtol=1e-12))
def write_csv(path,header,rows):
with open(path,'w',newline='') as f:
w=csv.writer(f);w.writerow(header);w.writerows(rows)
def main():
tic=time.time();m=Model();root=m.root()
# Tilt moves the central maximum slightly to the right of x=0.
barrier_x=brentq(lambda x:32*x*(x*x-1)+TILT,-.5,.5)
eq_E=float(m.pi@m.U)
temps=np.array([2.,root,root*.99,root*1.01,80.])
names=np.array(['warm','strong_hot','hot_minus_1pct','hot_plus_1pct','hot_80'])
t=np.unique(np.r_[0.,np.geomspace(1e-5,10000.,401),np.linspace(0,1.5,181)])
p=np.stack([m.distribution(T,t) for T in temps])
tv=.5*np.abs(p-m.pi).sum(axis=-1)
safe=np.maximum(p,np.finfo(float).tiny)
kl=(safe*np.log(safe/m.pi)).sum(axis=-1)
energy=p@m.U
left=p[:,:,m.x<barrier_x].sum(axis=-1)
cs=np.stack([m.coefficients(T) for T in temps])
cutoff=float(m.pi[m.x<barrier_x].sum())
states=[]
for i,(name,T) in enumerate(zip(names,temps)):
tt={format(e,'.0e'):m.threshold_time(T,e) for e in THRESHOLDS}
kt={format(e,'.0e'):m.threshold_time(T,e,'KL') for e in THRESHOLDS}
states.append(dict(name=str(name),temperature=float(T),initial_TV=float(tv[i,0]),initial_KL=float(kl[i,0]),initial_mean_energy=float(energy[i,0]),initial_left_population=float(left[i,0]),a2=float(cs[i,0]),a3=float(cs[i,1]),TV_threshold_times=tt,KL_threshold_times=kt))
convergence=[]
for n in [470,940,1880]:
mm=m if n==N_MAIN else Model(n,full=False)
rr=mm.root();warm=mm.threshold_time(2,1e-3);hot=mm.threshold_time(rr,1e-3)
fixeda=mm.slow(root);fixedt=mm.threshold_time(root,1e-3)
convergence.append(dict(n=n,dx=mm.h,root_temperature=rr,lambda2=float(mm.rate[0]),lambda3=float(mm.rate[1]),spectral_gap_ratio=float(mm.rate[1]/mm.rate[0]),warm_TV_1e3_time=warm,retuned_hot_TV_1e3_time=hot,fixed_main_root_temperature=root,fixed_main_root_a2=fixeda,fixed_main_root_TV_1e3_time=fixedt,fixed_main_root_TV_1e6_time=mm.threshold_time(root,1e-6),fixed_main_root_TV_1e8_time=mm.threshold_time(root,1e-8)))
# Independent matrix-exponential propagation on a smaller grid to control runtime.
small=Model(235);smallroot=small.root()
Q=diags([small.rp,small.diag,small.rm],[-1,0,1],format='csc')
expm_checks=[]
for T in [2.,smallroot]:
# One Krylov/Taylor expm_multiply call evaluates uniform time points.
ts=np.linspace(0,1,11);reference=expm_multiply(Q,small.p0(T),start=0,stop=1,num=11,endpoint=True)
spectral=small.distribution(T,ts)
expm_checks.append(dict(n=235,temperature=T,times=ts.tolist(),max_abs_probability_error=float(np.abs(reference-spectral).max()),max_TV_between_solvers=float(.5*np.abs(reference-spectral).sum(axis=1).max()),max_mass_error=float(np.abs(reference.sum(axis=1)-1).max()),minimum_probability=float(reference.min())))
stationary_res=np.r_[0,m.rp*m.pi[:-1]]+np.r_[m.rm*m.pi[1:],0]+m.diag*m.pi
balance=np.abs(m.rp*m.pi[:-1]-m.rm*m.pi[1:])
dbscale=np.maximum(m.rp*m.pi[:-1],m.rm*m.pi[1:])
check=dict(stationary_residual_infinity=float(np.abs(stationary_res).max()),detailed_balance_max_absolute=float(balance.max()),detailed_balance_max_relative=float((balance/dbscale).max()),mass_error_max=float(np.abs(p.sum(axis=-1)-1).max()),minimum_probability=float(p.min()),minimum_transition_rate=float(min(m.rp.min(),m.rm.min())),raw_stationary_eigenvalue=m.raw_zero,slow_left_equilibrium_mean=float(m.pi@m.L[:,0]),slow_left_variance=float(m.pi@(m.L[:,0]**2)),max_TV_upward_step=float(np.diff(tv,axis=1).max()),max_KL_upward_step=float(np.diff(kl,axis=1).max()),initial_reconstruction_max_error=float(np.max(np.abs(m.pi+m.R@m.coefficients(2)-m.p0(2)))),generator_boundary_outward_rates=0.)
Ts=np.unique(np.r_[np.geomspace(1,200,301),root,2.]);a2=np.array([m.slow(T) for T in Ts])
sweepTV=np.array([m.threshold_time(T,1e-3) for T in Ts])
# +/- 1% detuning and local derivative quantify why exact cancellation is tuned.
derivative=(m.slow(root*1.0001)-m.slow(root*.9999))/(.0002*root)
detuning=[]
for delta in [-.05,-.01,-.001,0.,.001,.01,.05]:
T=root*(1+delta)
detuning.append(dict(relative_temperature_change=delta,temperature=T,a2=m.slow(T),TV_1e3_time=m.threshold_time(T,1e-3),TV_1e4_time=m.threshold_time(T,1e-4),TV_1e6_time=m.threshold_time(T,1e-6)))
crossing_TV=brentq(lambda tt:float(m.metrics(2,tt)[0]-m.metrics(root,tt)[0]),.01,1.)
crossing_KL=brentq(lambda tt:float(m.metrics(2,tt)[1]-m.metrics(root,tt)[1]),.01,1.)
metric_caveat=dict(TV_distance_curve_crossing_time=crossing_TV,KL_distance_curve_crossing_time=crossing_KL,description='Hot starts farther and overtakes later. Arrival ranking depends on metric and tolerance: for example KL<=0.01 is reached first by warm, but TV<=0.001 is reached far sooner by strong_hot. No claim that hot has smaller distance at all times or at every possible cutoff.',energy_note='Mean potential energy is only one observable. It does not uniquely determine the full probability distribution and is not assigned an evolving temperature.')
# Symmetric geometry + zero tilt: parity eliminates odd slow mode for EVERY Boltzmann state.
sym=Model(600,xlo=-3,xhi=3,tilt=0,full=False)
symcontrol=dict(description='Reflecting symmetric [-3,3], same even core/tails, zero tilt. Every Boltzmann initial density is even, so odd slow mode has zero overlap; this is NOT an isolated strong-Mpemba root.',temperatures=[2.,root,80.],a2=[sym.slow(T) for T in [2.,root,80.]],lambda2=float(sym.rate[0]))
summary=dict(title='A constructed continuous strong-Mpemba numerical experiment',classification='Reproducible numerical demonstration; not a physical experiment, not a new discovery, not a result about freezing water.',model=dict(potential_core='U(x)=8*(x^2-1)^2+0.2*x, abs(x)<=1',potential_tails='U(x)=4*(abs(x)-1)^2+0.2*x, abs(x)>1',domain=[XLO,XHI],boundary='reflecting, zero probability flux',regularity='U and U_prime continuous; U_double_prime piecewise continuous with jumps at x=±1',bath_temperature=TB,mobility=MOBILITY,diffusion_coefficient=MOBILITY*TB,units='nondimensional: k_B=1, mobility=1; no mapping to seconds, kelvin or physical particles',fokker_planck='partial_t rho = partial_x(U_prime*rho + T_b*partial_x rho)',equilibrium='rho_eq(x) proportional to exp(-U(x)/T_b)',initial='rho(x,0;T_i) proportional to exp(-U(x)/T_i)',barrier_maximum_x=barrier_x,barrier_potential=float(potential(barrier_x)),equilibrium_mean_energy=eq_E,equilibrium_left_population=cutoff,left_population_definition='cell centers x < exact central barrier maximum; cell-mass quadrature'),discretization=dict(n=N_MAIN,dx=m.h,cell_centers=True,probability_convention='p_i is a bin probability mass, density=p_i/dx',q_i_plus_1_from_i='D/dx^2 * exp(-(U_(i+1)-U_i)/(2*T_b))',q_i_from_i_plus_1='D/dx^2 * exp(+(U_(i+1)-U_i)/(2*T_b))',generator_convention='column generator Q: dp/dt=Q p, columns sum to zero',accuracy='second-order consistent detailed-balance finite-volume nearest-neighbor discretization on uniform cells',symmetrization='S=diag(pi)^(-1/2) Q diag(pi)^(1/2); S symmetric tridiagonal',time_solution='p(t)=pi + sum_(k>=2) a_k exp(-lambda_k*t) r_k, with tiny numerical mean projected out of r_k',left_mode_normalization='sum_i pi_i*l_k_i=0 and sum_i pi_i*l_k_i^2=1; a_k=sum_i (p0_i-pi_i)*l_k_i'),temperature_root=dict(value=root,bracket=[5.,60.],a2_residual=float(m.slow(root)),d_a2_d_temperature=float(derivative),interpretation='Nontrivial finite-grid zero of the slowest nonstationary coefficient; trivial equilibrium zero at T_i=T_b excluded. Continuum behavior is supported by mesh refinement and independent validation.',temperature_tuning_note='Exact cancellation is idealized; finite precision, mesh error and temperature detuning reintroduce a small late slow tail.'),spectrum=dict(first_12_decay_rates=m.rate[:12].tolist(),lambda2=float(m.rate[0]),lambda3=float(m.rate[1]),tau2=float(1/m.rate[0]),tau3=float(1/m.rate[1]),rate_ratio=float(m.rate[1]/m.rate[0])),states=states,TV_speedup_warm_over_strong_hot={format(e,'.0e'):states[0]['TV_threshold_times'][format(e,'.0e')]/states[1]['TV_threshold_times'][format(e,'.0e')] for e in THRESHOLDS},threshold_definition='First time D[p(t),pi] <= epsilon. TV and KL contract under this time-homogeneous Markov semigroup, so first entry implies remaining below; no arbitrary finite observation-window cutoff. Exact equilibrium is asymptotic.',checks=check,independent_expm_checks=expm_checks,grid_convergence=convergence,symmetric_control=symcontrol,temperature_detuning=detuning,metric_and_energy_caveats=metric_caveat,reproducibility=dict(python=platform.python_version(),numpy=np.__version__,scipy=scipy.__version__,randomness='None; deterministic discretization and propagation',command='OPENBLAS_NUM_THREADS=1 python compute.py',runtime_seconds=time.time()-tic),sources=[dict(title='Kumar & Bechhoefer, Exponentially faster cooling in a colloidal system (2020)',url='https://doi.org/10.1038/s41586-020-2560-x',role='Experimental precedent and physics context; our landscape is explicitly constructed, not their exact numerical parameters.'),dict(title='Kumar & Bechhoefer, arXiv:2008.02373',url='https://arxiv.org/abs/2008.02373',role='Accessible primary manuscript; no claim of exact replication.'),dict(title='Klich, Raz, Hirschberg & Vucelja, Mpemba Index and Anomalous Relaxation (2019)',url='https://doi.org/10.1103/PhysRevX.9.021060',role='Strong Mpemba effect: zero projection on the slowest relaxation mode.'),dict(title='Lu & Raz, Nonequilibrium thermodynamics of the Markovian Mpemba effect and its inverse (2017)',url='https://doi.org/10.1073/pnas.1701264114',role='Spectral relaxation framework for the Markovian Mpemba effect.')])
(OUT/'results.json').write_text(json.dumps(summary,indent=2,ensure_ascii=False))
np.savez_compressed(OUT/'trajectories.npz',x=m.x,U=m.U,pi=m.pi,dx=m.h,t=t,temperatures=temps,state_names=names,p=p,TV=tv,KL=kl,energy=energy,left_population=left,slow_mode_profile=m.R[:,0],slow_left_mode=m.L[:,0],modal_coefficients=cs,decay_rates=m.rate,barrier_x=barrier_x,sweep_temperature=Ts,sweep_a2=a2,sweep_TV_1e3_time=sweepTV)
write_csv(OUT/'landscape_equilibrium.csv',['x','U','pi_bin_mass','rho_equilibrium','slow_right_mode_bin','slow_left_mode'],zip(m.x,m.U,m.pi,m.pi/m.h,m.R[:,0],m.L[:,0]))
write_csv(OUT/'time_series.csv',['state','initial_temperature','time','TV','KL','mean_energy','left_population'],((names[i],temps[i],t[j],tv[i,j],kl[i,j],energy[i,j],left[i,j]) for i in range(len(temps)) for j in range(len(t))))
write_csv(OUT/'temperature_sweep.csv',['initial_temperature','a2','TV_1e3_arrival_time'],zip(Ts,a2,sweepTV))
write_csv(OUT/'temperature_detuning.csv',list(detuning[0]),([r[k] for k in detuning[0]] for r in detuning))
write_csv(OUT/'convergence.csv',list(convergence[0]),([r[k] for k in convergence[0]] for r in convergence))
write_csv(OUT/'threshold_times.csv',['state','initial_temperature','metric','epsilon','first_entry_and_remain_time'],((s['name'],s['temperature'],metric,e,s[metric+'_threshold_times'][format(e,'.0e')]) for s in states for metric in ['TV','KL'] for e in THRESHOLDS))
print(json.dumps(dict(root=root,spectrum=summary['spectrum'],states=states[:2],speedup=summary['TV_speedup_warm_over_strong_hot'],checks=check,convergence=convergence,runtime=time.time()-tic),indent=2))
if __name__=='__main__':main()