J04: Project 4: a numerical simulator with analytic validation

Prerequisites and observable goals#

Level: Fourth integrated project; numerical correctness rather than animation quality.

Purpose, model, and nonrequirements#

Simulate the nondimensional initial-value problem

Introduce velocity $v=x'$ to obtain the first-order state equation $(x,v)'=(v,-x)$. The exact solution is $(\cos t,-\sin t)$. Energy $E=(x^2+v^2)/2$ is exactly $1/2$ for the continuous model. These facts allow independent validation of both state and invariant behavior.

The tool supports forward Euler and classical fourth-order Runge–Kutta (RK4), a positive integer step count, and a positive finite final time. It reports model, method, steps, step size, numerical final state, exact final state, Euclidean final-state error, and maximum relative energy error over recorded states. It does not model damping, forcing, collisions, real units, adaptive step selection, or safety-critical dynamics.

Architecture and milestones#

Keep the generic vector ODE integrator in A06 separate from the oscillator model and report generator. The integrator accepts a derivative function, initial state, time interval, number of steps, and method. The CLI parses arguments and serializes the report; it does not contain a second hidden implementation of Euler or RK4.

First verify the model by differentiation. Next trace one Euler step by hand. Then test the generic integrator on a known scalar decay problem. Add oscillator state and energy diagnostics. Finally compare step counts $N,2N,4N$ at the same final time for both methods, retaining numerical values rather than only drawing curves.

Acceptance criteria include rejecting zero or negative step counts, invalid time intervals, unknown methods, incompatible derivative dimensions, and nonfinite states under the integrator’s checks. For ordinary tested settings, RK4 errors should decrease substantially under refinement before roundoff dominates. Euler’s oscillator energy behavior should match the algebraic prediction rather than being mistaken for physical energy creation.

Derive the diagnostic before interpreting it#

Euler gives $x_{k+1}=x_k+hv_k$ and $v_{k+1}=v_k-hx_k$. Therefore

The cross terms cancel. Energy grows by factor $1+h^2$ per step for every nonzero step, even though the exact oscillator conserves energy. Over a fixed time with $N=T/h$, this growth tends toward one as $h\to0$, consistent with convergence, but can become unacceptable over many steps at a fixed $h$. Stability and convergence are related but distinct questions.

RK4 evaluates derivatives at four stages and combines them with weights $1,2,2,1$ divided by six. A06 gives the full implementation and derivation context. It has much smaller error on this smooth test at moderate small steps, but is not exactly energy-conserving and is not unconditionally stable for every step or problem.

Complete solutions: J04#

Compare the first step where your work differs. A corrected answer is not yet an independent success.

Open every solution

Complete project implementation#

The full generic Euler/RK4 integrator is printed in A06; this project adds the following complete model and CLI layer.

"""J04: oscillator simulation with analytical and energy validation."""
import argparse
import json
import math
import sys
from .ode_lab import integrate

def run(steps: int = 100, end: float = 1.0, method: str = "rk4") -> dict:
    data = integrate(lambda t, y: (y[1], -y[0]), (1., 0.), 0., end, steps, method)
    x, v = data[-1][1]
    exact = (math.cos(end), -math.sin(end))
    energies = [0.5*(state[0]**2 + state[1]**2) for _, state in data]
    return {"model": "x''=-x; x(0)=1, v(0)=0; dimensionless",
            "method": method, "steps": steps, "h": end/steps,
            "final_state": [x, v], "exact_final_state": list(exact),
            "final_l2_error": math.hypot(x-exact[0], v-exact[1]),
            "max_relative_energy_error": max(abs(e-0.5)/0.5 for e in energies)}

def main(argv=None):
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--steps", type=int, default=100)
    parser.add_argument("--end", type=float, default=1.0)
    parser.add_argument("--method", choices=["euler", "rk4"], default="rk4")
    args = parser.parse_args(argv)
    try:
        print(json.dumps(run(args.steps, args.end, args.method), indent=2, allow_nan=False))
    except (ValueError, ArithmeticError) as error:
        print(str(error), file=sys.stderr)
        return 2
    return 0

if __name__ == "__main__":
    raise SystemExit(main())

Run, compare, and interpret#

python -m solutions.simulation_cli --method euler --steps 100 --end 1
python -m solutions.simulation_cli --method rk4 --steps 100 --end 1
python -m unittest discover -s tests -v

The reference tests compare known solutions and verify selected refinement behavior. To create a convergence report, run $N=10,20,40,80$ at $T=1$ and record the state errors. For a method of order $p$ in its asymptotic regime, the error ratio $e_N/e_{2N}$ tends toward $2^p$: roughly 2 for first-order Euler and 16 for fourth-order RK4. The observed ratios are diagnostics, not replacements for the method’s mathematical error analysis.

If errors stop decreasing, possible causes include roundoff dominance, an inaccurate reference, a state-shape bug, or an unstable step regime. If energy appears constant while phase is wrong, inspect state error separately. If both methods agree but the initial conditions are wrong, agreement still does not validate the intended problem.

Correctness and limitations#

The reporting layer faithfully describes the generic integrator’s final state and the specified analytic reference. The integrator validates the state dimension and finite values at its evaluated stages. Passing these checks prevents some silent representation failures but does not prove every smoothness or stability condition needed by numerical theory. The simulator is deterministic for its fixed inputs in the tested environment, yet last-bit differences can occur on other platforms or implementations.

An analytical validation case is especially valuable because it can expose errors that two numerical implementations share. For example, accidentally returning $(v,+x)$ changes the oscillator into an unstable equation; an animation may still look like motion, but substitution into the original ODE and comparison with $(\cos t,-\sin t)$ immediately fail.

Rubric and unassigned extensions#

Highest reasoning credit derives the exact solution and energy diagnostic. Highest verification credit reports refinement at equal final times and includes invalid inputs. Highest communication credit distinguishes numerical evidence from physical fidelity and explains state units. Damping, forcing, event detection, adaptive methods, and symplectic integrators are unassigned extensions; each requires new validation cases. T05 supplies a separate complete orbit lab, not an automatic extension of this oscillator’s guarantees.

Retention and readiness#

Retrieve A05 error sources and A06 stability. Advance when you can explain a simulator’s error plot in terms of the actual equation and method rather than merely choosing the smoothest-looking curve.