Phugoid Oscillation

In the first lesson, Phugoid Motion, we described the physics of a glider’s oscillatory trajectory, seen as an exchange of kinetic and potential energy. This analysis goes back to Frederick Lanchester, who published his book “Aerodonetics” on aircraft stability in 1908. We concluded that first exposure to our problem of interest by plotting the flight paths predicted by Lanchester’s analysis, known as phugoids.

Here, we will look at the situation when an aircraft is initially moving on the straight-line phugoid (obtained with the parameters $C=2/3$, $\cos\theta=1$, and $z=z_t$ in the previous analysis), and experiences a small upset, a wind gust that slightly perturbs its path. It will then enter into a gentle oscillation around the previous straight-line path: a phugoid oscillation.

As in the first lesson, $z$ is a depth measured positive downward from the chosen energy-reference level, and $z_t$ is its equilibrium value for trimmed flight. Thus $z$ is not the aircraft’s altitude: an upward acceleration is $-d^2z/dt^2$. If we assume that the perturbation is small, then $\cos\theta=1$ is a good approximation and Newton’s second law in the vertical direction is:

We previously saw that the following relation holds for the ratio of lift to weight, in terms of the trim velocity $v_t$:

This will be useful: we can divide Equation 1 by the weight and use Equation 2 to replace $L/W$. Another useful relation from the previous lesson expressed the conservation of energy (per unit mass) as $v^2 = 2 gz$. With this, Equation 1 is rearranged as:

Look at Equation 3 for a moment. Does it ring a bell? Do you recognize it?

If you remember from your physics courses the equation for simple harmonic motion, you should see the similarity!

Take the case of a simple spring. Hooke’s law is $F=-kx$, where $F$ is a restoring force, $x$ the displacement from a position of equilibrium and $k$ the spring constant. This results in the following ordinary differential equation for the displacement:

which has the solution $x(t) = A \cos(\omega t- \phi)$, representing simple harmonic motion with an angular frequency $\omega=\sqrt{k/m}=2\pi f$ and phase angle $\phi$.

Now look back at Equation 3: it has nearly the same form and it represents simple harmonic motion with angular frequency $\omega=\sqrt{g/z_t}$ around the equilibrium depth $z_t$.

Think about this for a moment ... we can immediately say what the period of the oscillation is: exactly $2 \pi \sqrt{z_t/g}$ — or, in terms of the trim velocity, $\pi \sqrt{2} v_t/g$.

This is a remarkable result! Think about it: we know nothing about the aircraft, or the flight altitude, yet we can obtain the period of the phugoid oscillation simply as a function of the trim velocity. For example, if trim velocity is 200 knots, we get a phugoid period of about 47 seconds—over that time, you really would not notice anything if you were flying in that aircraft.

Next, we want to be able to compute the trajectory of the aircraft for a given initial perturbation. We will do this by numerically integrating the equation of motion.

Prepare to integrate¶

We want to integrate the differential equation and plot the trajectory of the aircraft. Are you ready?

Equation 3 is a second-order ordinary differential equation (ODE). Let’s represent the time derivative with a prime and write it like this:

There’s a convenient trick when we work with ODEs: we can turn this 2nd-order equation into a system of two 1st-order equations introducing an intermediate variable for the first derivative. Like this:

Here $b=z'$ is the rate of change of depth: $b>0$ means downward velocity, while $b<0$ means upward velocity.

Another way to look at a system of two 1st-order ODEs is by using vectors. You can make a vector with the two state variables,

and write the differential system as a single vector equation:

If you call the right-hand-side $\vec{f}(\vec{u})$, then the equation is very short: $\vec{u}'(t) = \vec{f}(\vec{u})$—but let’s drop those arrows to denote vectors from now on, as they are a bit cumbersome: just remember that $u$ and $f$ are vectors in the phugoid equation of motion.

Next, we’ll prepare to solve this problem numerically.

Initial value problems¶

Let’s step back for a moment. Suppose we have a first-order ODE $u'=f(u)$. You know that if we were to integrate this, there would be an arbitrary constant of integration. To find its value, we do need to know one point on the curve $(t, u)$. When the derivative in the ODE is with respect to time, we call that point the initial value and write something like this:

In the case of a second-order ODE, we already saw how to write it as a system of first-order ODEs, and we would need an initial value for each equation: two conditions are needed to determine our constants of integration. The same applies for higher-order ODEs: if it is of order $n$, we can write it as $n$ first-order equations, and we need $n$ known values. If we have that data, we call the problem an initial value problem.

Remember the definition of a derivative? The derivative represents the slope of the tangent at a point of the curve $u=u(t)$, and the definition of the derivative $u'$ for a function is:

If the step $\Delta t$ is already very small, we can approximate the derivative by dropping the limit. We can write:

With Equation 11, and because we know $u'(t)=f(u)$, if we have an initial value, we can step by $\Delta t$ and find the value of $u(t+\Delta t)$, then we can take this value, and find $u(t+2\Delta t)$, and so on: we say that we step in time, numerically finding the solution $u(t)$ for a range of values: $t_1, t_2, t_3 \cdots$, each separated by $\Delta t$. The numerical solution of the ODE is simply the table of values $t_i, u_i$ that results from this process.

Discretization¶

In order to execute the process described above and find the numerical solution of the ODE, we start by choosing the values $t_1,t_2,t_3 \cdots t_n$—we call these values our grid in time. The first point of the grid is given by our initial value, and the small difference between two consecutive times is called the time step, denoted by $\Delta t$. The solution value at time $t_n$ is denoted by $u_n$.

Let’s build a time grid for our problem. We first choose a final time $T$ and the time step $\Delta t$. In code, we’ll use readily identifiable variable names: T and dt, respectively. With those values set, we calculate num_steps, the number of updates needed to reach $T$. Because the grid includes both the initial and final times, it contains num_steps + 1 points.

Let’s write some code. The first thing we do in Python is load two scientific-Python libraries: NumPy for numerical functions and arrays, and Matplotlib for plotting.

Now initialize T and dt, calculate num_steps, and build a NumPy array containing the num_steps + 1 time points in the grid.

We have our time grid! Now it’s time to apply the numerical time stepping represented by Equation 11.

Euler’s method¶

The approximate solution at time $t_n$ is $u_n$, and the numerical solution of the differential equation consists of computing a sequence of approximate solutions by the following formula, based on Equation 11:

Equation 12 is called Euler’s method.

Applying Equation 12 to the phugoid system in Equation 6 gives the following algorithm that we need to implement in code:

And solve!¶

To implement Equation 13, we need to set things up in code: define the parameter values needed in the model, initialize a NumPy array to hold the two state variables, and initialize another array for the depth-coordinate values.

Now we can step in time using Euler’s method. range(num_steps) supplies the indices from 0 through num_steps - 1; each pass advances both state variables and stores the new depth at index n + 1.

Make sure you understand what this code is doing. This is a basic pattern in numerical methods: iterations in a time variable that apply a numerical scheme at each step.

Plot the vertical displacement¶

If the code is correct, we have stored the downward-positive depth coordinate in the array z. For a picture of the aircraft’s motion in the sky, however, we expect the vertical axis to point upward. We therefore define the vertical displacement from trim as

Equation 14 defines $\eta$ as positive above the trim level and negative below it. This change does not alter the numerical solution: Euler’s method still advances $z$ and $b$, and we convert from $z$ to $\eta$ only for visualization. In particular, a positive $b_0$ makes $z$ increase and $\eta$ decrease, so a downdraft appears as downward motion on the plot.

You should explore the Matplotlib Pyplot tutorial (if you need to) and familiarize yourself with the tools that control the size, labels, line style, and so on. Creating good plots is a useful skill: it is about communicating your results effectively.

Here, we set the figure size, vertical limits, grid, and axis labels. The final ax.plot() call draws a continuous red line.

Exact solution¶

The equation for phugoid oscillations is a 2nd-order, linear ODE and it has an exact solution of the following form:

where $A$ and $B$ are constants that we solve for using initial conditions.

Our numerical solution used the initial conditions:

Applying the initial conditions in Equation 16 to Equation 15 and solving for $A$ and $B$ gives:

We already defined all of the variables in Equation 17, so we can immediately compute the exact depth coordinate $z_{\rm exact}$. We then apply the sky-view transformation in Equation 14, $\eta_{\rm exact}=z_t-z_{\rm exact}$, for plotting.

Compare with the exact solution¶

Now we can plot both the numerical and exact vertical displacements to see how well Euler’s method approximated the phugoid oscillation. Both curves use $\eta$, so upward and downward on the graph match the aircraft’s motion in the sky.

To add another curve to a plot, call ax.plot() a second time on the same Axes. Giving each curve a label lets ax.legend() identify them.

The two curves agree fairly well at first, but the numerical oscillation grows toward the end. That growth is a numerical artifact: the ideal phugoid is undamped, so its exact amplitude remains constant. Re-run the previous steps with a smaller time step, say $dt=0.01$, and notice that the artificial growth becomes slower.

This brings up two different features of numerical methods. A method is convergent if, over a fixed time interval, its numerical solution approaches the exact solution as the time-step size tends to zero. If the global error behaves like $E=O(\Delta t^p)$, the method has order $p$; Forward Euler is first-order convergent.

Consistency is a related local property. Substitute the exact solution into one step of the numerical method and examine the residual. For Forward Euler, the unnormalized one-step defect is $O(\Delta t^2)$, so the residual per unit time is $O(\Delta t)$ and tends to zero. Forward Euler is therefore consistent. Consistency alone does not guarantee convergence: the method must also control how errors grow from one step to the next.

Convergence on a fixed interval does not guarantee faithful behavior over arbitrarily long times. For this undamped oscillator, Forward Euler adds a small amount of numerical energy at every step, causing the amplitude to grow. Reducing $\Delta t$ slows that growth, but no nonzero time step eliminates it completely. The refinement study below asks the fixed-interval question: does the global error tend to zero, and at what rate, as $\Delta t$ decreases?

Convergence¶

To compare the two solutions, we need to use a norm of the difference, like the $L_1$ norm, for example.

Equation 20 sums the individual differences between the exact and numerical solutions at the mesh points. In other words, $E$ is a discrete representation of the integral over the interval $T$ of the absolute difference between the computed $z$ and $z_{\rm exact}$:

Although Equation 21 is written using the depth coordinate $z$, we plotted the upward-positive displacement $\eta$. Subtracting the same $z_t$ and reversing the sign does not change an absolute difference: $|\eta-\eta_{\rm exact}|=|z-z_{\rm exact}|$.

We check for convergence by calculating the numerical solution using progressively smaller values of dt. We already have most of the code that we need. We just need to add an extra loop and an array of different $\Delta t$ values to iterate through.

Calculate the error¶

We now have a numerical depth-coordinate solution for each $\Delta t$ in the list z_solutions. The list method .append() added each completed NumPy array without requiring us to know its length in advance. To calculate the error corresponding to each $\Delta t$, we can write a function.

Now, we iterate through each $\Delta t$ value and calculate the corresponding error. In the following code cell, we use the built-in function zip to pair each array in z_solutions with its time-step size in dt_values. The trial-specific names keep this loop from overwriting the baseline variables t, dt, and z_exact used earlier in the lesson.

Remember, if the method is convergent then the error should get smaller as $\Delta t$ gets smaller. To visualize this trend across several scales, we use ax.loglog() to make both axes logarithmic. A relationship $E\propto\Delta t^p$ then appears as a straight line with slope $p$.

Our paper derivation predicts $p=1$, so we also plot a reference line proportional to $\Delta t$. Its vertical placement is fixed by the error on the finest grid; only its slope carries meaning. In Python, an index of -1 selects the last item, and converting dt_values to a NumPy array lets us scale all its entries at once. The equal aspect setting makes one decade occupy the same visual distance on each axis, so a slope-one line appears at 45 degrees. We are comparing the computed curve visually with the expected order, not calculating an observed value of $p$.

This is the kind of result we like to see: as $\Delta t$ shrinks toward the left, the error decreases. Over the finer step sizes, the computed curve is approximately parallel to the $O(\Delta t)$ reference line, visually supporting the first-order result derived on paper.

Refactor the code with an agent¶

Preserve the reference and define the checks first¶

Before asking the agent for code, preserve results from your direct implementation and establish independent expected values. Acceptance criteria written after seeing a proposal can drift toward whatever the proposal happens to do. Here, the original working loop provides regression evidence, while a hand-calculated state checks the equations and Euler update directly.

Add the following cell to your notebook now. The .copy() method creates a separate array whose values will not change if z is later reassigned. The test depth is 120, while the test trim depth is 80, and the test depth rate is -3. These choices make both components of the right-hand side nonzero. They also differ from the lesson parameters, so an implementation that silently hard-codes those values cannot pass accidentally.

The supplied audit will check three levels of the refactoring:

1. rhs_linear_phugoid() returns rhs_expected for the non-equilibrium state.

2. euler_step() returns u_next_expected for one step.

3. The refactored loop reproduces every value in z_direct.

Verify the two expected arrays on paper before continuing. From $u_{\text{test}}=[120,-3]$, the derivative is $[-3, 10(1-120/80)]=[-3,-5]$; one step of size 0.1 gives $[120,-3]+0.1[-3,-5]=[119.7,-3.5]$. These values come from the mathematical specification, not from agent-produced code.

These are the exact checks you will run after adding an accepted refactoring:

np.testing.assert_allclose(
    rhs_linear_phugoid(u_test, g_test, z_t_test), rhs_expected
)
np.testing.assert_allclose(
    euler_step(
        u_test, rhs_linear_phugoid, dt_test, g_test, z_t_test
    ),
    u_next_expected,
)
np.testing.assert_allclose(
    u_refactored[0], np.array([z_0, b_0])
)
assert u_refactored.shape == (num_steps + 1, 2)
np.testing.assert_allclose(u_refactored[:, 0], z_direct)

Read each assertion now and identify the mathematical or interface requirement it will check. Do not change the expected values to make a later proposal pass.

Record a compact task brief¶

Add the task brief below as a Markdown cell in your notebook. Lesson 1 displayed every specification field separately; here related decisions are combined under four headings. This shorter form is proportional to a local refactoring with supplied checks, while the function interfaces remain explicit because they are part of the numerical design.

Agent refactoring task brief¶

• Goal and scope: Propose a right-hand-side function for the linear phugoid and a separate Forward Euler step, then use both in a replacement for the single-grid time-marching loop. Put the proposal in new cells and leave the direct calculation, plots, exact solution, error calculation, and computed trajectory unchanged.

• Non-negotiables: Order the state as $u=[z,b]$, with derivative $f(u)=[b,g(1-z/z_t)]$, and calculate both components from the same old state. rhs_linear_phugoid(u, g, z_t) returns a NumPy array of shape (2,). euler_step(u, f, dt, g, z_t) returns u + dt * f(u, g, z_t) using vector operations. Store the history in u_refactored with shape (num_steps + 1, 2) and initialize u_refactored[0] = np.array([z_0, b_0]). Use NumPy and existing names only—z_direct exists, but b_direct does not. Pass g and z_t explicitly; do not use globals, SciPy, classes, in-place state updates, or a packaged solver.

• Allowed actions: Read the attached notebook and add your proposal in new cells, according to the Goal and scope above. Do not edit other cells or files, run code, use terminal tools beyond the allowed actions. Do not access the network.

• Done when: The proposal would pass every supplied check for the right-hand side, one Euler step, the initial state, history shape, and complete depth trajectory. The notebook is edited to add code cells containing exactly the two functions and replacement loop. Return at most 100 words explaining the state ordering and parameter forwarding. Propose no other changes.

This brief is the minimum sufficient specification for this task. The short prompt in the next step merely points the agent to it.

Invoke the agent, then review its proposal¶

Save your notebook so that the attached file contains the reference values, checks, and task brief. Attach it to the agent interface and send only this brief request:

Treat the response as a proposal. Before adding it to your notebook, check it against the task brief:

• Do the state ordering, signs, and right-hand side agree with Equation 6?

• Does the Euler step evaluate the complete right-hand side from the current state before updating either state variable?

• Do both functions follow their interfaces, use vector operations and explicit parameters, and avoid invented external names?

• Does the loop initialize and fill the specified two-column history while leaving the direct calculation intact?

• Did the response stay within the allowed actions, requested output, and task boundary?

Record one sentence saying whether you provisionally accept, revise, or reject the proposal and why. If you revise it, record the change. Only then add the accepted version in new cells.

Compare with a neighbor, then audit the refactoring¶

After recording your provisional decision, compare the agent proposal with a classmate’s proposal. Do not judge them by whether the lines are identical. Instead, use the task brief to identify any meaningful differences in the function interfaces, state ordering, parameter use, and time loop. If both proposals appear to satisfy the brief, explain whether the differences are merely stylistic or require additional evidence.

Add your accepted proposal to new cells in your notebook, leaving the direct implementation intact. Then run the supplied checks shown before the task brief. np.testing.assert_allclose() reports a failure when corresponding numerical values disagree beyond floating-point tolerances; silence means that the stated comparison passed. The first two checks inspect the local mathematical pieces, while the shape check and final comparison ask whether their repeated use produced the requested state history and preserved the complete depth trajectory.

If a check fails, do not change its expected value. Compare the proposed code with Equation 6, the Forward Euler update, and the specified interfaces; revise or reject the proposal only after locating the discrepancy.

Reuse the accepted functions in the refinement study¶

After the checks pass, return to the refinement study in your notebook. Keep its outer loop over dt_values, but replace the hand-written derivative and update inside its time loop with the accepted functions:

for n in range(num_steps_trial):
    u_trial = euler_step(
        u_trial, rhs_linear_phugoid, dt_trial, g, z_t
    )
    z_trial[n + 1] = u_trial[0]

Re-run the error calculation and compare it with the values saved before the refactoring:

np.testing.assert_allclose(error_values, error_values_direct)

This final regression check is stronger than deciding that the new error plot merely looks similar. It confirms that the structural change left every reported error value unchanged.

Debrief the activity¶

The agent contributed a proposed code structure, but the model equations, interfaces, constraints, and evidence were fixed before delegation. You inspected the proposal and tested it against both independent calculations and the original implementation.

What’s next?¶

This lesson marched a two-component state forward in time using Euler’s method and checked the result against an exact solution. In the next lesson, we will return to the full nonlinear phugoid model, include aerodynamic drag, and expand the state to track speed, trajectory angle, and position. The right-hand side will change, but the vector Euler step you developed here will remain the same—and, without an exact solution, refinement evidence will become even more important.