Rocket Flight Capstone

The four lessons in this module steered through a complete journey from a physical model to a numerical result that can support an engineering judgment. Writing code was one necessary but not sufficient component. You reach a trustworthy computation over a thoughtful build-up of technical steps.

So far you have:

• derived a model based on ordinary differential equations,

• reconstructed time-stepping methods, and explored their order of accuracy

• located events that fall between stored times,

• studied convergence with numerical experiments, and

• compared methods at a required accuracy.

You also practiced specifying tasks for an AI agent and auditing the code it proposed. This is a key new technical skill in the era of agentic AI, and like any other, it requires practice. That is the purpose of this exercise.

This capstone brings all these ideas together in a new setting: the vertical flight of a small rocket. The physical model involves gravity, aerodynamic drag, changing propellant mass, and known transitions from powered flight to coasting flight, reaching apogee, and falling back to hit the ground.

You will first do a fixed-step RK2 calculation and then compare it with an adaptive solver from the SciPy library. Neither result is authoritative simply because the code runs or comes from a standard library, or produces correct-looking results. Your task as always is to decide what evidence makes the reported flight quantities trustworthy.

Engineering brief¶

A simulation team needs reference results for an idealized vertical-flight model. These results will be used to check later implementations of the same mathematical model. The quantities of interest are the rocket’s maximum speed, apogee, and impact conditions.

The engineering question is:

Can we make accurate predictions of the rocket’s apogee and impact conditions to the required numerical accuracy, using a fixed-step RK2 calculation and/or a properly configured adaptive SciPy solver?

These are the accuracy targets:

• maximum speed to within $0.1\ \mathrm{m/s}$;

• apogee altitude to within $1\ \mathrm{m}$;

• apogee time to within $0.02\ \mathrm{s}$;

• impact time to within $0.05\ \mathrm{s}$; and

• impact speed to within $0.1\ \mathrm{m/s}$.

The targets apply to numerical error within the stated model. They do not claim that this simplified model predicts a real rocket to the same accuracy.

Your final verdict will identify the evidence supporting an RK2 step size and a SciPy tolerance configuration, explain how burnout and the flight events were handled, and separate numerical uncertainty from limitations of the model. Code generation may be assisted, but you remain responsible for recording the specification of delegated work, auditing the result, and deciding what the evidence supports.

The vertical-flight problem¶

We model the rocket as a point mass constrained to move vertically. Altitude $h$ is measured upward from the launch point, and velocity $v=dh/dt$ is positive during ascent and negative during descent. The rocket launches from rest at ground level with $100\ \mathrm{kg}$ of propellant. Here are all the problem settings:

The idealized model assumes constant gravity, air density, drag coefficient, reference area, exhaust speed, and burn rate during powered flight. It neglects wind, lateral motion, atmospheric variation, rocket attitude, and any change in aerodynamic properties. The engine switches off when the propellant is exhausted.

Propellant burn rate¶

A positive burn rate $\mu(t)$ removes propellant mass from the rocket, thus $dm_p/dt=-\mu(t)$. With the stated initial mass and constant powered-flight burn rate, burnout occurs at

The prescribed mass-flow history is a step function:

The value changes discontinuously at the known burnout time, $t_b$. You should treat burnout as an explicit boundary between two integrations, rather than allowing one numerical step to cross it unnoticed.

Let’s begin by loading the numerical Python libraries, and setting up the problem. Then make a plot of the burn rate.

Integrating Equation 2 gives the exact remaining propellant mass:

This exact history will provide a simple but important check of the numerical state.

Derivation of the equations of motion¶

Mass balance¶

The rocket’s instantaneous mass is the sum of its constant dry mass and remaining propellant mass:

During a short interval $dt$, a positive mass $dm_e=\mu(t)dt$ leaves the rocket. Consequently,

Using a separate symbol $\mu$ for the positive outflow rate avoids a common sign ambiguity: $\mu$ is positive while propellant is burning, whereas $dm_p/dt$ is negative.

Newton’s second law and thrust¶

The engine expels propellant downward at the positive rate $\mu(t)$ and effective exhaust speed $v_e$ relative to the rocket. In this model, the corresponding upward thrust is

Once the exhaust momentum flux is represented by the thrust, Newton’s second law for the remaining rocket is simply

where $F_{\mathrm{ext}}$ is the sum of gravity and aerodynamic drag. This is an effective exhaust-speed model and other engine-flow details are ignored.

Gravity and aerodynamic drag¶

Gravity acts downward with force $-mg$. Aerodynamic drag has magnitude $\tfrac12\rho AC_Dv^2$ and always opposes the motion. In one signed velocity component, both directions are represented by

When $v>0$, the drag force is negative and opposes ascent. When $v<0$, $v|v|<0$, so the drag force is positive and opposes descent. Replacing $v|v|$ by $v^2$ would silently give the wrong direction during descent.

The total external force is therefore

The initial-value problem¶

Substituting Equation 6 and Equation 11 into Equation 7, using $m=m_s+m_p$, and adding the altitude and propellant equations gives

The three initial conditions are

In vector form, with $u=[h,v,m_p]^T$, the model has the non-autonomous form $u'=f(t,u)$. The explicit appearance of time through $\mu(t)$ means that the RK2 step used later must evaluate its second stage at both the midpoint state and the midpoint time.

Flight regimes and reported events¶

The solution passes through three physical regimes and three transitions:

Burnout is a breakpoint known in advance, while apogee and impact are events discovered from the computed solution. Impact must mean the downward ground crossing after the rocket has been aloft; the initial condition $h(0)=0$ is the launch point, not an immediate impact.

Predictions and checks before computing¶

A numerical trajectory should be judged against expectations that do not come from the same computation. Three useful checks follow directly from the model.

• Exact propellant history. Use Equation 3 to calculate the remaining propellant at $t=3.2\ \mathrm{s}$ and confirm that it reaches exactly zero at burnout.

• No-drag burnout-speed bound. If drag is removed from the model during powered flight, integrating the velocity equation gives

During ascent, drag can only reduce the velocity from this ideal value. Evaluate Equation 14 at burnout and use it as an upper bound for the powered-flight speed.

• Empty-shell terminal speed. During descent after burnout, the mass is $m_s$ and a steady downward velocity satisfies $dv/dt=0$. Its magnitude is

The rocket begins its descent from rest at apogee, so its impact-speed magnitude should approach this value from below.

Finally, every term in Equation 7 has units of force: $mg$, $\mu v_e$, and $\rho Av^2$ each have units $\mathrm{kg\,m/s^2}$. Take notice of dimensional agreement and any sign conventions.

Build a time-aware RK2 flight solver¶

The numerical calculation needs three distinct layers:

Keeping these responsibilities separate makes each part easier to inspect and test.

The state equation $u'=f(t,u)$ depends explicitly on time because the mass-flow rate changes at burnout. The RK2 step will need to adopt this time dependence explicitly.

Compute the model’s derivatives¶

The right-hand-side function returns the derivatives in the same order as the state-vector components $u=[h,v,m_p]^T$. Naming the three forces separately makes their signs inspectable before they are combined into the acceleration. The mass-flow function is passed as an argument so that the time dependence is an explicit dependency of the model.

Advance one time-aware midpoint step¶

Lesson 4 introduced explicit-midpoint RK2 for an autonomous system (meaning, time does not appear explicitly on the right-hand side). For $u'=f(t,u)$, both the midpoint state and midpoint time must be used. Rewrite the RK2 algorithm in this textbook form:

The function below is still a single numerical step. It knows nothing about burnout, apogee, or impact. In particular, shortening a step to land exactly at burnout belongs in the driver; it is not part of the RK2 formula.

Interpolate a bracketed event¶

A step brackets an event when the monitored component is positive at the beginning and nonpositive at the end. The same linear interpolation used for touchdown in Lesson 4 estimates the event time and the complete state. This helper function assumes that the driver has already checked the bracket.

Drive the calculation through burnout and flight events¶

The driver normally advances by the requested fixed step $\Delta t$. It shortens a step only when necessary to end exactly at the known burnout time or at the time limit. This is breakpoint handling, not adaptive time stepping (whose goal is error control).

After each completed step, the driver looks for the first downward crossing of $v=0$ and then the first downward crossing of $h=0$. It stores interpolated apogee and impact states and returns a status rather than assuming that impact must have occurred.

Make one provisional flight¶

Use $\Delta t=0.6\ \mathrm{s}$ for a first diagnostic run. This value deliberately does not divide the $5\ \mathrm{s}$ burn time, so the result exposes whether the driver lands on the breakpoint as intended. It is a provisional step, not a claim that the engineering accuracy targets have been met. A $60\ \mathrm{s}$ time limit comfortably allows the anticipated trajectory to reach the ground.

Before running the cell, predict the result status and the order of burnout, maximum speed, apogee, and impact.

Plot the three state components against time. Use the figure as a diagnostic: look for the expected linear decrease of propellant during the burn, continuity of altitude and velocity at burnout, one apogee, and a return to $h=0$. A smooth-looking curve cannot establish numerical accuracy.

Decide what the first run supports¶

Develop preliminary confidence by combining evidence that would fail in different ways:

1. Code-to-equation inspection: point from each returned derivative to the matching term in Equation 12, including its sign.

2. Exact mass check: compare the reported propellant at $3.2\ \mathrm{s}$ and at burnout with your paper values. Confirm that burnout occurs at exactly $5\ \mathrm{s}$ even though $0.6\ \mathrm{s}$ does not divide it.

3. Event logic: confirm that the status is impact, the event order is physically sensible, the apogee state has $v=0$, and the impact state has $h=0$.

4. Independent bounds: compare the maximum speed with the no-drag burnout-speed bound and the impact-speed magnitude with the empty-shell terminal speed from your paper record.

5. Trajectory inspection: use the plot to look for wrong signs, jumps, repeated events, or propellant used after burnout.

Record any discrepancy before changing code. Passing these checks means that the calculation has survived useful attempts to expose an obvious defect. It does not show that $\Delta t=0.6\ \mathrm{s}$ meets the acceptance targets, and agreement of printed digits from one run is not an accuracy argument. Keep this result unchanged as the baseline for the verification and refinement work that follows later.

Verify and refine RK2 with Richardson extrapolation¶

The provisional trajectory was physically plausible and passed several exact and qualitative checks. We now ask a narrower numerical question: how much discretization error remains in each reported quantity, and what time step is supported by the evidence?

Following the terminology introduced by Roache (1998) (and now standard), this is solution verification: estimating how accurately the differential equations have been solved for a particular problem setting. It is not validation of the rocket model against reality. Refining one RK2 implementation also cannot expose every shared coding or modeling error, which is why the independent checks already performed—and the different solver introduced later—remain necessary.

Estimate discretization error from a grid family¶

Let $Q$ be one quantity of interest, such as apogee altitude, and let $Q_{\Delta t}$ be its value computed with time step $\Delta t$. In the asymptotic refinement regime, a method of order $p$ has an error model of the form

where $Q$ is the exact value of the quantity for the stated continuous model and $C$ does not depend on $\Delta t$. For a refinement ratio $r$, compare a coarse result $Q_{\Delta t}$ with a finer result $Q_{\Delta t/r}$. Eliminating the leading error term gives the Richardson-extrapolated value

The magnitude of the correction is an estimate of the discretization error in the finer result:

For midpoint RK2 with step halving, $p=2$, $r=2$, and the denominator is $2^2-1=3$. The difference between two runs is therefore not itself the estimated error in the finer result; the leading-order estimate is one third of that difference.

The assumed order should also be tested. Three successive grids give the observed order

Values approaching 2 support the RK2 error model. Do not expect exact agreement on every triplet: coarse runs may not yet be asymptotic, while event interpolation and very small successive differences can make an observed-order estimate irregular. The absolute value makes the logarithm usable, but a reversal in the signs of successive differences is still a warning to investigate. Inspect a sequence rather than selecting one convenient value.

Plan the refinement study before generating code¶

Use the halving family

Because every integration stops explicitly at $t_b=5\ \mathrm{s}$, the discontinuity is a common boundary for every grid. The refinement study will examine these targets:

The maximum-speed target matches the impact-speed target and corresponds to about $0.04\%$ of the anticipated maximum speed. It is tight enough to discriminate a coarse result while remaining proportionate to the other requirements.

Record the agent task brief¶

You can use the brief below, copied into a Markdown cell in your working notebook, but feel free to improve it as you see fit. Above it, record the date and the agent or persona used.

Before attaching your working notebook, confirm that it contains the NumPy and Matplotlib imports, model parameters and mass-flow function, all four solver function definitions, and the setup cell defining u_0, time_limit, and rhs_args. Include these cells even if you have not copied the surrounding lesson narrative. The agent needs their definitions in the saved notebook; values that exist only in your running kernel are not available in the attachment.

Agent RK2-refinement task brief¶

Goal and scope— Your goal is to provide the driver code to perform solution verification of a vertical rocket-flight numerical solver. Add two readable, unexecuted Python cells immediately below this brief. Cell 1 runs the RK2 rocket calculation over the prescribed time-step family. Cell 2 computes and plots the Richardson quantities. Generate code that will report numerical evidence when I run it; leave the supported time-step choice and all conclusions to me.

Non-negotiables– Before adding code, inspect the attached notebook for the NumPy (np) and Matplotlib (plt) imports; mass_flow_rate(), rocket_rhs(), rk2_step(), interpolate_zero_crossing(), and integrate_rocket_rk2(); and the model parameters and setup values u_0, time_limit, t_burn, rhs_args, m_p0, and mu_0. Confirm that the driver accepts integrate_rocket_rk2(u_0, dt, time_limit, t_burn, *rhs_args) and returns the fields needed below: status, times, states, steps, rhs_evaluations, maximum_speed, maximum_speed_time, apogee_time, apogee_state, impact_time, and impact_state. State components are ordered as altitude, velocity, and propellant mass; states has one row per stored time. If a required definition is missing or the interface differs, stop before adding cells and identify what I need to supply or clarify. Do not invent or reconstruct missing dependencies.

Reuse the existing functions and setup values exactly as defined. Do not redefine or modify the model, RK2 step, event interpolation, or flight driver.

• Cell 1 must use dt_values = [0.6 / 2**level for level in range(7)], in that order. Run one flight for each value and retain one dictionary per run in a list named rk2_refinement_runs.

• Each run record must contain the time step, status, completed steps, RHS evaluations, maximum speed and its time, apogee altitude and time, impact time, positive impact-speed magnitude, and maximum absolute propellant error relative to $m_p(t)=\max(m_{p,0}-\mu_0t,0)$ at the stored times. If a run does not end with status == 'impact', retain its status and work but use None for unavailable flight quantities.

• Print one readable row per run. Do not hide failed or unfinished runs.

• Cell 2 must define these five quantities and targets: maximum speed, 0.1 m/s; apogee altitude, 1.0 m; apogee time, 0.02 s; impact time, 0.05 s; and positive impact-speed magnitude, 0.1 m/s.

• A pair or triplet is usable only when every run has status == 'impact' and the runs occupy consecutive entries in the original dt_values family. Do not filter out unsuccessful runs and then compare across the resulting gaps.

• For each usable adjacent coarse–fine pair, let Q_coarse be the result at step 2*h and Q_fine the result at step h. With $p=2$ and $r=2$, calculate Q_ext = Q_fine + (Q_fine - Q_coarse) / 3, estimated_error = abs(Q_fine - Q_coarse) / 3, and normalized_error = estimated_error / target. Preserve the signed correction in Q_ext; the error estimate applies to the fine result.

• For every usable three-grid sequence, calculate the observed order as $p_{\mathrm{obs}}=\log_2\left|\left(Q_{\Delta t}-Q_{\Delta t/2}\right)/\left(Q_{\Delta t/2}-Q_{\Delta t/4}\right)\right|$, using three consecutive results in coarse-to-fine order. If either required difference is zero, record None instead of dividing by zero. Attach the observed order from steps (4*h, 2*h, h) to the Richardson record for the pair (2*h, h). Record None when the required usable triplet is absent, including for the first pair.

• Retain one dictionary per quantity and grid comparison in a list named richardson_records. Print the fine time step, fine result, extrapolated result, estimated error, normalized error, and observed order.

• Plot normalized Richardson error against the fine time step for all five quantities on logarithmic axes. Include a horizontal line at 1 for the acceptance target. Use the natural numerical ordering of the time-step axis: fine steps on the left and coarse steps on the right.

• Use ordinary loops, dictionaries, NumPy, and Matplotlib. Keep units visible in printed headings or labels. Do not select a time step, label a run accepted, or write an engineering conclusion.

Allowed actions– Read the attached notebook and add only the two requested cells. Do not run code, change existing cells, create or modify other files, access the network, install packages, import scipy or pandas, or replace the requested calculations with a new helper library.

Done when– Both cells are present without outputs or execution counts. Stop and report where they were inserted, what variables they create, and any ambiguity you encountered. Do not report numerical results.

Invoke the agent¶

Save your working notebook, attach it to a new agent conversation, and send only this request:

An agent may produce correct-looking formulas with the wrong grid orientation, compare the wrong rows, or silently redefine a function. Its completion report is a claim to inspect, not evidence that the study is correct.

Audit before running¶

Inspect the proposed experiment. Confirm that earlier cells are unchanged and the two new cells have no outputs. Then inspect the code against the brief:

• Are the seven time steps generated by repeated halving, with the same initial state, model parameters, event logic, and time limit in every run?

• Does every run remain visible, including a non-success status?

• Are impact speed and both error measures treated as positive magnitudes, while the Richardson extrapolation retains the signed fine–coarse correction?

• Is the fine-grid error estimate divided by $2^2-1=3$?

• Does each observed-order calculation use three consecutive grids in the correct coarse-to-fine order?

• Does the propellant comparison use the exact history at the stored times?

• Does the plot divide each estimated error by the matching target, without comparing quantities that have different units directly?

• Is there any new model, solver, import, file operation, execution command, hidden conclusion, or code you cannot explain?

If a line of code is too compact or dense to audit, ask the agent to rewrite only that cell with ordinary loops and intermediate names. Read the revision before running it. Do not use Run All.

Run the study and judge the evidence¶

Run Cell 1 first. Confirm that the $0.6\ \mathrm{s}$ row reproduces your saved provisional result, every status is impact, the work grows as expected under halving, and the exact propellant errors remain near floating-point roundoff. Investigate any discrepancy before continuing.

Run Cell 2. Choose one quantity and one three-grid sequence, then reproduce its two differences, observed order, extrapolated value, and error estimate with a calculator. Check these against the generated table before interpreting the whole study.

Look for the region where successive changes decrease and the observed orders tend toward 2. A candidate RK2 step is supported when all five estimated errors are below their targets and the next finer run gives the same outcome with smaller estimates. Prefer the coarsest step supported by that evidence; a smaller step costs more without improving the stated decision.

Do not conceal an irregular observed order. If an output is already changing by much less than its target but its observed order is unsettled, record both facts and retain the next finer run. Event-derived quantities can show a less regular sequence because the event falls at a different fraction of each step. The later comparison with another solver will provide complementary evidence.

In your notebook, record:

• the RK2 step you provisionally support;

• the controlling quantity and its Richardson estimate;

• the observed-order evidence;

• the supporting finer run;

• the total RHS work for the selected run; and

• any irregularity or limitation that remains.

This is the end of the RK2 verification and refinement stage. Do not yet treat agreement within this one implementation as an independent confirmation of the answer.

Bridge to adaptive ODE integration in SciPy¶

Our RK2 solver uses a step size chosen before the run and repeats the calculation on finer grids to estimate its error. Scientific libraries also provide adaptive solvers that estimate local error while integrating and change their step size accordingly.

SciPy’s solve_ivp() documentation describes a common interface to several ODE methods. We will use RK45, an explicit Runge–Kutta method suitable as a starting point for this nonstiff model. A library implementation reduces the code we must write, but it does not choose meaningful tolerances, handle known discontinuities automatically, or decide whether the reported quantities are trustworthy.

The solve_ivp() contract¶

The solver call has this shape:

from scipy.integrate import solve_ivp

result = solve_ivp(
    fun=right_hand_side,
    t_span=(t_start, t_end),
    y0=initial_state,
    method='RK45',
    rtol=relative_tolerance,
    atol=absolute_tolerances,
    events=event_functions,
    dense_output=True,
    max_step=maximum_step,
)

The right-hand side model function must have the signature fun(t, y) and return derivatives with the same shape as y. The pair t_span=(t_start, t_end) defines the integration interval, and y0 is the complete state vector at its beginning. Extra model data can be supplied with args or captured in a small wrapper function.

SciPy’s ODE solvers can track events, which are the zeros of a continuous function of time and the state vector. The argument events expects a function name, or list of functions, with signature event(t, y) and additional arguments passed via args. The solver looks for the root of event(t, y(t)) = 0 in each time step, and returns the values of t and y corresponding to the event taking place.

What RK45 adapts¶

RK45 uses an embedded pair: formulas of two different orders estimate the local error, and the solver accepts, rejects, and resizes steps to satisfy a tolerance test. SciPy compares the estimated error in each state component with a scale of the form

The relative tolerance rtol is dimensionless. Because altitude, velocity, and propellant mass have different units and scales, a component-wise array such as atol=np.array([atol_h, atol_v, atol_mp]) is more interpretable than one absolute tolerance for all three components.

These are controls on estimated local state error, not direct guarantees on apogee or impact error. We must still tighten the tolerances and examine changes in the quantities of interest.

Two options are easy to confuse. max_step places a ceiling on an internal step; it does not make the method fixed-step. t_eval requests times at which results are stored; it does not prescribe the solver’s internal steps. With dense_output=True, result.sol(t) evaluates the solver’s continuous interpolant at other times.

Events are roots; burnout is a boundary¶

A SciPy event function should return a scalar that defines the event when it becomes zero. These functions also have optional attributes that specify whether integration stops and which crossing direction counts. For example, we might assign:

def apogee_event(t, y):
    return y[1]  # velocity

apogee_event.direction = -1   # detect positive-to-negative crossings
apogee_event.terminal = False # record the event and continue

def impact_event(t, y):
    return y[0]  # altitude

impact_event.direction = -1
impact_event.terminal = True

The assignments after each function definition do not call the function. In Python, a function is an object, so it can carry attributes accessed with dot notation. solve_ivp() reads the direction and terminal attributes as configuration metadata when it receives the function through events.

The negative direction selects downward zero crossings. It distinguishes apogee from any upward crossing of $v=0$ and impact from the upward departure at $h(0)=0$. SciPy locates a detected root within an accepted step and returns event times and states in t_events and y_events. Because detection begins with a sign change over a step, multiple crossings within one step can be missed; max_step can provide an additional safeguard when the event behavior demands it.

The attribute terminal tells the solver whether to stop at the event. We record apogee and continue, but stop at impact.

Burnout is different: its time is known, and the right-hand side changes discontinuously there. The SciPy calculation should therefore use two calls—powered flight ending exactly at $t_b$, followed by unpowered flight beginning from that state. Adaptive stepping is not a reason to integrate blindly across a known discontinuity.

Inspect the returned result¶

The value returned by solve_ivp() is an object that contains both solution data and diagnostic information in fields as follows:

Notice that result.y stores state components by rows, whereas our RK2 driver’s states array stores them by columns. A successful status means that the numerical algorithm reached its requested endpoint or terminal event; it does not establish that the tolerances were adequate or that the physical model is valid.

The next stage will ask an agent to propose the two-stage (burn, coast/fall) SciPy implementation. Your task will be to audit it against this interface before running or trusting it.

Agent-proposed SciPy implementation¶

You now know enough about solve_ivp() to inspect its interface, but not enough to be expected to design a complete two-stage implementation from scratch. An agent will propose that code under a detailed numerical specification. You will control its editing authority, inspect the proposal before execution, test the important interfaces, and decide what the results support.

First revisit Lesson 4’s agent task brief and pre-run audit. You will emulate the parts of that workflow that are already familiar:

• record the date, agent or persona, and working-notebook filename;

• authorize edits only in the requested location;

• require unexecuted cells and a change record;

• inspect every proposed line before running it; and

• restrict any correction to the cell that needs it.

The SciPy-specific requirements are supplied below. You are not being asked to invent an unfamiliar library interface.

Preserve the accepted RK2 comparison¶

In your working notebook, add the following cell immediately before the agent brief. Replace the blank with the RK2 step supported by your Richardson study. This creates a fresh result with the standard driver interface, independent of any key names chosen by the earlier agent-generated table.

dt_rk2_accepted = _____
rk2_accepted = integrate_rocket_rk2(
    u_0, dt_rk2_accepted, time_limit, t_burn, *rhs_args
)
assert rk2_accepted['status'] == 'impact'

Run this cell and confirm that it reproduces the row you selected. The SciPy agent may read rk2_accepted but may not modify it.

Record the agent task brief¶

Copy the brief below into a Markdown cell in your working notebook. Fill in the notebook filename and place your provenance note above it, following the Lesson 4 pattern.

Agent SciPy-implementation task brief¶

Goal and scope

Add three readable, unexecuted Python cells immediately below this brief in [notebook filename]. Cell 1 defines a two-stage SciPy RK45 flight integrator. Cell 2 runs a prescribed tolerance study. Cell 3 compares successive SciPy results, the accepted RK2 result, and one maximum-step sensitivity run. Report evidence only; do not select a tolerance configuration or write a verdict.

Non-negotiables

Do not change any existing cell or redefine rocket_rhs(), rk2_step(), interpolate_zero_crossing(), or integrate_rocket_rk2().

Cell 1 — SciPy functions

• Import only solve_ivp from scipy.integrate. Define readable powered_rhs(t, y) and coast_rhs(t, y) wrappers that call the existing rocket_rhs() with the existing parameters. The powered wrapper must use the constant rate mu_0 throughout the powered phase, including its endpoint; the coast wrapper must use zero mass flow.

• Define apogee_event(t, y) returning velocity and impact_event(t, y) returning altitude. Both have direction = -1; apogee is nonterminal and impact is terminal.

• Define integrate_rocket_scipy(u_0, rtol, atol, max_step, time_limit). It must call solve_ivp() with method='RK45', dense_output=True, no t_eval, and the supplied rtol, component-wise atol, and max_step.

• Use two solver calls. Powered flight integrates from 0.0 to t_burn without event functions. Check its success before continuing. Save its final propellant value as burnout_mass_residual, then copy the burnout state and set its propellant component exactly to zero from the known mass balance. Unpowered flight starts from that state at t_burn and integrates to time_limit with events ordered as (apogee_event, impact_event).

• Do not integrate across burnout with the piecewise mass_flow_rate() function. Do not locate SciPy events with the RK2 linear-interpolation helper.

• Combine the two histories without duplicating the burnout point and return them in the RK2 convention: times has shape (n_times,) and states has shape (n_times, 3). Inspect the event arrays rather than assuming events exist.

• Return a dictionary containing status ('impact', 'time_limit', or 'solver_failure'), times, states, total rhs_evaluations from both nfev values, both raw SciPy solution objects, burnout time and state, burnout_mass_residual, apogee time and state, impact time and state, maximum speed, and its time. Determine maximum speed from the combined stored states.

Cell 2 — tolerance study

• Use max_step = 1.0 s for every tolerance run. Use these four configurations, in this order:

  • rtol=1e-4, atol=np.array([1e-4, 1e-6, 1e-6]);

  • rtol=1e-6, atol=np.array([1e-6, 1e-8, 1e-8]);

  • rtol=1e-8, atol=np.array([1e-8, 1e-10, 1e-10]);

  • rtol=1e-10, atol=np.array([1e-10, 1e-12, 1e-12]).

• Run integrate_rocket_scipy() once per configuration. Retain one dictionary per run in scipy_tolerance_runs, including configuration, status, five quantities of interest, burnout-mass residual, and RHS evaluations. The five quantities are maximum speed, apogee altitude, apogee time, impact time, and positive impact-speed magnitude. A run is usable for this complete-flight comparison only when status == 'impact'; otherwise record None for all five quantities. Keep every status visible.

• Print one readable row per run. Do not claim that a configuration is adequate.

Cell 3 — comparison evidence

• Define targets of 0.1 m/s for maximum speed, 1.0 m for apogee altitude, 0.02 s for apogee time, 0.05 s for impact time, and 0.1 m/s for impact speed.

• For every adjacent pair of usable SciPy runs (status == 'impact'), compute the absolute change in each quantity and divide it by that quantity’s target. Compare every usable SciPy run with the matching quantity from rk2_accepted in the same way. Retain and print all records; do not turn these comparisons into pass/fail labels.

• Repeat the tightest tolerance configuration once with max_step=0.5 s and store it as scipy_max_step_check. Compare it with the tightest max_step=1.0 s result using the same quantities and targets.

• Plot normalized successive-tolerance changes against rtol for all five quantities on logarithmic axes, with a horizontal line at 1. Keep the numerical rtol ordering rather than reversing the axis.

• Do not use Richardson extrapolation on the tolerance sequence: tightening rtol by a factor of 100 does not define a uniform grid-refinement ratio or a known output-error order.

Allowed actions

Read the attached notebook and add only the three requested cells. Do not run code, alter existing cells, create or modify other files, access the network, install packages, or write my interpretation.

Done when

All three cells are inserted without outputs or execution counts. Stop and report where they were added, which functions and variables they create, and any ambiguity encountered. Do not report numerical results.

Invoke the agent¶

Save and attach your working notebook as in Lesson 4. Send this short request:

Do not ask the agent to run or debug the cells yet. First inspect what it proposed.

Audit the proposal before execution¶

Confirm that the agent changed only the requested location and that all three cells are unexecuted. Then trace the implementation against the brief:

• Do the powered and coast wrappers use constant, phase-specific mass-flow rates, including at $t_b$?

• Are there exactly two solve_ivp() calls, with the second initialized from the first call’s final state after the exact propellant reset?

• Is the burnout residual saved before the reset rather than hidden by it?

• Do both event functions have the correct component, direction, and terminal setting, and are they passed in the stated order?

• Does the code use SciPy’s event results rather than the RK2 interpolation helper?

• Are rtol, component-wise atol, max_step, and dense_output passed explicitly, with no t_eval?

• Is the duplicate burnout column removed when histories are combined, and is result.y transposed into the RK2 state-history convention?

• Does a run enter the five-quantity comparisons only when status == 'impact', with None recorded for every quantity otherwise?

• Are work, successive changes, RK2 differences, and the maximum-step check computed exactly as specified?

• Has the agent avoided Richardson extrapolation on the tolerance sequence and avoided choosing a configuration?

If a line is too compact to explain, ask the agent to rewrite only the affected cell using intermediate variables and comments. If you find a defect, describe the violated requirement and authorize a correction only to that cell. Record the exchange and inspect the revision before execution.

Test the interfaces, then run the study¶

Run only Cell 1. Before starting a trajectory, call the two phase wrappers once at $t_b$ with the same physically valid state. Confirm that the powered wrapper returns $dm_p/dt=-\mu_0$ while the coast wrapper returns $dm_p/dt=0$. Evaluate both event functions on simple states and confirm that each returns the intended component. These focused checks can expose a phase or state-index error before it is buried in a full trajectory.

Next run Cell 2. Require successful powered and unpowered solves, one apogee, one terminal impact, event states close to their defining zeros, a small burnout-mass residual, and physically ordered events. Inspect message on either raw solution object if the status is not impact.

Finally run Cell 3 and reproduce at least one normalized change and one RK2 comparison from the printed values. Choose the loosest tolerance configuration supported by all of the following:

• the next tighter configuration changes every quantity by less than its target;

• agreement with the accepted RK2 result is consistent with both calculations’ estimated numerical uncertainty;

• the tighter `max_step=0.5\ \mathrm{s}$ check does not change the conclusion; and

• all solver statuses, events, and physical checks remain valid.

A smaller rtol is not automatically a better engineering choice if it increases nfev without affecting any required quantity. Conversely, agreement between two solvers is not persuasive if they share the same incorrect model translation or mishandle the same breakpoint.

Check the coast against an analytical solution¶

The RK2 and SciPy calculations are independent numerical implementations, but both translate the same model. We can add a different kind of evidence after burnout. During the upward coast, the propellant is gone, the mass is the constant $m_s$, and $v>0$. Define

so the velocity equation becomes

Let $h_b$ and $v_b>0$ be the altitude and velocity at burnout. Separating variables and integrating from $v_b$ to zero gives the coast time to apogee:

The speed $v_T$ is the empty-shell terminal-speed magnitude already found in Equation 15. Using $dv/dt=v\,dv/dh$ gives the altitude gained during the coast:

These expressions are exact for the stated coast model. They are conditional on the supplied burnout state: they can expose a wrong coast equation, drag sign, phase switch, or apogee calculation, but they cannot reveal an error already present in $h_b$ or $v_b$.

Reconstruct exact_coast_apogee() in your working notebook and map each expression to Equation 25 or Equation 26. Then identify the tolerance configuration you accepted in the preceding study and run it once more to retain its complete result:

selected_scipy_index = _____
selected_configuration = (
    scipy_tolerance_runs[selected_scipy_index]['configuration']
)
scipy_accepted = integrate_rocket_scipy(
    u_0,
    selected_configuration['rtol'],
    selected_configuration['atol'],
    selected_configuration['max_step'],
    time_limit,
)
assert scipy_accepted['status'] == 'impact'

For both accepted calculations, use that method’s own burnout state as input to exact_coast_apogee(). Reconstruct and run this comparison cell:

coast_check_records = []
for method_name, numerical_result in (
    ('RK2', rk2_accepted),
    ('SciPy', scipy_accepted),
):
    exact_result = exact_coast_apogee(
        numerical_result['burnout_state'],
        numerical_result['burnout_time'],
        m_s, g, rho, A, C_D,
    )
    record = {
        'method': method_name,
        'apogee_time_difference': abs(
            numerical_result['apogee_time']
            - exact_result['apogee_time']
        ),
        'apogee_altitude_difference': abs(
            numerical_result['apogee_state'][0]
            - exact_result['apogee_altitude']
        ),
    }
    coast_check_records.append(record)
    print(record)

Reproduce at least one difference from the printed quantities, then compare both records with the apogee-time and apogee-altitude acceptance targets.

A disagreement points specifically to the numerical treatment between burnout and apogee. Agreement does not verify the powered phase, so also compare the two methods’ burnout altitudes and velocities and keep the earlier powered-flight checks in your final chain of evidence.

Write the capstone verdict¶

State whether the fixed-step RK2 and adaptive SciPy calculations support the requested predictions. Your verdict should include:

• the selected RK2 step and its Richardson evidence;

• the selected SciPy tolerances, max_step, and supporting tolerance change;

• the two methods’ reported maximum speed, apogee, and impact conditions;

• their differences relative to the numerical-error targets;

• the analytical coast-to-apogee differences and what that check cannot test;

• RHS evaluations for each selected calculation;

• how burnout and events were handled and checked;

• what agent-produced code you audited or corrected; and

• which uncertainties belong to numerical solution and which remain limitations of the physical model.

Do not report more digits than the accepted numerical uncertainty supports.