Check a heat-diffusion solver

Enforce the explicit stability limit, compare a known sine-mode solution and measure grid convergence.

Physics · Lesson 04🔋 Offline core analysis⏱ 50 minNumerical teaching model📦 numpy · pandas · scipy · matplotlib · IPython

Goal and how to run

This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.

Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.

Setup and parameters

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 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='ph04-heat'
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

Load the embedded snapshot

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={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['models']

Read the method functions

In [3] · Python

def heat_diffusion_explicit(nodes=51,length_m=1,alpha=1e-4,duration_s=600,ratio=.4,ambient=20,amplitude=50):
    """FTCS heat equation with fixed-temperature ends and a known sine-mode initial condition."""
    if nodes<3 or length_m<=0 or alpha<=0 or duration_s<=0 or not 0<ratio<=.5:
        raise ValueError('Require >=3 nodes, positive parameters and FTCS ratio <=0.5')
    x=np.linspace(0,length_m,nodes);spacing=x[1]-x[0]
    count=int(np.ceil(duration_s/(ratio*spacing**2/alpha)));dt=duration_s/count;actual_ratio=alpha*dt/spacing**2
    times=np.linspace(0,duration_s,count+1);temperature=np.empty((count+1,nodes))
    temperature[0]=ambient+amplitude*np.sin(np.pi*x/length_m);temperature[:,0]=ambient;temperature[:,-1]=ambient
    for j in range(count):
        temperature[j+1,1:-1]=temperature[j,1:-1]+actual_ratio*(temperature[j,2:]-2*temperature[j,1:-1]+temperature[j,:-2])
    return dict(x=x,t=times,temperature=temperature,spacing_m=spacing,step_s=dt,ratio=actual_ratio)

print("Method ready: heat_diffusion_explicit")

Method ready: heat_diffusion_explicit

Data and model boundary

Pendulum, oscillator, projectile and heat trajectories are declared numerical models, not laboratory measurements. The gravitational acceleration is the conventional standard g0 when used. GW150914 uses actual frozen detector strain and published event timing; a filtered curve is not an independent significance or source-parameter analysis. Blackbody curves assume an ideal emitter.

Steps

1. Respect the explicit stability limit

For FTCS, r=alpha·dt/dx² must not exceed 0.5 in one dimension. End temperatures are fixed at the ambient value. The initial sine profile has a closed solution and provides a meaningful solver reference.

In [4] · Python

p=snapshot("models")["heat"]
solution=heat_diffusion_explicit(51,p["length_m"],p["diffusivity_m2_per_s"],600,.4,p["ambient_degC"],p["amplitude_degC"])
print("Actual dx (m):",solution["spacing_m"],"; dt (s):",solution["step_s"],"; r:",solution["ratio"])
fig,ax=plt.subplots(figsize=(8,4))
for time in [0,150,300,600]:
 index=int(np.argmin(abs(solution["t"]-time)));ax.plot(solution["x"],solution["temperature"][index],label=f"t={solution['t'][index]:.1f} s")
ax.set(xlabel="Position (m)",ylabel="Model temperature (°C)",title="Fixed-temperature ends and diffusing sine profile");ax.legend();plt.tight_layout();plt.show()

Actual dx (m): 0.02 ; dt (s): 1.6 ; r: 0.4

2. View the computed temperature field

Color is temperature in °C; the horizontal coordinate is metres and vertical coordinate is seconds. Cooling occurs because fixed end reservoirs remove heat. An isolated-system energy-conservation claim would use different boundary conditions.

In [5] · Python

fig,ax=plt.subplots(figsize=(8,5));im=ax.pcolormesh(solution["x"],solution["t"],solution["temperature"],shading="auto",cmap="inferno")
ax.set(xlabel="Position (m)",ylabel="Model time (s)",title="Explicit heat-equation solution");fig.colorbar(im,ax=ax,label="Temperature (°C)");plt.tight_layout();plt.show()

3. Measure joint space/time grid convergence

Keep r fixed while refining x, which also refines dt. Compare the final profile with T0+A sin(pi x/L) exp(−alpha(pi/L)²t). This controlled model benchmark measures numerical error, not material-property uncertainty.

In [6] · Python

nodes=[21,41,81,161];records=[]
for count in nodes:
 trial=heat_diffusion_explicit(count,p["length_m"],p["diffusivity_m2_per_s"],600,.4,p["ambient_degC"],p["amplitude_degC"])
 exact=p["ambient_degC"]+p["amplitude_degC"]*np.sin(np.pi*trial["x"]/p["length_m"])*np.exp(-p["diffusivity_m2_per_s"]*(np.pi/p["length_m"])**2*600)
 records.append({"nodes":count,"dx_m":trial["spacing_m"],"dt_s":trial["step_s"],"max_error_degC":float(np.max(abs(trial["temperature"][-1]-exact)))})
convergence=pd.DataFrame(records)
fig,ax=plt.subplots(figsize=(7,4));ax.loglog(convergence.dx_m,convergence.max_error_degC,"o-")
ax.set(xlabel="Spatial spacing (m)",ylabel="Maximum final error (°C)",title="Known-solution grid audit");plt.tight_layout();plt.show()

4. Test invalid stability settings and export

The implementation rejects r>0.5. Check fixed ends, bounded temperatures and decreasing refinement error. Record both requested r and actual dt so rounding the number of steps does not hide the actual scheme.

In [7] · Python

try:heat_diffusion_explicit(ratio=.7)
except ValueError:print("Unstable FTCS ratio rejected")
else:raise AssertionError("Unstable ratio accepted")
assert np.allclose(solution["temperature"][:,[0,-1]],p["ambient_degC"])
assert solution["temperature"].min()>=p["ambient_degC"]-1e-9
assert np.all(np.diff(convergence.max_error_degC)<0)
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame(solution["temperature"],index=solution["t"],columns=solution["x"]).to_csv(OUTPUT_DIR/"temperature-field.csv")
convergence.to_csv(OUTPUT_DIR/"grid-convergence.csv",index=False)
print("Boundary, stability and refinement checks passed; exports:",OUTPUT_DIR)

Unstable FTCS ratio rejected
Boundary, stability and refinement checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ph04-heat

Exercises and reference explanation

Explain the dependence of allowable dt on dx.
Why does integrated thermal excess decrease with these boundaries?

Reference explanation. The stable dt scales with dx²/alpha. Fixed-temperature endpoints exchange heat with reservoirs, so the system is not isolated.

Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.

Optional AI coding exercises

Prompt 1

Using NumPy/Pandas/Matplotlib, compare r=.2 and .4 at the same spatial grid, recording actual dt.

Prompt 2

With the same packages, add final analytic and numerical profiles plus an error subplot.

Human acceptance. Reject r>.5, preserve fixed ends and verify refinement against the sine-mode solution.

Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.

Optional online extension

In [8] · Python

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.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

Sources, snapshots and next steps

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.

Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.