用可核验区间寻找物理根

追踪二分迭代，对照参考求解器，并从冻结的 SI 常数推导维恩位移。

数学与统计 · 第 02 课🔋 离线核心分析⏱ 50 分钟参考常数与模型📦 numpy · pandas · scipy · matplotlib · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.optimize import brentq
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ma02-bracketed-roots'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

读取内嵌快照

In [2] · Python

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'constants': '2533802a1a546c612d06ba2e931e4f32a418be854c75a5e9c327ff57d5b56126'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['constants']

方法函数

In [3] · Python

def bisection(function,lower,upper,tolerance=1e-10,max_iterations=100):
    """Bracket a continuous scalar root; retain each actual interval update."""
    if not np.isfinite([lower,upper,tolerance]).all() or lower>=upper or tolerance<=0:raise ValueError('Invalid bracket or tolerance')
    f_lower=float(function(lower));f_upper=float(function(upper))
    if not np.isfinite([f_lower,f_upper]).all():raise ValueError('Non-finite endpoint values')
    if f_lower==0:return lower,[]
    if f_upper==0:return upper,[]
    if np.sign(f_lower)==np.sign(f_upper):raise ValueError('The bracket does not change sign')
    history=[]
    for iteration in range(max_iterations):
        midpoint=(lower+upper)/2;value=float(function(midpoint))
        if not np.isfinite(value):raise ValueError('Non-finite function value')
        history.append([iteration,lower,upper,midpoint,value,upper-lower])
        if value==0 or (upper-lower)/2<=tolerance:return midpoint,history
        if np.sign(value)==np.sign(f_lower):lower=midpoint;f_lower=value
        else:upper=midpoint;f_upper=value
    raise RuntimeError('Requested root tolerance not reached')

print("Method ready: bisection")

Method ready: bisection

数据与模型边界

解析函数与 Monte Carlo 点属于数学教学模型，不是实验测量。使用 Iris 观测时，可交换性和采样设计仍是统计假设；方便样本不自动代表总体。数值离散误差与采样不确定性应分别解释。

分步分析

1. 排除平凡根

普朗克波长形式给出 5(1−exp(−x))−x=0。x=0 是数学根，但不是本课所求的正谱峰。区间 [1,10] 存在变号，可隔离所需正根。

In [4] · Python

constants=snapshot("constants")["constants"]
def equation(value):return 5*(-np.expm1(-value))-value
root,history=bisection(equation,1.0,10.0,tolerance=1e-10)
reference=brentq(equation,1.0,10.0,xtol=1e-13)
trace=pd.DataFrame(history,columns=["iteration","lower","upper","midpoint","f_midpoint","width"])
print("Positive root:",root,"; independent brentq:",reference)
display(trace.head(8))
grid=np.linspace(0,10,501)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(grid,equation(grid));ax.axhline(0,color="black",lw=1);ax.axvline(root,color="red",linestyle="--")
ax.set(xlabel="x",ylabel="5(1−exp(−x))−x",title="Trivial zero and the positive Wien root");plt.tight_layout();plt.show()

Positive root: 4.9651142317961785 ; independent brentq: 4.965114231744276

	iteration	lower	upper	midpoint	f_midpoint	width
0	0	1.0000	10.000000	5.500000	-0.520434	9.000000
1	1	1.0000	5.500000	3.250000	1.556129	4.500000
2	2	3.2500	5.500000	4.375000	0.562059	2.250000
3	3	4.3750	5.500000	4.937500	0.026637	1.125000
4	4	4.9375	5.500000	5.218750	-0.245820	0.562500
5	5	4.9375	5.218750	5.078125	-0.109283	0.281250
6	6	4.9375	5.078125	5.007812	-0.041240	0.140625
7	7	4.9375	5.007812	4.972656	-0.007280	0.070312

2. 检查区间证书

每次迭代保留两端异号区间，宽度逐次减半。终止条件将中点位置误差约束为当前区间半宽。对于平坦函数，仅有较小残差不保证同等位置精度。

In [5] · Python

fig,ax=plt.subplots(figsize=(7,4));ax.semilogy(trace.iteration,trace.width/2,label="Certified half-width")
ax.semilogy(trace.iteration,np.maximum(np.abs(trace.midpoint-reference),1e-16),label="Observed error to reference")
ax.set(xlabel="Iteration",ylabel="x error bound / observed error",title="Bisection convergence certificate");ax.legend();plt.tight_layout();plt.show()

3. 将无量纲根转换为物理量

使用冻结的 h、c、k_B，b=hc/(k_B x) 的单位为米·开尔文。曲线展示假定黑体温度 K 对应的峰值波长 nm。这些是理想模型值，不是恒星实测光谱。

In [6] · Python

b=constants["h"]["value"]*constants["c"]["value"]/(constants["k_B"]["value"]*root)
temperatures=np.linspace(2500,10000,101);peak_nm=1e9*b/temperatures
fig,ax=plt.subplots(figsize=(7,4));ax.plot(temperatures,peak_nm)
ax.set(xlabel="Assumed temperature (K)",ylabel="Peak wavelength (nm)",title="Wien displacement derived from a verified root");plt.tight_layout();plt.show()
print("Derived displacement constant (m K):",b)

Derived displacement constant (m K): 0.002897771955154881

4. 检查符号、容差及导出单位

将根与参考求解器核对，并检查全部记录区间。不包含变号的区间必须报错。数值容差是计算设置，不是实验测量不确定性。

In [7] · Python

assert abs(root-reference)<=trace.width.iloc[-1]/2+1e-12
assert all(equation(lo)*equation(hi)<=0 for lo,hi in zip(trace.lower,trace.upper))
assert np.isclose(b,2.897771955e-3,rtol=1e-8)
try:bisection(equation,6,10)
except ValueError:print("Non-bracketing interval rejected")
else:raise AssertionError("Missing sign change accepted")
OUTPUT_DIR.mkdir(parents=True,exist_ok=True);trace.to_csv(OUTPUT_DIR/"bisection-trace.csv",index=False)
pd.DataFrame({"temperature_K":temperatures,"peak_wavelength_nm":peak_nm}).to_csv(OUTPUT_DIR/"wien-curve.csv",index=False)
print("Root certificate and unit checks passed; exports:",OUTPUT_DIR)

Non-bracketing interval rejected
Root certificate and unit checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ma02-bracketed-roots

自主练习与参考答案

为何求维恩谱峰时排除 x=0？
区间减半能证实什么？

参考解释。 平凡根不是正的内部谱峰。对于保留的连续函数变号区间，半宽约束中点位置误差。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

仅用 NumPy、SciPy、Matplotlib 比较 1e-6、1e-8、1e-10 容差，保留全部区间证书。

提示词 2

使用同样的包标注 3000、5000、8000 K 的预测峰值波长及单位。

人工验收。 各保留区间变号，与 brentq 的根误差不超过半宽，波长为正。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [8] · Python

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://physics.nist.gov/cuu/Constants/',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged

来源、快照与下一步

NIST CODATA physical constants — NIST public reference data; credit NIST/CODATA; retrieved 2026-10-09. SHA-256 77fb90e66c40db3e6eb16630bc9c88e4c7c8beddbe5e71be406f2f26e3f67e67. Preserve official ASCII table; extract seven explicitly named SI constants and their units.

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。