codexangstrom-periodic-heating-amplitude-phase.py
angstrom-periodic-heating-amplitude-phase.py
"""Original solved rod example, with independent complex finite-difference checks.
Python 3.14, NumPy 2.3.5, matplotlib 3.10.7; no SciPy.
Run in _media; writes its basename PNG to the current working directory.
"""
from pathlib import Path
import json
import matplotlib.pyplot as plt
import numpy as np
ALPHA = 1e-4
PERIOD = 100.0
OMEGA = 2*np.pi/PERIOD
LOSS = 0.01
THETA0 = 4.0
SENSORS = np.array([0.02, 0.05])
def gamma(loss):
return np.sqrt((loss + 1j*OMEGA)/ALPHA)
def finite_rod(x, length, loss=LOSS):
"""Exact phasor: prescribed left temperature, insulated right end."""
g = gamma(loss)
return THETA0*np.cosh(g*(length-x))/np.cosh(g*length)
def finite_difference(length, dx, loss, insulated=False):
"""Solve discretized αz'' − (loss+iω)z=0 by complex Thomas elimination.
This assembles the differential equation, without using its analytic roots.
The final node is eliminated via z_N = z_(N-1) for insulation at the
half-grid point L-dx/2; otherwise the right temperature is zero at L.
"""
n = round(length/dx)
x = np.arange(n+1)*dx
m = n-1
diag = np.full(m, 2*ALPHA/dx**2 + loss + 1j*OMEGA, dtype=complex)
off = -ALPHA/dx**2
rhs = np.zeros(m, dtype=complex)
rhs[0] = -off*THETA0
if insulated:
diag[-1] += off
for j in range(1, m):
factor = off/diag[j-1]
diag[j] -= factor*off
rhs[j] -= factor*rhs[j-1]
z = np.zeros(n+1, dtype=complex)
z[0] = THETA0
z[n-1] = rhs[-1]/diag[-1]
for j in range(m-2, -1, -1):
z[j+1] = (rhs[j]-off*z[j+2])/diag[j]
if insulated:
z[-1] = z[-2]
return x, z
def infer(z, spacing=0.03):
attenuation = np.log(abs(z[0]/z[1]))
lag = np.angle(z[0]/z[1])
return OMEGA*spacing**2/(2*attenuation*lag)
def checks():
out = {"parameters": {"alpha": ALPHA, "period_s": PERIOD, "loss_rate_s_inverse": LOSS}, "cases": []}
for loss in [0.0, LOSS]:
g = gamma(loss)
exact = THETA0*np.exp(-g*SENSORS)
errors = []
for dx in [0.001, 0.0005, 0.00025]:
x, z = finite_difference(0.8, dx, loss)
obs = z[np.rint(SENSORS/dx).astype(int)]
errors.append(float(np.max(abs(obs/exact-1))))
assert errors[-1] < 3e-6
assert errors[0]/errors[1] > 3.9 and errors[1]/errors[2] > 3.9
# Recover phasors from real time traces by an independent least-squares fit.
t = np.linspace(0, 3*PERIOD, 1500, endpoint=False)
design = np.column_stack([np.ones_like(t), np.cos(OMEGA*t), np.sin(OMEGA*t)])
trace = 0.7 + (exact[:, None]*np.exp(1j*OMEGA*t)).real
fit = np.linalg.lstsq(design, trace.T, rcond=None)[0]
recovered = fit[1]-1j*fit[2]
assert np.max(abs(recovered-exact)) < 1e-13
assert abs(infer(recovered)/ALPHA-1) < 1e-13
assert abs(2*g.real*g.imag-OMEGA/ALPHA) < 1e-12
out["cases"].append({"loss_rate": loss, "a": float(g.real), "b": float(g.imag), "amplitudes_K": abs(exact).tolist(), "lag_rad": float(g.imag*0.03), "time_lag_s": float(g.imag*0.03/OMEGA), "inferred_alpha": float(infer(recovered)), "finite_difference_relative_errors": errors, "refinement_ratios": [errors[0]/errors[1], errors[1]/errors[2]], "trace_fit_max_error_K": float(np.max(abs(recovered-exact)))})
# Finite-end check uses a different boundary condition and exact cosh solution.
finite_errors = []
for dx in [0.0004, 0.0002, 0.0001]:
x, z = finite_difference(0.08, dx, LOSS, insulated=True)
exact = finite_rod(x[:-1], 0.08-dx/2)
finite_errors.append(float(np.max(abs(z[:-1]-exact))/THETA0))
assert finite_errors[-1] < 3e-7
assert finite_errors[0]/finite_errors[1] > 3.8
out["finite_insulated_end"] = {"finite_difference_max_relative_to_drive_errors": finite_errors, "alpha_apparent_L_0_08": float(infer(finite_rod(SENSORS, 0.08))), "alpha_apparent_L_0_3": float(infer(finite_rod(SENSORS, 0.3)))}
return out
def main():
checked = checks()
plt.rcParams.update({"font.size": 10, "figure.facecolor": "white"})
fig, axes = plt.subplots(2, 2, figsize=(10.5, 7.6), dpi=120, constrained_layout=True)
x = np.linspace(0, 0.12, 400)
colors = ["#146c94", "#b34f23"]
for loss, col in zip([0, LOSS], colors):
g = gamma(loss)
label = "No side loss" if loss == 0 else "Linear loss λ = 0.010 s⁻¹"
axes[0, 0].plot(x*100, THETA0*np.exp(-g.real*x), color=col, label=label)
axes[0, 1].plot(x*100, g.imag*x, color=col, label=label)
xn, zn = finite_difference(0.8, 0.0005, loss)
idx = np.rint(SENSORS/0.0005).astype(int)
axes[0, 0].plot(SENSORS*100, abs(zn[idx]), "o", color=col, ms=5)
axes[0, 1].plot(SENSORS*100, -np.angle(zn[idx]), "o", color=col, ms=5)
axes[0, 0].set(title="Spatial amplitude: loss makes decay steeper", xlabel="Distance x (cm)", ylabel="Amplitude A(x) (K)")
axes[0, 0].legend(fontsize=9)
axes[0, 1].set(title="Unwrapped lag: loss makes this slope shallower", xlabel="Distance x (cm)", ylabel="Lag from driven end (rad)")
for ax in axes[0]:
for sx in SENSORS*100:
ax.axvline(sx, color="#999999", ls=":", lw=0.8)
t = np.linspace(0, 2*PERIOD, 600)
for sx, col in zip(SENSORS, ["#146c94", "#7c3f96"]):
z = THETA0*np.exp(-gamma(LOSS)*sx)
axes[1, 0].plot(t, (z*np.exp(1j*OMEGA*t)).real, color=col, label=f"x = {sx*100:.0f} cm, A = {abs(z):.3f} K")
axes[1, 0].set(title="Steady periodic traces with linear loss", xlabel="Time t (s)", ylabel="Oscillating temperature θ (K)")
axes[1, 0].legend(fontsize=9, loc="lower right")
axes[1, 0].text(0.04, 0.95, "Δφ = 0.4912 rad; Δt = 7.818 s", transform=axes[1, 0].transAxes, va="top", fontsize=9)
lengths = np.linspace(0.06, 0.32, 300)
apparent = np.array([infer(finite_rod(SENSORS, length)) for length in lengths])
axes[1, 1].plot(lengths*100, apparent/ALPHA, color="#35754c")
axes[1, 1].axhline(1, ls="--", color="#444444", lw=1)
axes[1, 1].plot([8, 30], [infer(finite_rod(SENSORS, length))/ALPHA for length in [0.08, 0.3]], "o", color="#35754c")
axes[1, 1].set(title="Insulated far end: bias from using one exponential", xlabel="Total rod length L (cm)", ylabel="Apparent α / true α", ylim=(0.7, 1.7))
for ax in axes.flat:
ax.grid(alpha=0.2)
fig.suptitle("Original model: α = 10⁻⁴ m²/s, period = 100 s, driven amplitude = 4 K", fontsize=13)
fig.savefig(Path(__file__).stem + ".png", facecolor="white")
plt.close(fig)
print(json.dumps(checked, indent=2))
if __name__ == "__main__":
main()