Lecture 02: Initial-Value Problems and Error Control — Lecture Notes

Terms Developed in Context

The terms below connect the mathematical model, implementation choices, computed diagnostics, and limits of the scientific claim.

Initial-value problem

A differential equation paired with a state specified at a starting time. In this lecture it is used to address How do local error estimates control accumulated error over an astrophysical timescale?

Model connection: Read the term against Local truncation error and the governing relation \(\dot{\mathbf y}=\mathbf f(t,\mathbf y),\quad e_{\rm local}\leq atol+rtol\,|\mathbf y|\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: A tiny local error estimate cannot cure unstable equations, unresolved events, poor scaling, or a wrong right-hand side. This distinction separates numerical behavior from physical accuracy and model adequacy.

Local truncation error

The one-step defect made by a numerical integrator assuming the incoming state is exact. In this lecture it is used to address How do local error estimates control accumulated error over an astrophysical timescale?

Model connection: Read the term against Global error and the governing relation \(\dot{\mathbf y}=\mathbf f(t,\mathbf y),\quad e_{\rm local}\leq atol+rtol\,|\mathbf y|\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: A tiny local error estimate cannot cure unstable equations, unresolved events, poor scaling, or a wrong right-hand side. This distinction separates numerical behavior from physical accuracy and model adequacy.

Global error

The accumulated difference between the numerical and exact trajectories at a later time. In this lecture it is used to address How do local error estimates control accumulated error over an astrophysical timescale?

Model connection: Read the term against Adaptive timestep and the governing relation \(\dot{\mathbf y}=\mathbf f(t,\mathbf y),\quad e_{\rm local}\leq atol+rtol\,|\mathbf y|\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: A tiny local error estimate cannot cure unstable equations, unresolved events, poor scaling, or a wrong right-hand side. This distinction separates numerical behavior from physical accuracy and model adequacy.

Adaptive timestep

A step size changed by an error estimator to meet relative and absolute tolerance targets. In this lecture it is used to address How do local error estimates control accumulated error over an astrophysical timescale?

Model connection: Read the term against Stiffness and the governing relation \(\dot{\mathbf y}=\mathbf f(t,\mathbf y),\quad e_{\rm local}\leq atol+rtol\,|\mathbf y|\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: A tiny local error estimate cannot cure unstable equations, unresolved events, poor scaling, or a wrong right-hand side. This distinction separates numerical behavior from physical accuracy and model adequacy.

Stiffness

A separation of dynamical timescales that makes explicit integration inefficient or unstable. In this lecture it is used to address How do local error estimates control accumulated error over an astrophysical timescale?

Model connection: Read the term against Initial-value problem and the governing relation \(\dot{\mathbf y}=\mathbf f(t,\mathbf y),\quad e_{\rm local}\leq atol+rtol\,|\mathbf y|\). The relevant computational evidence is a controlled output, residual, invariant, or refinement sequence rather than visual smoothness alone.

Interpretive limit: A tiny local error estimate cannot cure unstable equations, unresolved events, poor scaling, or a wrong right-hand side. This distinction separates numerical behavior from physical accuracy and model adequacy.

Problem and physical meaning

An ODE solver advances a state by approximating a local trajectory. Adaptive Runge-Kutta methods estimate error from embedded formulas, but tolerances constrain an error estimator rather than guaranteeing a global physical error bound.

How do local error estimates control accumulated error over an astrophysical timescale? Begin from the scientific observable or state variable, and name what is directly specified versus inferred. The model should be simple enough to validate, but not quietly presented as the entire astrophysical system.

Model and derivation

The governing expression for this lecture is:

Define the symbols and the independent variable before numerical evaluation. Dimensional analysis checks that both sides have matching units; if the model is dimensionless, state the adopted scale factors and the transformation back to physical units.

An ODE solver advances a state by approximating a local trajectory. Adaptive Runge-Kutta methods estimate error from embedded formulas, but tolerances constrain an error estimator rather than guaranteeing a global physical error bound.

Numerical implementation

Choose a solver from stiffness and conservation structure; scale variables and tolerances; compare at a fixed output grid; halve the maximum step or tighten tolerances; compare with an analytic solution.

Keep physical quantities and numerical arrays conceptually separate. Convert units at well-defined boundaries, preserve raw inputs, and record dtype, solver, tolerances/timestep, grid, and software versions. A numerical method is part of the model evidence, not an invisible implementation detail.

Worked, programmatically computed example

At rtol=1e-03, endpoint error is 0.0002106; at rtol=1e-12, it is 7.577e-14. The analytic endpoint is 0.36787944.

These values are generated in src/astr474_computations.py, then shared with lecture/lab/assignment renderers. Re-run the generator after changing an input; do not hand-edit a displayed result without recomputing dependent artifacts.

At least one independent check is required: an analytic solution, a second numerical method, a conservation law, or a controlled convergence sequence. When these disagree, diagnose units, frame, discretization, source terms, and model scope before selecting a result.

Error budget and limitations

A tiny local error estimate cannot cure unstable equations, unresolved events, poor scaling, or a wrong right-hand side.

Classify each uncertainty before combining terms: local/global solver error, timestep or spatial discretization, Monte Carlo sampling, parameter uncertainty, measurement uncertainty, and model discrepancy are not interchangeable. Conservation can test internal consistency but does not prove the model is physically complete.

Connection to evidence and research practice

Lab 01 compares SciPy RK45 outputs with the analytic exponential-decay solution across tolerances.

For real parameters, the source record distinguishes the fetched product and epoch from derived values. For synthetic experiments, the initial conditions, generator, seed, and exact file are recorded. A simulation is not an observation, even when its parameters resemble a real object.

Misconceptions to avoid

• More displayed digits do not imply more physical accuracy.
• A successful solver exit or smooth plot does not establish convergence.
• Agreement between two solvers with the same wrong model is not independent physical validation.
• A formal covariance or Monte Carlo interval is conditional on the adopted model and noise process.

Study and review questions

1. How do local error estimates control accumulated error over an astrophysical timescale?
2. Which numerical parameter would you vary first, what observable would you track, and how would you know the sequence reached an asymptotic regime?
3. Design one independent check and state what failure would change your scientific interpretation.

Preparation for the next lecture: explain how this method's assumptions change when the problem becomes more coupled, stochastic, or data-driven.

Lecture-specific derivation and assumptions

An embedded Runge-Kutta pair advances y_n with a high-order estimate and forms a lower-order estimate from the same stage evaluations. Their difference estimates local truncation error. A controller chooses h_new roughly as h*(tol/error)^(1/(q+1)), clipped to avoid extreme changes; this controls a local weighted norm, not the global error directly. For y'=-y, separation gives y(t)=y(0)e^-t, so at t=1 the exact reference is e^-1. Global error reflects accumulated local defects and stability of the propagation.

Algorithm trace and computed-output interpretation

For each tolerance, record solver success, accepted output points, endpoint value, and |y(1)-e^-1|. Plot actual endpoint error against requested rtol on log axes. Tightening rtol should initially reduce error, then plateau when roundoff, method behavior, or absolute-tolerance scaling dominates. Repeat with a different method and fixed physical output time; compare errors rather than the number of returned samples alone.

Computed result: At rtol=1e-03, endpoint error is 0.0002106; at rtol=1e-12, it is 7.577e-14. The analytic endpoint is 0.36787944.

A solver's tolerance is an internal control target, not a guarantee that the physical answer is accurate to that many digits. Scaling matters when state components have different units or magnitudes; use componentwise absolute tolerances or nondimensionalize. Event timing, stiffness, and discontinuous forcing require additional checks not supplied by this smooth decay benchmark.

Two-meeting teaching arc

Each weekly module is taught in two 75-minute meetings. The sequence below includes a formative check in each meeting and is intended to use the lecture's worked computation, not merely present it.

1. Meeting 1 (75 min): 10 min IVP examples; 20 min local truncation derivation; 25 min embedded RK stages and controller; 15 min analytic decay solution; 5 min exit prediction.
2. Meeting 2 (75 min): 10 min recap; 20 min tolerance scaling and dimensional state variables; 25 min execute/interpret the convergence sweep; 15 min stiffness/events and method choice; 5 min write one falsifiable solver check.

Targeted study questions

1. Why can reducing rtol stop improving the endpoint error?
2. What changes in the error norm if the state has both meters and kilograms?
3. Which independent test would expose a sign error in the right-hand side?