Derive the heat time-step restriction

Try

Approximate heat on a uniform grid with spacing h and time step τ. Let Uⱼᵐ represent the value at position jh and time mτ. Forward time and centered space give Uⱼᵐ⁺¹=Uⱼᵐ+r(Uⱼ₋₁ᵐ−2Uⱼᵐ+Uⱼ₊₁ᵐ), where r=κτ/h².

This is a discrete model of the PDE. It introduces approximation error even with exact arithmetic, and its stability depends on the relation between the two step sizes.

Predict

If h is halved while τ stays fixed, what happens to r? Predict whether finer spatial resolution alone can make an explicit computation unstable. This is a useful counterexample to the idea that any refinement automatically improves a numerical solution.

Experiment

First predict the result for a rod with fixed endpoint values 0 and 1. All variables here are dimensionless and κ = 1. With u(x,0) = x + sin(πx), the exact answer is u(x,t) = x + exp(−π²t) sin(πx). Subtracting the affine lift x gives a zero-boundary problem. The initial maximum exceeds 1, so the maximum principle bounds the solution by the initial and boundary extrema, not by the interval [0,1].

Predict, run, retain and compare

Dimensionless position x, time t and value u; diffusivity κ = 1. Endpoints stay 0 and 1. The initial profile can exceed 1: the maximum principle uses its initial extrema, not just the boundary values. Controls prepare the next run; the plot always labels the submitted settings.

Intervals N81632 Time step1/1281/5121/2048 Initial profileFixed initial data: x + sin(πx)Grid-dependent probe: smooth + 2⁻¹⁰ high modeGrid-dependent probe: x + high mode Maximum committed steps014832128

Run and retain evidence (0/12) Reset controls and preview Download retained JSON

Records are temporary; download before reload or leaving. Reset retains records. The initial preview is not a retained run. The cap is 12 runs, 128 steps and 3,968 interior updates per run. A rejected candidate above |u| = 1,000,000 retains the preceding grid.

Copy retained JSON evidence

This is exactly the report used by Download. Select and copy it into an external file before leaving.

Retained heat-grid report{ "version": "heat-grid-investigation-v1", "limitations": "Synthetic dimensionless teaching model. Floating arithmetic; exact-arithmetic bounds exclude rounding. Supplied runs are not independent learner implementation.", "records": [] }

Submitted: N = 8, τ = 1/128, profile smooth, budget 128. Status: target-reached; completed 8/8 steps toward T = 1/16; attempted 8.

r = 0.5000000; maximum-norm positivity certificate: available. Fixed-grid spectral threshold 0.5197831: stable. These are different claims.

Inspect committed step 0: t = 0.0000001: t = 0.0078125002: t = 0.015625003: t = 0.023437504: t = 0.031250005: t = 0.039062506: t = 0.046875007: t = 0.054687508: t = 0.06250000

Solid red: continuous PDE sampled at grid points. Dashed navy: exact-arithmetic discrete-mode formula evaluated in floating point. Dots: actual stencil. Lines join samples; they do not resolve between-grid behavior. All sampled values fit the shared vertical range.

Step 8, t = 0.06250000: stencil minimum 0.000000, maximum 1.125325. Observed maximum PDE error 0.008851443; stencil versus discrete formula 1.110223e-16; discrete versus continuous formula 0.008851443.

Exact-arithmetic discretization bound: 0.03170869. This rounded display excludes floating rounding and physical modeling error. The formula comparison is not a rigorous roundoff bound because both calculations use floating arithmetic.

Interior updates: 56 committed, 0 rejected, 56 attempted. Stencil arithmetic: 168 products and 112 additions; excludes oracle, plotting, allocation and verification work.

Mode factors: k=1, amplitude=1, G=0.9238795.

Evidence challenge

Compare the same smooth data at (N, τ) = (8,1/128), (16,1/512), (32,1/2048). Then hold τ = 1/128 while refining N. Predict before adding the small high-mode probe. The high-mode profiles depend on N: they test robustness, not convergence for one fixed initial function. Retain each manifest, proof, numerical trace and revision in the written investigation below.

Rewrite the update as a weighted average with weights r, 1−2r and r.

r	Neighbor weights	Center weight	Convex average?
1/4	1/4 each	1/2	Yes
1/2	1/2 each	0	Yes
3/4	3/4 each	−1/2	No

For a center value 1 and zero neighbors, the last row produces −1/2. The exact homogeneous heat equation with nonnegative data cannot do that. Calculate a second step for an alternating grid profile to see why this is an amplification issue, not just an isolated sign surprise.

Understand

On a periodic or infinite grid, substitute a Fourier mode Uⱼᵐ=Gᵐe^(ijθ). The amplification factor is G=1+r(e^(−iθ)−2+e^(iθ))=1−4r sin²(θ/2). Requiring its modulus at most 1 for all θ gives 0≤r≤1/2. The alternating highest-frequency mode θ=π is the limiting case. At r=1/2 it changes sign without decaying; stability permits that neutral behavior even though exact heat damps every nonconstant mode.

For smooth solutions, Taylor expansion gives a consistency residual O(τ+h²) after dividing the update error by τ. The one-step defect itself is O(τ²+τh²). Confusing these quantities leads to incorrect order statements. Stability, compatible boundary treatment and regularity estimates are also needed to turn consistency into convergence. Under standard linear well-posed settings, a stability-plus-consistency theorem supplies that connection; it is not a universal rule for every nonlinear scheme.

Backward Euler replaces the spatial values by their new-time values. Its factor is G=1/[1+4r sin²(θ/2)], so it is stable for every nonnegative r. It requires solving a coupled linear system at each step. Large steps can still misrepresent decay and timing: unconditional stability does not imply unconditional accuracy.

Prove the finite-boundary claim

For N intervals on [0,1], h = 1/N and interior indices j = 1,…,N−1. Keep both endpoints exact at every step. The error therefore has zero endpoints. If 0 ≤ r ≤ 1/2, the three update weights are nonnegative and sum to one. Each new error magnitude is at most the previous maximum magnitude. This proves maximum-norm nonexpansiveness directly, including the rows adjacent to a boundary. It also proves preservation of the initial-and-boundary range for the numerical solution.

The finite-grid modes are sin(kπj/N), k = 1,…,N−1. Substitution gives Gₖ = 1 − 4r sin²(kπ/(2N)). These modes form an orthogonal basis of the interior grid: the discrete sine orthogonality identity follows by summing the finite geometric series for complex exponentials and taking real parts. Thus the update is diagonal in an orthogonal basis. Every mode is nonamplifying exactly when r ≤ 1/[2cos²(π/(2N))]. This is the fixed-grid Euclidean-norm spectral condition. It is slightly weaker than r ≤ 1/2; the stronger condition also supplies positivity and the maximum-norm proof above. For N = 8, r = 0.51 is a useful counterexample to “every r above 1/2 is spectrally unstable.” No single-grid condition alone proves convergence across a family of grids.

Accumulate the local defect instead of quoting convergence

Assume a classical solution with continuous second time and fourth spatial derivatives on the closed time-space rectangle, exact initial samples and compatible exact boundary values. Let Mtt and Mxxxx bound the magnitudes of those derivatives. Taylor’s theorem in time bounds the forward step remainder by τ²Mtt/2. The centered second difference approximates uxx with error at most h²Mxxxx/12. Since ut = uxx, the one-step defect has magnitude at most d = τ²Mtt/2 + τh²Mxxxx/12. This is a bound before division by τ.

Let Eₙ be the largest interior absolute error after n steps. Subtract the exact sampled equation from the numerical update. The convex-weight argument gives Eₙ₊₁ ≤ Eₙ + d. Induction gives Eₙ ≤ E₀ + nd = E₀ + tₙ(τMtt/2 + h²Mxxxx/12). This proves convergence for a fixed smooth initial-boundary problem when both steps shrink and r stays at most 1/2. Nonzero initial error adds E₀; inaccurate boundary values need additional estimates and are excluded here.

For the finite sine profiles in the workbench, Mtt and Mxxxx are bounded by π⁴ times the sum of |aₖ|k⁴, where aₖ is each sine amplitude. Exponential decay is at most one for nonnegative time, and the affine lift contributes zero to these derivatives. The workbench displays this exact-arithmetic bound rounded to floating precision. It is not a guarantee about accumulated floating-point error, and the stencil/discrete-oracle difference includes rounding in both computations.

At T = 1/16, the smooth pairs (N,τ) = (8,1/128), (16,1/512), (32,1/2048) give observed maximum errors approximately 0.00885144, 0.00215684 and 0.000535841. Their respective proved discretization bounds are approximately 0.0317087, 0.00792717 and 0.00198179. Error reduction agrees with the joint refinement argument; numerical agreement does not replace that argument.

Repair a misleading refinement experiment

Hold τ = 1/128 and change N from 8 to 16. Now r changes from 1/2 to 2. The smooth low mode can still look reasonable: its final observed error is about 0.0122358. Add the grid’s highest mode with amplitude 2⁻¹⁰. Its factor at N = 16 is approximately −6.92314, and after eight steps the observed maximum PDE error exceeds 5,153. The negative factor alternates signs while magnitudes grow. At N = 32 the same time step can trigger the resource ceiling before T; that incomplete run is not a final-time error measurement.

The high-mode probe changes its frequency with N. It is a stability stress test, not a convergence comparison for one fixed continuous initial function. To repair the smooth-data convergence experiment, reduce τ by four when doubling N and retain identical smooth initial data. To repair the robustness report, retain the changed profile explicitly and analyze its mode factor rather than hiding the perturbation or discarded steps.

Apply

Choose κ and h, compute the largest stable explicit τ, then halve h and recompute. Explain the resulting increase in time-step count over a fixed duration. Compare an implicit large-step result with an explicit refined result and an exact sine mode, separating spatial truncation, time truncation and floating-point error.

Continue

Transport has a different restriction, proportional to h rather than h². Following the direction of characteristics explains upwinding and why a small step alone cannot rescue every stencil.