测量数值微分与积分误差
将高斯曲线表、有限差分和梯形积分与解析结果比较。
目标与运行方法
本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。
学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。
环境与参数
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
import math
from scipy.integrate import quad
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='ma01-numerical-calculus'
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
读取内嵌快照
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={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['models']
方法函数
def trapezoid_integral(x,y):
"""Integrate a tabulated curve with explicitly increasing abscissae."""
x=np.asarray(x,dtype=float);y=np.asarray(y,dtype=float)
if x.ndim!=1 or x.shape!=y.shape or len(x)<2 or not np.isfinite(x+y).all():raise ValueError('Invalid curve')
widths=np.diff(x)
if np.any(widths<=0):raise ValueError('x must be strictly increasing')
return float(np.sum(widths*(y[:-1]+y[1:])/2))
print("Method ready: trapezoid_integral")
Method ready: trapezoid_integral
数据与模型边界
解析函数与 Monte Carlo 点属于数学教学模型，不是实验测量。使用 Iris 观测时，可交换性和采样设计仍是统计假设；方便样本不自动代表总体。数值离散误差与采样不确定性应分别解释。
分步分析
1. 声明函数与参考积分
使用 [−2,2] 区间上的 f(x)=exp(−x²)。本基准的 x 与函数值无量纲。有限区间积分为 sqrt(pi)·erf(2)，不能与无限区间积分 sqrt(pi) 混淆。
model_inputs=snapshot("models");seed=model_inputs["seed"]
LOWER,UPPER=-2.0,2.0;POINTS=161
x=np.linspace(LOWER,UPPER,POINTS);y=np.exp(-x*x)
reference=math.sqrt(math.pi)*math.erf(2)
quad_reference,quad_error=quad(lambda value:np.exp(-value*value),LOWER,UPPER)
assert np.isclose(reference,quad_reference,rtol=1e-12)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(x,y,label="exp(−x²)");ax.fill_between(x,y,alpha=.2)
ax.set(xlabel="x (dimensionless)",ylabel="f(x)",title="Finite-interval integration problem");ax.legend();plt.tight_layout();plt.show()
2. 对曲线表求导
平滑函数的内部中心差分为二阶误差。端点公式可能具有不同误差阶，因此单独计算内部误差。np.gradient 使用真实 x 网格，不假定间距为一。
derivative=np.gradient(y,x,edge_order=2);exact_derivative=-2*x*np.exp(-x*x)
fig,ax=plt.subplots(figsize=(7,4));ax.plot(x,exact_derivative,label="Analytic",lw=2);ax.plot(x,derivative,"--",label="Tabulated finite difference")
ax.set(xlabel="x",ylabel="df/dx",title="Derivative from explicit grid spacing");ax.legend();plt.tight_layout();plt.show()
print("Maximum interior derivative error:",float(np.max(np.abs(derivative[1:-1]-exact_derivative[1:-1]))))
Maximum interior derivative error: 0.00040651967031690717
3. 细化积分网格
舍入误差较小时，梯形积分误差应近似按 h² 下降。收敛图双轴取对数。这里的误差描述已知平滑模型，不是有噪声测量的置信区间。
intervals=np.array([20,40,80,160,320]);estimates=[];widths=[]
for count in intervals:
grid=np.linspace(LOWER,UPPER,count+1);widths.append((UPPER-LOWER)/count)
estimates.append(trapezoid_integral(grid,np.exp(-grid*grid)))
error=np.abs(np.array(estimates)-reference)
fig,ax=plt.subplots(figsize=(7,4));ax.loglog(widths,error,"o-",label="Observed absolute error")
ax.loglog(widths,error[0]*(np.array(widths)/widths[0])**2,"--",label="h² comparison")
ax.set(xlabel="Grid spacing h",ylabel="Absolute integral error",title="Trapezoid grid refinement");ax.legend();plt.tight_layout();plt.show()
4. 核对并导出数值审计
将可读实现与 NumPy 梯形规则对照，并检查网格顺序约束。解析值、求积容差和离散误差分别记录。
assert np.isclose(trapezoid_integral(x,y),np.trapezoid(y,x))
assert error[-1]<error[0] and np.all(error>0)
try:trapezoid_integral(x[::-1],y[::-1])
except ValueError:print("Reversed abscissae rejected")
else:raise AssertionError("Invalid grid accepted")
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"x":x,"f":y,"finite_difference":derivative,"analytic_derivative":exact_derivative}).to_csv(OUTPUT_DIR/"curve-and-derivative.csv",index=False)
pd.DataFrame({"intervals":intervals,"h":widths,"integral":estimates,"absolute_error":error}).to_csv(OUTPUT_DIR/"grid-refinement.csv",index=False)
(OUTPUT_DIR/"references.json").write_text(json.dumps({"finite_interval":[LOWER,UPPER],"analytic_integral":reference,"quad_reference":quad_reference,"quad_error_estimate":quad_error},indent=2))
print("Reference checks passed; exports:",OUTPUT_DIR)
Reversed abscissae rejected Reference checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ma01-numerical-calculus
自主练习与参考答案
- 说明高斯函数有限区间和无限区间积分的区别。
- 比较内部与端点微分误差。
参考解释。 [−2,2] 外的尾部被排除。内部中心差分与端点公式需要分别分析。
先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。
可选的 AI 编程练习
提示词 1
仅用 NumPy、SciPy、Matplotlib 将每个积分网格加密一倍，报告误差比并保持解析参考值。
提示词 2
使用同样的包比较 edge_order=1、2，分别画端点和内部误差。
人工验收。 核对严格递增网格、NumPy 梯形规则及有限区间 erf 解析表达式。
检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。
可选联网扩展
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://docs.scipy.org/doc/scipy/tutorial/integrate.html',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
来源、快照与下一步
IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256
a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256
4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.
下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。