From one numerical step to a method choice you can defend.
Same initial state and step size, three methods. Lesson 1A explains the difference. Figure inputs and rendering.
These notes teach gravitational simulation through experiments you can run in Python: analytic two-body orbits, interacting and restricted three-body systems, planetary N-body dynamics, and softened three-dimensional clusters. Every lesson comes back to three questions:
- What does an integrator assume?
- How do we test what it gives us?
- When does the evidence support using it?
Higher order is a useful starting point, but it does not settle which method to choose. Two second-order methods can cause about 300 times differ in the same planetary test. The difference comes from how diff methods are approximate: the whole orbit, or the smaller interaction between nearly Keplerian orbits. This course unpacks that difference, then tests where it stops helping.
You can read the whole course on GitHub without running anything. Each lesson ties a physical question to the code that answers it, the measured figures, and the numbers behind them.
- New to simulation? Read Setup and your first run, then Chapter 0: States, units, and gravity.
- Want the methods themselves? Open Nine integrators, step by step, or go straight to the code map.
- Choosing a method for your own problem? Use the integrator selection guide. Each suggestion names the conditions it was tested under.
- Checking a result? Open its lesson, read the comparison rules, and follow the evidence links down to the original tables.
- Quick start: install once, run one lesson
- Chapter 0. States, code units, and Newtonian gravity
- Method notes. Euler, RK4, Leapfrog, RK45, Hermite4, BS, GL4, WH2, and WH4
- 1A. Circular orbits: what does one step do?
- 1B. Step size and convergence order
- 1C. Eccentric orbits and adaptive steps
- 1D. Long-term energy, angular momentum, and phase
- 1E. Kepler splitting and the two-body limit
- 1F. From 2D to 3D: rotation consistency
- 2A. A third body: physical change or numerical error?
- 2B. Convergence and efficiency with real perturbations
- 2C. How perturbation strength changes method choice
- 2D. Long-term three-body dynamics
- 2E. Close encounters: resolving a fast event
- 2F. Spatial three-body motion and angular momentum
- 2G. When does a numerical reference stop being useful?
- 2R. The circular restricted three-body problem
- 3A. Five bodies: reuse the solver, widen the diagnostics
- 3B. N and the cost of direct gravity
- 3C. One system, many orbital timescales
- 3D. Self-gravitating clusters without a central star
- 3E. Beyond trajectories: ensembles and uncertainty
- 3F. An evidence-based integrator selection guide
- Part 4 reading path and measured scope
- 4A. Correct the mapping coordinates
- 4B. Resolve the changing motion
- 4C. Pay for an encounter
- 4D. Read the cost of a prediction
- 4E. Follow an evolving spatial system
Run python runpart4.py for 4A–4E, python runpart4.py --only 4E for the secular-cycle study, or python runpart4.py --only 4D to inspect the old regime map without new integrations. The original nine-method registry is unchanged. 4A adds one processed WH2 construction through a dedicated fixed-step loop.
- Which integrator should I try?
- Measured evidence for each regime
- Reading error, cost, and output files
- Find a saved run and follow an experiment
- Code map and design notes
- Glossary
- Mathematical sources and acknowledgments
- Reproducibility and environment, tests, and verification scope
From the folder that contains pyproject.toml, with Python 3.10 or newer:
python -m pip install -e ".[plots,reference]"
python runpart1.py --only 1AThen read Lesson 1A. Figures are saved without opening windows. Add --show if you want to see them as they are made.
Installation, runner options, and troubleshooting
Every lesson and every saved figure can be read on GitHub without Python. To run the experiments, clone or extract the complete folder and open the folder that contains pyproject.toml in your editor. We call this folder the project root.
Use the same Python 3.10-or-newer environment for all four parts. Install once with python -m pip install -e ".[plots,reference]"; the dependency ranges are declared in pyproject.toml. Matplotlib draws the figures; MATLAB is not needed. Part 4 does not require a separate installation or an exact REBOUND version.
Saved runs record the versions that produced them. Their labels describe historical measurements, not required versions for your own run. Compatible versions can change step histories, timings or image pixels. Read your run's reference checks before interpreting small differences. Reproducibility explains the supported ranges and verification reports the combinations actually checked.
In a terminal at the project root:
python -m venv .venvOn Windows PowerShell:
.\.venv\Scripts\Activate.ps1On macOS or Linux:
source .venv/bin/activateThen install the package with the plotting and reference extras:
python -m pip install -e ".[plots,reference]"If PowerShell will not activate the environment, call .\.venv\Scripts\python.exe wherever these notes write python. In VS Code, select the same interpreter. The editable install imports the code from this folder, so your edits take effect without reinstalling.
NumPy handles arrays, Matplotlib draws the plots, and the reference extra installs REBOUND for the IAS15 and external BS references. Most early two-body lessons run without REBOUND. Lesson 1E and a full reproduction need it. The nine course integrators themselves live in the local core and need nothing extra.
python runpart1.py --only 1AThe runner prints the output folder. Open it next to Lesson 1A: you should find three PNG figures, summary and history tables, and saved arrays.
Every numbered experiment writes a fresh timestamped folder, including 1A. Repeating the command keeps the bundled results intact. The lesson links still point to the supplied snapshots; the runner prints the location of your new run.
| Task | Part 1 | Part 2 | Part 3 | Part 4 |
|---|---|---|---|---|
| List lessons | python runpart1.py --list |
python runpart2.py --list |
python runpart3.py --list |
python runpart4.py --list |
| Run one | python runpart1.py --only 1B |
python runpart2.py --only 2B |
python runpart3.py --only 3B |
python runpart4.py --only 4E |
| Run the whole part | python runpart1.py |
python runpart2.py |
python runpart3.py |
python runpart4.py |
| Show figures | Add --show |
Add --show |
Add --show |
Add --show |
To run the four main experiment sequences in order:
python runpart1.py
python runpart2.py
python runpart3.py
python runpart4.pyFull scans and long integrations take longer than the first lesson. GL4 is the slowest, because every step solves a dense Newton system. If you use --show, close each figure window to let the run continue.
Part 1 lessons accept a few parameters through the runner. Select exactly one lesson when you pass them:
python runpart1.py --only 1D --help
python runpart1.py --only 1D --orbits 50Part 1 modules use package-relative imports, so run them through runpart1.py rather than as files. Part 2 and Part 3 scripts can also run directly:
python experiments/three_body/exp2f_spatial_three_body.py --orbits 2 --no-show
python experiments/n_body/exp3e_ensemble_statistics.py --helpFor the complete Part 4 sequence, use the root runner. With no --only argument it runs all five experiments, validates 4E, then builds the lesson analysis views:
python runpart4.pyTo run only the validated final study, use python runpart4.py --only 4E. The labelled B/C/E follow-ups remain a separate command: python experiments/part4/supplement_part4.py --only 4B 4C 4E. Full Part 4 postprocessing reads the supplied follow-up measurements alongside the new primary runs; the batch input record identifies those sources.
The bundle already holds the earlier results that 3F needs. To rebuild the summary from them without integrating again:
python runpart3.py --only 3FInside a full Part 3 batch, 3F reads that batch's 3A–3E plus the existing Part 1 and Part 2 results. Run alone, it searches the bundled result folders for completed runs. You can add other folders explicitly:
python runpart3.py --only 3F --results-root "D:/AstroSimBench-results"Then read the new source_manifest.json. It lists exactly which runs fed the summary. If missing happend, then it stays unmeasured.
| Symptom | What to check |
|---|---|
No module named astrosimbench |
Install from the folder that contains pyproject.toml, with the same interpreter that runs the experiment. |
| The wrong copy of the core is imported | Run python -c "import astrosimbench; print(astrosimbench.__file__)" and check that the path points into this project. |
| Part 1 relative-import error | Use python runpart1.py --only ..., not the experiment file directly. |
| No plot window | Figures are saved by default. Add --show, or open the printed output folder. For 4E, the first figures follow its reference checks; stage progress appears in the terminal and 4E.log. |
| A method row says failed | Read its reason and settings. A numerical failure and an import or programming error mean different things. |
| Orange cells in the regime chart | The scanned settings did not reach that target. Read the target and the reference policy before judging the method. |
The core tests run after installing the test extras:
python -m pip install -e ".[test]"
python -m unittest discover -s tests -vEvery lesson links to the exact functions it uses. If you would rather start from the code, these are the files to open.
| What it does | Where it lives |
|---|---|
| Shared state, force-model, and cost contracts | core.py |
| Direct gravity and analytic jerk | forces.py |
| The integration driver (fixed and adaptive) | integrate() |
| Euler, RK4, RK45 | euler_step(), rk4_step(), rk45_trial() |
| Leapfrog and GL4 | leapfrog_step(), gl4_step() |
| Hermite4 | hermite4_step() |
| Bulirsch–Stoer | bs_trial() |
| WH2, WH4, and Jacobi coordinates | planetary_wh2_step(), planetary_wh4_step(), jacobi.py |
| Processed WH2 (Part 4 addition) | processed.py |
| Part 4 experiments, external adapters and secular diagnostics | Code guide, dynamics.py, zlk.py |
| Part 4 orchestration and optional current-project checks | runpart4.py, check_project.py |
| Kepler propagation | kepler.py, kepler_drift() |
| Restricted three-body model | restricted.py |
| Initial conditions | systems.py |
| Energy, momentum, elements, distributions | diagnostics.py |
| External references and validity windows | references.py |
| The method registry | __init__.py |
The experiments live in experiments/, one main script per lesson. 4E also has an output-sampling validation stage. While editing, python check_project.py checks current links, syntax and command-line entry points. Add --results to recompute saved Part 4 measurements and tables without integrating. These optional checks permit normal source edits and write no report files unless requested. See verification scope and measured limits.
After a lesson, a short example is a good place to try your own question. Each module's docstring states its purpose, assumptions, command, outputs, and limits. The examples call the existing core methods. Their main() and recording functions only organize the run. They add no integrators.
| Example | What it lets you try | Command from the project root |
|---|---|---|
| two_body.py | Connect an initial state, an analytic reference, errors, and costs | python -m examples.two_body |
| nbody.py | Change N and the system geometry through the same interface | python -m examples.nbody --n 5 --method leapfrog |
| efficiency.py | Collect a full cost table and draw a small error/work plot | python -m examples.efficiency |
| reference_window.py | Check when two numerical references stop agreeing | python -m examples.reference_window |
The first two only print diagnostics. The efficiency example overwrites results/efficiency.png and results/efficiency_costs.csv, so copy those first if you want to keep them. The reference example prints a sampled window and needs REBOUND. For a controlled comparison, use the numbered experiments rather than the examples.
| Area | Contents |
|---|---|
| Core | Nine original course integrators plus the Part 4 processed WH2 construction, gravity models, coordinates, diagnostics, and external-reference adapters |
| Experiments | Twenty original lesson scripts plus five Part 4 additions |
| Results | A reading guide to the saved PNG, CSV, JSON, and NPZ outputs |
| Lessons | The reading path, with code and evidence links |
| Documentation | Methods, glossary, sources, design notes, the Part 4 overview, and the evidence pages |
| Examples | Short scripts to adapt after the relevant lessons |
| Tests | Functional checks for the core, external adapters, runners, sampling, output protection and current-project consistency |
| Part 1, Part 2, Part 3, Part 4 runners | Command-line entry points |
Scope These are teaching implementations and finite, measured experiments. IAS15 and the external BS are checked references, not exact trajectories. A run that finishes, a small energy change, and a reached accuracy target are three different outcomes. Every method choice in these notes stays tied to its model, time window, accuracy target, and cost.
This is AstroSimBench v1.0.0. CITATION.cff has the name of author, Jingyi(Seph) Zhang (ORCID), and GitHub reads it for the "Cite this repository" button.
Original code is under MIT; original teaching prose, project-generated figures and result collections are under CC BY 4.0, where rights apply. Licensing and attribution explains the boundaries and third-party dependencies. See the release guide for tagging, Git byte preservation and optional Zenodo archiving.
The chapter-style Markdown navigation is inspired by CSC207-UofT/207-course-notes. The physics lessons and explanations are written for AstroSimBench, and no CSC207 course content is reproduced. See Sources for the mathematical attributions and for what separates a published algorithm from our implementation of it.
I started AstroSimBench expecting a higher-order method to be more accurate and therefore a better choice. But during my experiments, they kept splitting that expectation apart. At one reported setting in Lesson 1D, Leapfrog conserved angular momentum to round-off while its state error was about a hundred times RK4’s. In a matched planetary test in Lesson 2C, two second-order methods differed in state error by a factor of about 300. Order was one variable among several. The project became a way to locate the others and to measure where each one stops helping.
My part was the design, the direction, the review and the interpretation. I chose the questions and the order of the lessons, from a two-body orbit with an exact answer to planetary splitting methods and their correctors. I decided what each experiment holds fixed, what counts as evidence, and which conclusions the measurements support. I enjoyed most was watching one small question open into the logic underneath it, until I could explain why a method behaves the way it does.
I wrote the first versions of the nine course integrators myself, together with much of the code they rest on, including the gravity routines in forces.py, the energy diagnostics in diagnostics.py, the universal-variable Kepler drift in universal_kepler.py, the Jacobi coordinates in jacobi.py, and the Wisdom–Holman steps in wisdom_holman.py, wh2.py and wh4.py. I also wrote early versions of the 1A–1D experiments, though that is a small share of the experiment code. As the experiments grew, many methods were revised to fit them, and that later engineering was AI-assisted. So were most of the experiment code, the plotting and the Part 4 additions: ChatGPT and Claude drafted and revised much of that code and the prose, built the test and verification tooling, and helped check the methods against external references such as REBOUND.
I compared implementations with their published equations, investigated changes in the numerical results, and rejected conclusions the data did not support. Published algorithms and coefficients are credited to their sources. The scripts, saved measurements and verification records are included so that readers can repeat the checks and examine the reasoning. I began by asking which method was best. I now ask what the structure of a problem allows a method to leave out.
If you find a mistake, tell us the lesson, the exact command or paragraph, what you expected, and the figure or data that shows the problem. A proposal for a new method should name the physical model and the measurable objective it targets.