Butterfly Effect Simulation: Weather Predictability

Explore the butterfly effect with two nearly identical Lorenz-63 trajectories. Animate their divergence and learn what the model reveals about weather predictability.

#butterfly effect #Lorenz attractor #numerical convergence #Taylor method #weather predictability

A butterfly effect simulation makes a difficult forecasting idea visible: two states that are almost indistinguishable now can develop very different futures under exactly the same equations. This calculator runs the Lorenz 1963 convection model from two initial states and measures when their separation reaches a threshold you choose. It is intended for engineers, students, teachers, and anyone exploring why weather can't be predicted indefinitely in precise detail.

Start with the defaults. Run A begins at (1, 1, 1); Run B differs only by 0.001 in its x coordinate. Press play or drag the playback position to follow the coloured points around the Lorenz attractor. Blue represents Run A and orange represents Run B. Grey gives the numerically agreed portion of Run A as context; coloured trails show the preceding three model-time units. The drawing is an orthogonal projection of three-dimensional state space, not a map of air motion.

The first result is the earliest crossing of the chosen full three-dimensional separation threshold. A second integration at a different Taylor order checks the trajectory and crossing. A greater-than result is shown only when both integrations agree throughout the simulation window and neither finds a crossing. If agreement is lost, the relevant result reads numerically unresolved instead of claiming a lower bound. Move playback earlier or shorten the window to examine a resolved part of the experiment. The trajectory agreement card identifies the available span.

The cursor separation uses the ODE locally, rather than a straight line between plotted distances. Charts omit the trajectory tail after refinement agreement is lost, and the animated states stop at the last agreed time. Identical starts remain an exact zero-separation control, even if their late-time coordinates cannot be resolved.

Use the logarithmic separation chart to compare tiny differences with larger ones. A vertical increase of one means ten times the separation. Error can shrink temporarily, so separation need not increase steadily. The second chart compares the two x histories. Change the initial state, perturbation, threshold, or window to explore the model.

The dimensionless model time has no calibrated conversion to days. Zhang and colleagues' 2019 global-model study supports a roughly two-week horizon for detailed midlatitude weather under its experimental assumptions. That estimate does not come from this three-variable model and is not a universal expiry date for forecasts. Probabilistic outlooks, large-scale averages, and slower climate signals can retain useful skill for longer.

This tool excludes geography, moisture, observations, data assimilation, model error, and forecast probabilities. Numerical refinement agreement is an accuracy check, not a rigorous error bound or meteorological confidence level. Use the calculator to study sensitivity and numerical dynamics, rather than judge a particular weather forecast.

How it is calculated

Lorenz (1963), equations (25)–(27), printed p. 135, supplies the dimensionless convection equations:

x˙=σ(y−x),y˙=x(ρ−z)−y,z˙=xy−βz.\dot{x}=\sigma(y-x),\qquad\dot{y}=x(\rho-z)-y,\qquad\dot{z}=xy-\beta z.

The classic parameters are σ=10\sigma=10, ρ=28\rho=28, β=8/3\beta=8/3 (printed p. 136). Run A starts at (x0,y0,z0)(x_0,y_0,z_0) and Run B at (x0+δ0,y0,z0)(x_0+\delta_0,y_0,z_0).

The library's fixed-tolerance ode45 is insufficient for the long-time benchmarks. The tool therefore integrates each run using an explicit order-24 Taylor method, with an independently evolved order-20 comparison. The recurrence is derived here by substituting power series into the published ODE. For coefficients already scaled by the candidate step hh:

ak+1=hσk+1(bk−ak),a_{k+1}=\frac{h\sigma}{k+1}(b_k-a_k),
bk+1=hk+1(ρak−bk−∑j=0kajck−j),b_{k+1}=\frac{h}{k+1}\left(\rho a_k-b_k-\sum_{j=0}^k a_jc_{k-j}\right),
ck+1=hk+1(∑j=0kajbk−j−βck).c_{k+1}=\frac{h}{k+1}\left(\sum_{j=0}^k a_jb_{k-j}-\beta c_k\right).

Candidate steps are at most 0.1 model-time units. The last two retained terms control step size: targets are 10−14(1+max⁡∣state∣)10^{-14}(1+\max|state|) for order 24 and 10−11(1+max⁡∣state∣)10^{-11}(1+\max|state|) for order 20. Both orders advance independently on the accepted shared steps. The polynomial gives dense output on a grid spaced by at most 0.02. No curve-fitting calibration to reference answers is used.

The distance is δ(t)=(xB−xA)2+(yB−yA)2+(zB−zA)2\delta(t)=\sqrt{(x_B-x_A)^2+(y_B-y_A)^2+(z_B-z_A)^2}. The first sampled crossing brackets an event; an order-18 local Taylor expansion and bisection refine its time. The cursor uses the same local ODE evaluation.

Numerical agreement requires the 3D coordinate gap between orders to remain at most 0.001 up to the reported time. The two crossing statuses must agree, and crossing times must differ by at most 0.001. A no-crossing lower bound additionally needs agreement through the whole window and a threshold margin exceeding twice the largest refinement gap plus 0.001. Failing results are withheld as numerically unresolved. This is an agreement test, not rigorous interval arithmetic.

The scene uses ((x+y)/2,z)((x+y)/\sqrt{2},z) at equal display scale. Grey context and charts stop at the first loss of trajectory agreement. The log chart clips separation below 10−1210^{-12} for display only.

Assumptions

Both runs use identical fixed parameters and deterministic equations; only the initial x coordinate is perturbed. All states, parameters, and time are dimensionless, so the metric/US switch has no numerical effect. The governing equations and constants come from Lorenz (1963). Taylor recurrences, adaptive control, numerical agreement criteria, display projection, and crossing threshold are defined here.

The grey path includes later times than the playback cursor, within the trajectory agreement span. No observations or random forcing enter the model. Identical initial conditions imply exactly zero separation by ODE uniqueness; this remains valid independently of late-time coordinate accuracy.

Where it applies

Educational Lorenz-63 experiment, not a weather forecast or a conversion from model time to days. Crossing depends on the initial state, perturbation, norm, and chosen threshold. Separation can contract and recross; the first crossing does not mean permanent loss of agreement or measure a largest Lyapunov exponent.

Use perturbations of at least 0.001, or exactly zero for a deterministic control. Chaotic amplification also magnifies roundoff and integration errors. A result shown as numerically unresolved must not be interpreted as a crossing time, a no-crossing claim, or a lower bound. Move playback earlier or shorten the window. Agreement between orders estimates reliability but does not prove an absolute error bound. A very brief between-sample threshold excursion can still escape the initial event search.

The default window is 26 model-time units with the cursor at 22. Longer windows up to 40 are available, but may lose numerical agreement. The long-time regression cases retain the reviewer's original 30-unit window and compare against independent SciPy DOP853 and RK45 integrations, including the threshold-46 crossing near 25.2824 and cursor distance near 14.9372.

Zhang et al. (2019), printed pp. 1077 and 1089, supplies separate real-weather context: a roughly two-week estimate for detailed midlatitude weather in those experiments. Some probabilistic, larger-scale, and slower signals retain useful skill for longer. Lorenz-63 does not reproduce the atmosphere's multiscale error cascade.

References

  • Lorenz, E. N. (1963). Deterministic Nonperiodic Flow. Journal of the Atmospheric Sciences 20, 130–141. Equations (25)–(27), section 7 and Table 1. https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2 . Open scan: https://cdanfort.w3.uvm.edu/research/lorenz-1963.pdf . The attached scan has a manually transcribed text layer for only the two cited Table 1 rows; the page images are retained.
  • Zhang, F., et al. (2019). What Is the Predictability Limit of Midlatitude Weather? Journal of the Atmospheric Sciences 76, 1077–1091. https://doi.org/10.1175/JAS-D-18-0269.1 . Open author-hosted copy: https://texmex.mit.edu/ftp/ftp/pub/emanuel/PAPERS/Zhang_etal_2019.pdf
  • SciPy solve_ivp official documentation (DOP853 and RK45, tolerance and event definitions): https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html . Long-time benchmark values were supplied by the independent reviewer and reproduced with both solvers at rtol=2.3e-14; they are numerical verification of the Lorenz equations, not published table entries.

Inputs

InputSymbolUnitDefault
Initial x perturbation δ₀
Run B starts at (x₀ + δ₀, y₀, z₀). Set zero for an identical-state control. For quantitative comparisons, use δ₀ ≥ 0.001.
delta00.001
Initial state x₀
Dimensionless convection amplitude.
x01
Initial state y₀
Dimensionless horizontal temperature-mode amplitude.
y01
Initial state z₀
Dimensionless vertical temperature-mode amplitude.
z01
Separation threshold
First crossing of the full 3D Euclidean state distance. This is a chosen model criterion, not a meteorological skill threshold.
D_limit1
Simulation window (model time)
Dimensionless time. No conversion to hours or days is implied.
T_end26
▶ Playback position (0–1)
Press play to follow both states; drag from 0 (start) to 1 (end). The ODE solutions are cached while this moves.
progress0.8461538461538461

Results

With the default inputs above, the calc gives:

ResultValue
First threshold crossing (model time)18.33
Crossing statusThreshold reached; refinement agrees
3D separation at cursor14.94
Amplification from initial separation14900 ×
Playback time (dimensionless)22
Trajectory refinement agreementThrough model time 26

Checked against sources

Each case below was calculated by this calc and agrees with the value from its source (within the tolerance where one is shown).

Lorenz Table 1 — iteration 5 (t = 0.05)

Source: Lorenz (1963), Table 1, printed p. 136 (PDF p. 7): N=0005, X=0004, Y=0012, Z=0000. Divide printed state integers by 10, using Δτ=0.01. The source truncates tenths and uses its older second-order double-approximation method; tolerance 0.1 allows this coarse precision. No attempt is made to reproduce its numerical discretization. Original DOI: https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2

Inputs: Initial state x₀ = 0; Initial state z₀ = 0; Initial x perturbation δ₀ = 0; Simulation window (model time) = 0.05; ▶ Playback position (0–1) = 1; the others as above.

ResultSourceThis calc
Run A final x0.4 ± 0.10.4291
Run A final y1.2 ± 0.11.256
Run A final z0 ± 0.10.01189

Lorenz Table 1 — iteration 20 (t = 0.20)

Source: Lorenz (1963), Table 1, printed p. 136 (PDF p. 7): N=0020, X=0030, Y=0066, Z=0007. Divide state integers by 10; t=N×0.01. Tolerance 0.1 accommodates the source's truncated tenths and second-order integration. A modern accurate solution is not expected to match every printed digit of that older numerical method. Original DOI: https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2

Inputs: Initial state x₀ = 0; Initial state z₀ = 0; Initial x perturbation δ₀ = 0; Simulation window (model time) = 0.2; ▶ Playback position (0–1) = 1; the others as above.

ResultSourceThis calc
Run A final x3 ± 0.13.055
Run A final y6.6 ± 0.16.638
Run A final z0.7 ± 0.10.7671

Lorenz no-convection equilibrium

Source: Lorenz (1963), section 7, printed p. 136: no convection corresponds to the origin (0,0,0), an exact equilibrium. Independently known limiting solution of the published system; verifies zero-distance and no-crossing handling. https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2

Inputs: Initial state x₀ = 0; Initial state y₀ = 0; Initial state z₀ = 0; Initial x perturbation δ₀ = 0; Simulation window (model time) = 1; ▶ Playback position (0–1) = 1; the others as above.

ResultSourceThis calc
Run A final x00
Run A final y00
Run A final z00
Estimated cursor separation — see numerical status00
Threshold reached (1 = yes)00
Candidate crossing time or simulation end11

Independent long-time benchmark — default threshold

Source: Independent numerical-review benchmark supplied in this task (separation near 14.9372 at t=22), reproduced separately in review-benchmarks.py with SciPy solve_ivp DOP853 and RK45. At rtol=2.3e-14, DOP853 gives D(22)=14.9372040066 and t(D=1)=18.3310908327; RK45 gives D(22)=14.9372029862 and t(D=1)=18.3310908328. This checks the published Lorenz equations using independent integration implementations; these numbers are not claimed to be printed in Lorenz's paper. Solver documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Inputs: Simulation window (model time) = 30; ▶ Playback position (0–1) = 0.7333; the others as above.

ResultSourceThis calc
Threshold reached (1 = yes)11
Candidate crossing time or simulation end18.33118.331
Estimated cursor separation — see numerical status14.93714.937
Crossing result passes numerical agreement (1 = yes)11
cursor_ok11

Independent long-time benchmark — threshold 46

Source: Independent reviewer supplied a converged D=46 crossing near t=25.2824 (not >30). Independently reproduced in review-benchmarks.py using SciPy solve_ivp: DOP853 at rtol=2.3e-14 gives 25.2824443449; RK45 at the same tolerance gives 25.2824440276. Both runs use (1,1,1) and (1.001,1,1), sigma=10, rho=28, beta=8/3. These are numerical benchmarks of the original published equations, not values reported in Lorenz's table. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Inputs: Simulation window (model time) = 30; ▶ Playback position (0–1) = 0.7333; Separation threshold = 46; the others as above.

ResultSourceThis calc
Threshold reached (1 = yes)11
Candidate crossing time or simulation end25.28225.282
Estimated cursor separation — see numerical status14.93714.937
Crossing result passes numerical agreement (1 = yes)11
cursor_ok11

Related Resources

Responses

Your email address will not be published. Required fields are marked *

L ading...