Accuracy, Cost, and Engineering Judgment

The previous lessons are a great foundation in numerical methods for studying dynamical systems governed by ordinary differential equations. You learned to apply Euler’s method and studied its rate of convergence using numerical experiments. An exercise on paper in Lesson 2 using Taylor expansions showed that convergence to be first order. We found numerical evidence of this behavior in Lesson 3 using the full nonlinear phugoid model.

Euler’s method, however, showed a clear weakness for the undamped oscillator in Lesson 2: every step artificially increased its oscillation amplitude. You can understand this graphically. Each Euler step advances a variable (for example, position) to the next time interval using:

Recall that the derivative of a function corresponds to the slope of the tangent at a point. Euler’s method uses the slope at the initial point in an interval, and advances the numerical position with that initial velocity. In the sketch below we exaggerate the result: the Euler estimate overshoots the exact solution and the amplitude of the oscillation grows unphysically, adding a small amount of energy with every step.

First-order methods in general require very small time steps to achieve acceptable accuracy, and thus they are rarely used for real-world dynamical calculations. We now investigate a method with a higher order of accuracy. Among the most popular higher-order methods are the Runge-Kutta methods, developed around 1900: more than 100 years after Euler published his book containing the method now named after him.

Paper-airplane challenge¶

Throughout this lesson, we use this motivating challenge: how far can you make a paper airplane fly?

Once released, a paper airplane has no propulsion or active control; its flight depends strongly on the launch (angle and speed) and aerodynamic response. If we use computation to recommend a launch, we need to distinguish a genuine improvement from a change caused by time-stepping error.

Our engineering question is:

Find a competitive launch within specified bounds, then determine which method supports its predicted range to the required accuracy with less computational work.

You will use the following common problem data:

• Model parameters: $C_L=1$, $C_D=0.2$ (so $L/D=5$), $v_t=4.9\ \mathrm{m/s}$, and $g=9.81\ \mathrm{m/s^2}$. The lift-to-drag ratio is motivated by measurements reported by Feng & others (2009).

• Release point: $x_0=0$ and $y_0=h=2\ \mathrm{m}$. The release height is fixed.

• Launch-search bounds: $4\leq v_0\leq12\ \mathrm{m/s}$ and $-30^\circ\leq\theta_0\leq30^\circ$. Convert angles to radians before using the model. These bounds define a model exercise, not experimentally validated launch limits.

• Required numerical range accuracy: $e_R=0.01\ \mathrm{m}$ (1 cm). This is a target for the numerical calculation within the model, not a claim of centimeter accuracy for a real paper airplane.

Define range as the net horizontal displacement $R=x_{\mathrm{ground}}-x_0$ at the first crossing from positive altitude to nonpositive altitude. A run that has not reached the ground by the common time limit of 15 s, or encounters an invalid state, does not supply a usable range. This challenge is adapted from the computational phugoid exercise of Simanca & Sutherland (2002).

The investigation has two stages:

1. Controlled comparison: keep the launch fixed at $v_0=6.5\ \mathrm{m/s}$ and $\theta_0=-0.1\ \mathrm{rad}$. Compare Forward Euler and explicit midpoint RK2 using the same ground-crossing procedure, then determine the work each needs to support the range-accuracy target. The comparison brief follows the convergence study.

2. Agent-supported launch investigation: search within the stated bounds of launch speed and angle, check that the best launches are not sensitive to the time step, test four nearby launches, then compare the methods for one competitive candidate. Treat ranges within 1 cm of the best tested range as equally useful for this exercise; the search does not prove a global optimum. The agent activity follows the controlled comparison.

Work through these stages in order. First derive and inspect the midpoint update, then reconstruct and verify the shared touchdown evaluator and complete the controlled comparison. In Stage 2, an AI agent may write the repetitive parameter-search code into your notebook, but you will inspect and execute each new cell yourself.

Runge–Kutta methods¶

A method’s order of convergence describes how the accumulated, or global, discretization error changes as we refine the time step over a fixed time interval. For a sufficiently smooth solution, a convergent method of order $p$ has global error ${\mathcal O}(\Delta t^p)$. When the leading discretization-error term dominates, we expect approximately

where $e$ measures the global error and $C$ depends on the problem, method, and time interval, but not on $\Delta t$. First-order error ($p=1$) scales linearly with the step size; second-order error ($p=2$) scales quadratically. These are small-step trends, not a guarantee that a higher-order method has a smaller error at every chosen step size.

One idea for improving on Euler’s method is to estimate the derivative at an intermediate point, like the midpoint, which results in the so-called explicit midpoint method or modified Euler method. The scheme has two steps and is written as:

Notice that this step evaluates the right-hand side, $f(u)$, twice: at the current state and at a predicted midpoint state. The notation $u_{n+1/2}$ represents an estimate of the state at $t_n + \frac{\Delta t}{2}$. All Runge–Kutta methods use such intermediate evaluations, called stages. Their order depends on how the stage states are constructed and how the derivative evaluations are combined—not simply on how many evaluations are performed.

Explicit midpoint is a second-order Runge–Kutta method (RK2). Its higher order does not eliminate all of Euler’s limitations. For the undamped oscillator in Lesson 2, sufficiently small time steps give less artificial amplitude growth than Forward Euler, but the amplitude still grows for any fixed, nonzero step size. Convergence as the step size shrinks over a fixed time interval is different from stable long-time behavior at a fixed step size.

An interesting historical connection to our phugoid problem: Carl Runge’s daughter Iris—an accomplished applied mathematician in her own right—worked assiduously over the summer of 1909 to translate Lanchester’s “Aerodonetics.” She also reproduced his graphical method to draw the phugoid curves Tobies, 2012, p. 73.

Phugoid model with second-order RK¶

Let’s begin the paper-airplane investigation by computing a baseline flight under the full phugoid model using both Forward Euler and second-order Runge–Kutta. We will use the same launch conditions for both methods and examine the horizontal distance traveled before the airplane touches the ground.

We will first compare the methods over a common interval while the airplane is still aloft. That fixed-time study checks the behavior of RK2, but it cannot yet justify a touchdown range or a time step for the engineering challenge. Afterward, we will build and verify the event calculation needed to measure range.

Start by importing NumPy and Matplotlib with the aliases used in the previous lessons. We also set the font family and size through Matplotlib’s rcParams dictionary.

Use the challenge’s lift-to-drag ratio $L/D=5.0$ and trim speed of $4.9\ \mathrm{m/s}$. What do you think will happen if you make $L/D$ higher?

The initial speed is a little higher than the trim speed, the launch angle is negative, and the release height is 2 meters. We will use the same initial state and parameters for both numerical methods.

Reuse functions from a Python module¶

We have already written and inspected these three functions in Lesson 3:

• rhs_full_phugoid() gives the four model derivatives.

• euler_step() advances the complete state by one Forward Euler step, in the *args form.

• discrete_l1_difference() compares histories on nested time grids.

These functions are now reused often enough to save in a separate file. A module is an ordinary Python source file with the extension .py. To make one, create a plain-text file, copy the function definitions into it, and include the imports those definitions need—these functions only need import numpy as np. Keep parameter choices, integration loops, plots, and notebook-only commands in the notebook.

The course provides these definitions in a file named phugoid.py. Read the module source: its three functions are the ones already presented in previous lessons. We can download this file directly using urlretrieve() from Python’s standard library; add the code below to your notebook and execute it.

Here, url points to the raw Python file, and fname names the local copy. The file is saved in the kernel’s current working directory, normally the folder containing your notebook. Keep phugoid.py alongside your working notebook. You only need to download it once; rerunning the download overwrites that file, so save a separate copy of any local edits first.

Downloading saves the file; importing makes its functions available in the notebook. In the import statement, use the filename without .py. Python executes a module’s top-level statements when it first loads it, so you should only import code from sources you trust and have inspected. Our file imports NumPy and defines functions; it does not run a simulation.

Functions in the module do not inherit notebook variables, so we continue to pass the state and model parameters explicitly. The difference function requires both time-step sizes to check grid nesting and endpoint alignment.

If you edit the local module, restart the kernel and rerun your imports and calculations, skipping the download cell to preserve your edits. Rerunning an import alone does not reload an already imported module. For more detail, see the Python modules tutorial.

Define the RK2 step¶

The reused functions now come from the module, but the new numerical method stays visible here and in your notebook, as it is the first time you encounter it. Define rk2_step() to implement the modified Euler method in Equation 3, also known as second-order Runge–Kutta or RK2. The time loop will call this function once per step.

Compare on a common time interval¶

Begin with an illustrative time step of $\Delta t=0.01\ \mathrm{s}$. We have not yet shown that this step is accurate enough for touchdown range, so do not treat it as an engineering recommendation. Integrate only to $T=2\ \mathrm{s}$, when the baseline flight is still above ground, and compare the two approximations on exactly the same time grid.

As in Lesson 3, the variable num_steps counts updates; each history has num_steps + 1 rows to include the initial state. Both arrays have four columns ordered as [v, theta, x, y].

Extract time, horizontal position, and altitude for the exploratory plot. The final altitudes will also confirm that this entire comparison ends before impact.

Plot for exploration, not an accuracy test¶

A trajectory plot is a useful first diagnostic: it can expose a wrong sign, a discontinuity, or an implausible path. It cannot show that a range is accurate to 1 cm. Two curves can overlap visually while differing by more than the required tolerance.

For now, inspect the two paths and their separation. This is exploratory evidence, not a verdict.

Fixed-time trajectory convergence¶

Just like in Lesson 3, we want to check whether RK2 exhibits its expected convergence rate on this problem. Every history below ends at the same pre-impact time, $T=2\ \mathrm{s}$, so values at matching indices describe the same physical times.

The for loop computes the solution on several nested time grids, with the coarsest and finest step sizes differing by a factor of 100. We then compare horizontal position, a single quantity with units of meters. Mixing all four state columns would combine speed, angle, horizontal position, and altitude in one number with no coherent physical units.

For the general nonlinear trajectory computed here, we do not have an exact reference solution available. The finest-grid history is a numerical reference, not an exact solution, so the differences below are not exact errors. They are evidence about refinement behavior over this fixed interval.

Once those runs are complete, compare each horizontal-position history with the finest-grid history. Pass both time-step sizes to the imported discrete_l1_difference() function so it can align the grids.

Plot the differences against time-step size on logarithmic axes, as in the previous lessons.

The decreasing differences suggest convergence, but their size is measured relative to a numerical reference. In Lesson 3, the observed order for Euler’s method was close to 1. For RK2, the expectation is $p\approx2$ once the leading discretization-error term dominates. Let us test that expectation.

To compute the observed order of convergence, we use three grid resolutions that are refined at a constant rate, in this case $r=2$.

An observed order close to 2 is consistent with the expected second-order behavior for this grid family. In the regime where the leading error is proportional to $\Delta t^2$, halving the step size reduces that error to approximately one quarter of its previous size. The raw differences and their ratio make that conclusion inspectable rather than reporting only the final value of $p$.

This study concerns horizontal-position histories over a fixed, airborne interval. It does not verify touchdown range or justify $\Delta t=0.01\ \mathrm{s}$ for the challenge. Touchdown occurs at a step-dependent time, and locating it introduces a separate event-calculation error. We need to define and test that measurement before studying range convergence.

Locate touchdown and evaluate range¶

Touchdown is an event: it happens when the altitude first crosses the ground, not necessarily at one of our stored times. Starting from positive altitude, one numerical step brackets touchdown when

Assume the state changes linearly between those two computed endpoints. Let $\alpha$ be the fraction of the step needed to reach $y=0$. Linear interpolation gives

Because the two altitudes bracket zero, $0<\alpha\leq1$. The horizontal component of $u_{\mathrm{g}}$ supplies the interpolated touchdown position, and the net range is $R=x_{\mathrm{g}}-x_0$. We retain the two computed endpoints as the crossing bracket, but stop immediately: no later below-ground states are calculated.

A usable result also needs an explicit completion status. We will report touchdown, time_limit, or invalid_state; the third arises when a calculation breaks, and we will see below why it needs a status of its own. Failed or unfinished runs receive no range.

Write the interpolation as a function¶

The crossing formula is small enough to be its own function, so that we can test it by hand before it is buried inside a loop. locate_touchdown() receives the two times and states that bracket the crossing and returns the interpolated touchdown time and state.

Check the function against your hand calculation before using it. The bracket below uses the values from the On paper exercise; the expected values come from that calculation, not from the code.

Start with a direct loop¶

Before adding any checks, write the range calculation directly: take a step, test for a crossing, interpolate, and stop. If the flight is still aloft when the steps run out, report that instead. A for loop over int(np.floor(time_limit / dt)) takes only complete steps within the safety cap; ending less than one step early is harmless here. The function returns a status and a range, and the status time_limit means that no crossing was found.

Apply this function to the baseline launch, then to a launch with zero initial speed. The model contains $g/v$, so zero speed is outside its domain. Predict what the function will report in the second case before running the cell.

A silent wrong answer¶

The first result is the interpolated baseline range. The second is wrong in a way that is easy to miss. NumPy printed a warning about division by zero, but a warning does not stop a calculation. The first step produced an infinite trajectory angle, and the following step turned the state into nan, which stands for “not a number.” Every comparison with nan is False, including u_next[3] <= 0.0, so the crossing test never fires and the loop runs to the time limit. The draft function then reports that the airplane never landed, when in fact the calculation broke on its first step. Confirm the behavior of the comparison:

A result that looks like a legitimate outcome, time_limit, but comes from a broken calculation is a type of failure that needs to be prevented with defensive checks. The remedy is a validity test on every completed step, placed before the crossing test, so that a broken state is reported as such rather than silently compared. A state is valid when every entry is finite and the speed is positive: negative speed is outside the state convention, and zero speed is outside the model’s domain. A nonpositive altitude reached from a positive altitude is the event, not an invalid state.

We can write this test as a small function and check it on a valid state, the zero-speed state, and a state containing nan.

Map the possible outcomes¶

Before any more coding, you should understand the possible outcomes. The complete evaluator has four jobs:

1. Check the setup. A nonpositive time step, an initial state with the wrong number of entries, or a nonpositive initial altitude is a caller error: no crossing from positive altitude could ever be found, and the draft function above would silently report time_limit. The full function raises ValueError rather than starting a calculation with a broken setup.

2. Advance and count steps. The evaluator counts accepted completed steps, including the step that brackets touchdown. A step that first exposes an invalid state is rejected and not counted; failed runs do not enter the work comparison. For successful runs, multiplying by one right-hand-side evaluation per Euler step or two per midpoint step gives the computational work.

3. Recognize a failed flight calculation. A completed step whose state is non-finite or has nonpositive speed produces the status invalid_state and no range. The check inspects completed steps only; a stricter evaluator would also inspect RK2’s intermediate stage.

4. Recognize completion. The first positive-to-nonpositive altitude bracket produces touchdown; using up the steps without such a bracket produces time_limit. Only touchdown supplies a range.

Here is the same control flow as pseudocode. Read this outline first, then find each part in the Python function below.

check the setup

for each step up to the time limit:
    take one numerical step
    if the new state is invalid:
        stop with invalid_state
    store the new state
    if altitude crossed zero:
        interpolate touchdown
        stop with touchdown

report time_limit

At the start of the function, np.asarray(u_0, dtype=float) normalizes a list or array into floating-point NumPy data. We call that working state u because it changes as the integration advances, while u_0 names the original initial condition. The initial horizontal position is saved separately so the final range can always be measured from the launch point.

Check the non-success outcomes¶

A failure path deserves a test just as much as a successful calculation. Use a deliberately short time limit to produce an unfinished run, then supply the zero-speed launch that fooled the first draft of the evaluator. In both cases, inspect the status and confirm that range is unavailable. The zero-speed run now rejects its first attempted step instead of running to the time limit; its reported steps value is therefore zero.

The initial-state guard also deserves a test. The cell below deliberately passes a release point below ground (-1.0) and confirms that the function refuses to start. It is a small test asking: “Does the function refuse to simulate an aircraft that already starts below ground?” The draft evaluator would have treated the first update as touchdown because its altitude was already nonpositive. With both endpoints below ground, locate_touchdown() would extrapolate to a fictitious crossing outside the step rather than interpolate a real bracket.

The invalid input is created with -1.0 on the last element of u_below_ground. The test begins by assuming the guard has not run:

guard_triggered = False

Then it tries to call the function:

try:
    integrate_until_touchdown(...)

Inside the function, this condition detects the negative altitude:

if u.shape != (4,) or u[3] <= 0.0:
    raise ValueError(...)

Normally that would stop the cell. However, the except block catches that particular error, and guard_triggered = True records that the expected rejection occurred. The final assert checks that the function really did reject the input.

The complete flow is:

Negative altitude supplied
        ↓
Function raises ValueError
        ↓
except catches it
        ↓
guard_triggered becomes True
        ↓
assert passes

Reread this little test and explanation as many times as you need. Check your understanding by answering: what would happen if the altitude check were accidentally removed?

Verify the evaluator with steady glide¶

The general nonlinear flight lacks an exact solution we can use to verify the evaluator, but the model does have an exact steady-glide special case.

Use $[v_s,\theta_s,0,h]$ as a separate verification initial state, not as a replacement for the baseline launch. Because the vertical velocity is constant, the exact touchdown time is $t_{\mathrm{g}}=-h/(v_s\sin\theta_s)$. Both Euler and RK2 advance this straight-line motion exactly apart from floating-point effects.

The step sizes below place touchdown strictly between stored times. For every trial, require the touchdown status, confirm that the interpolated time lies inside the final bracket, and compare the range with Equation 8. This checks the event calculation; it does not establish accuracy for a curved trajectory or validate the physical model.

The assertions are deliberately demanding because this equilibrium produces straight-line motion: the time-stepping methods and the linear event model are exact for this case apart from floating-point rounding. Passing this test at several off-grid touchdown times gives us a focused check of the bracketing and interpolation logic. It does not tell us how small $\Delta t$ must be for the curved baseline flight; that requires range refinement.

Compare methods at a required accuracy¶

This is Stage 1 of the paper-airplane challenge. Keep the baseline launch and all model parameters fixed. The question is not which method runs faster on the same grid, but which needs fewer right-hand-side evaluations to support a range accurate to 1 cm.

We now have a verified range evaluator with explicit completion statuses and a common event treatment for both methods. Our next goal is to apply it first to the baseline launch at the illustrative step size, then establish a numerical reference by refinement.

Apply the evaluator to the baseline launch¶

Return to the baseline initial state and use the illustrative step size $\Delta t=0.01\ \mathrm{s}$. The next cell demonstrates the evaluator’s interface and confirms that each successful run stops at its first crossing bracket. The printed ranges are interpolated, but they are still provisional numerical results: we have not yet established their errors by refinement.

For a successful run, states contains the two endpoints that bracket touchdown. The plotting cell keeps the last valid computed point but replaces the nonpositive endpoint with touchdown_state, so the displayed physical path ends at $y=0$. Work is counted in right-hand-side evaluations: the completed steps, including the bracketing step, multiplied by one evaluation per Euler step or two per midpoint step. The dictionary rhs_calls_per_step records those factors so that the conversion stays visible.

The two methods are not required to use the same time step in the final engineering comparison. Equal $\Delta t$ is a useful diagnostic, but it gives RK2 about twice the computational work while its computational accuracy is superior. Our primary question is thus: what is the least RHS work each method needs to support the same range accuracy?

Establish a numerical range reference¶

The baseline ranges at one time step are provisional. Without an exact solution for the curved flight, the reference is a refined RK2 result, accepted when the range changes by less than one tenth of the accuracy target, 1 mm, in each of two successive step halvings. A reference must be resolved well beyond the tolerance it will be used to judge; this rule provides evidence of that, not a rigorous error bound. The cell below applies the rule starting from $\Delta t=0.01\ \mathrm{s}$ and stops at the first step size that satisfies it, or reports insufficient evidence after eight halvings. It also adds up the work spent on the reference, which is reported separately from the comparison itself, and it stores the accepted range as range_reference and its step size as dt_search, the names used by the comparison below and by Stage 2.

The two successive step halvings changed the RK2 range by less than 1 mm each, so we accept the last result as our working numerical reference. The cell stores that range as range_reference and its step size as dt_search; use these variables, rather than retyped values, in the accuracy-and-work comparison and in Stage 2.

In your notebook, also refine Euler and check that its ranges move toward range_reference, continuing until their difference is well inside the 1 cm target. Euler is first-order, so expect more halvings. The results need not agree in every printed digit: both retain discretization error. Do not assign to range_reference or dt_search in the Euler check, so that the RK2 values are not overwritten.

Plan the accuracy-and-work comparison¶

Which method needs fewer right-hand-side (RHS) evaluations for a range accurate to 1 cm? Reuse your existing functions, keeping the baseline launch, model parameters, touchdown interpolation, and 15 s time limit fixed.

1. Keep one shared reference. Use range_reference, the accepted RK2 range, as $R_{\mathrm{ref}}$ for both methods; your Euler cross-check plays no part in it.

2. Build a table of runs. For each method, start at $\Delta t=0.1\ \mathrm{s}$ and halve the step at least three times: 0.05, 0.025, and $0.0125\ \mathrm{s}$. Continue halving as needed. Starting coarse reveals where each method becomes accurate enough.

3. Find and check a qualifying run. A run qualifies when $|R-R_{\mathrm{ref}}|\leq0.01\ \mathrm{m}$, the accuracy target. After finding one, halve the step at least once more and check that the finer run still qualifies and its difference from the reference decreases. Otherwise, investigate and refine further.

4. Compare work. For each method, select the least-work qualifying run supported by the finer-step check. The selected runs may use different time steps.

Record one row per run, replacing the dashes and extending this table as needed:

For each row, compute work as result['steps'] * rhs_calls_per_step[method_name]. Count only that run. For failed or unfinished runs, record None for both range and difference, keep status and work, and exclude them from selection. Report the accumulated RK2 reference work and Euler cross-check work separately below the table.

Plot the absolute difference from the reference against RHS evaluations, using a logarithmic scale on both axes and a horizontal line at $0.01\ \mathrm{m}$. These differences are not exact errors.

Conclude with each method’s selected step size, work count, and supporting refinement. State which uses less work and by what factor, for this baseline and tested steps. If evidence is insufficient, say so.

Together with your earlier derivations, touchdown verification, and reference checks, this table, plot, and conclusion complete Stage 1.

Investigate launches with an agent¶

We return to the motivating question: which launch flies farthest, and which method supports its predicted range with less work? You now have a verified touchdown evaluator, an accepted search step dt_search, and a procedure for comparing the two methods at 1 cm accuracy. What remains is repetitive: hundreds of runs over a grid of launch speeds and angles. An agent will write that code into your notebook. You will read each cell before running it, run it yourself, and decide what the results mean.

Two ideas guide the whole activity. Keep them in mind while you read the agent’s output.

Treat launches within 1 cm of the best tested range as tied for our recommendation. We choose not to distinguish smaller improvements for this paper-airplane exercise. This is an engineering judgment about useful precision, even when the numerical calculation resolves smaller differences; it does not establish the accuracy of the model for a real airplane.

A time step is validated for a flight, not for a method. Stage 1 found the time step each method needs for the 6.5 m/s launch. The best launches fly faster and farther, so you can expect they could need smaller steps. The last part of the activity measures this.

The search uses dt_search, the step accepted for the range reference, rather than the coarser step that met the 1 cm target in Stage 1. This better-resolved step helps keep numerical error small compared with our 1 cm decision threshold.

Record the launch-search brief¶

Copy the brief below into a Markdown cell in your notebook and fill in the two bracketed entries. Above it, note the date and the agent or persona you use.

Agent launch-search task brief¶

Goal and scope

Add two readable, unexecuted Python cells below this brief in [notebook filename]. Cell 1 searches a grid of launches. Cell 2 checks the leading launches. Leave the choice of launch and all conclusions to me.

Non-negotiables

Reuse rk2_step(), rhs_full_phugoid(), integrate_until_touchdown(), the model parameters, the release point, and the 15 s time limit already in the notebook. Use the notebook variable dt_search, whose value is [accepted RK2 step]. Use ordinary loops and comments. Report angles in degrees and convert to radians only when building a state.

• Cell 1, the search:

  • Run RK2 at dt_search for all 17 equally spaced speeds from 4 to 12 m/s and all 25 equally spaced angles from $-30^\circ$ to $30^\circ$, endpoints included: 425 launches.

  • Keep one record per launch: speed, angle in degrees, range, status, and step count. A failed or unfinished run keeps range=None and is never ranked.

  • Print the count of each status, then the five largest ranges with their speed and angle. Keep all records in a list for Cell 2. If fewer than five runs succeed, print that instead and stop.

• Cell 2, the checks:

  • Time step: rerun the five best launches at dt_search/2. For each, print the range and status at both steps and the change between successful results.

  • Nearby launches: around the best launch from Cell 1, change only the angle by $-1.25^\circ$ and $+1.25^\circ$, then only the speed by -0.25 and +0.25 m/s. Skip any point outside the original bounds. Run these checks at dt_search/2. Print each range and its signed change from the center launch at that same step; positive means farther.

  • Retain all check records and statuses. Failed or unfinished checks keep None for range and change.

  • Print results only. Do not rank, choose, or comment on them.

Allowed actions

Read the attached notebook and add only these two cells. The following are disallowed: changing existing cells, redefining the model or numerical functions, running code, creating files, downloading anything, accessing the network, or installing packages.

Done when

Both cells are inserted, without outputs or execution counts. Stop and report where they were added and which variables they create. I will inspect and run them.

Invoke the agent¶

Save the notebook with your completed brief and attach it to your agent interface. If the interface lets you allow edits but not execution, set it that way. Send this short request:

Audit and run the search¶

Run Cell 1. Check that the status counts add up to 425. Look at the five best ranges and apply the tie rule: which of them lie within 1 cm of the best? Expect the best few to tie. If many runs failed, or the results look strange, stop and find out why before going on.

Run Cell 2. Check that each time-step change is below reference_tolerance, and inspect whether the order changes. If the changes are too large, halve again for the same launches before relying on the comparison. The four nearby launches show how small changes in speed or angle affect range. They may improve on the grid best, but they do not map a boundary or establish what happens at untested launches.

If a cell fails with an error, send the error message to the agent and allow it to change only that cell. Read the revised cell before running it, and note what changed.

v_candidate = _____            # m/s
theta_candidate_deg = _____    # degrees
u_candidate = np.array([
    v_candidate, np.deg2rad(theta_candidate_deg), x_0, y_0
])

print('dt (s)         range (m)    change (m)')
range_previous = None
reference_work_candidate = 0
range_reference_candidate = None
for dt_trial in (dt_search, dt_search / 2, dt_search / 4):
    result = integrate_until_touchdown(
        rk2_step, u_candidate, rhs_full_phugoid,
        dt_trial, time_limit, C_L, C_D, g, v_t,
    )
    assert result['status'] == 'touchdown'
    reference_work_candidate += result['steps'] * rhs_calls_per_step['RK2']
    if range_previous is None:
        print(f'{dt_trial:<12.6f} {result["range"]:12.6f}')
    else:
        change = abs(result['range'] - range_previous)
        assert change < reference_tolerance
        print(f'{dt_trial:<12.6f} {result["range"]:12.6f} {change:12.3e}')
    range_previous = result['range']

range_reference_candidate = range_previous
print(f'Candidate reference work: {reference_work_candidate} RHS evaluations')

Record the comparison brief¶

Now repeat the Stage 1 method comparison for your candidate launch. Copy the brief below into a Markdown cell under the cell that defines u_candidate. It has no blanks.

Do not expect the Stage 1 time steps to pass again. Your launch flies faster and farther than the 6.5 m/s launch, so both methods could need smaller steps, and a coarse Euler run may fail with invalid_state. Keep a failed run in the table and investigate what caused it.

Agent candidate-comparison task brief¶

Goal and scope

Add one readable, unexecuted Python cell immediately below this brief. It applies a method comparison (Euler vs. midpoint Runge-Kutta) to u_candidate, using range_reference_candidate as the reference. Report the table; leave conclusions to me.

Non-negotiables

• Reuse the existing model, step functions, touchdown evaluator, parameters, release point, and 15 s time limit. Use ordinary loops and comments. Leave range_reference, dt_search, and range_reference_candidate unchanged.

• For each method, start at a time step of 0.1 s and halve it repeatedly, at most 12 times. A run qualifies when its computed range differs from range_reference_candidate by at most range_target. Stop halving the time step after the first qualifying run and one finer run that also qualifies with a smaller difference. If that does not happen within 12 halvings, report insufficient evidence for that method.

• Print one row per run: method, time step, range, absolute difference from the reference, completed steps, RHS evaluations, and status. A failed or unfinished run keeps its status and work, shows None for range and difference, and never qualifies.

• Count computational work as result['steps'] * rhs_calls_per_step[method_name].

Allowed actions

Read the attached notebook and add only this cell. Disallowed: change existing cells, redefine numerical functions, run code, create files, download anything, access the network, or install packages. Do not say which method is better.

Done when

The cell is inserted, without output or execution count. Stop and report where it was added and which variables it creates. If anything in this brief is unclear, say so instead of guessing.

Invoke, audit, and run the comparison¶

Save the notebook, attach it as before, and send:

Run the cell. Pick one successful row and check its RHS count and its difference from the reference by hand from the printed values. Then set the table beside your Stage 1 table: do the Stage 1 time steps still qualify for this launch?

Reflect on a coarse Euler failure

For the launch at 10 m/s and $-17.5^\circ$, Euler at 0.1 s turns the computed trajectory upward and then produces a nonpositive speed. The evaluator reports invalid_state because the model requires positive speed and contains $g/v$.

Compare this outcome with the finer Euler runs and RK2. Does the failure persist under refinement? Explain why its disappearance supports a numerical artifact. This model uses fixed lift and drag coefficients and has no stall criterion; invalid_state is not a prediction of a physical stall.

Completion check¶

Stage 2 is complete when your notebook contains both briefs, the search and check results with your audit notes, the competitive launches and your candidate with its reference range, the comparison table, and your verdict with its agent record. Each numerical claim must point to the run that supports it.

Why we work this way¶

We wrote the touchdown checks by hand so that status, event location, and work counting would be understandable. Engineers rarely rewrite such checks for every study: they reuse a tested tool, and a coding agent may draft the repetitive experiment code, as it did here. What does not change is who specifies the runs, reads the code, checks representative results, and owns the conclusion. The lessons from this activity carry over to any computational study:

• numerical refinement supports the chosen accuracy,

• engineering judgment determines which differences matter, and

• a time step validated for one case must be checked again when the case changes.