💻 Programming simple ODE solvers (like A0)

Practical numerical integration: Order of accuracy

In order to achieve higher accuracy and better speed of execution, evaluate the r.h.s. twice and average the f's at the beginning and end of subinterval dt.
  f0 = f(x,t)         ; high accuracy starting-point f
  x_end = x + dt * f0   ; low accuracy (1st order) prediction of endpoint x
  f1 = f(x_end,t+dt)   ; quite high accuracy endpoint value of f
  x = x + dt*(f0 + f1)/2   ; high-accuracy (2nd order) prediction of endpoint x

The above is a 2nd order accurate integrator, one of many. That means the error is of order O(dt²), i.e. faster convergence with diminishing dt than 1st order scheme can provide. It's called the trapezoid method.

There is also a midpoint scheme, where first r.h.s. evaluation is used to get an accurate x-midpoint at time t+dt/2, and then r.h.s. is re-evaluated to get an accurate estimate of f(x_mid,t+dt/2), i.e. the average f over time from t to t+dt (this f is used for the final update of x). Midpoint scheme is also of 2nd order accuracy.

The diagram shows the accuracy as a function of dt (on log-log plot) for the trapezoid and the 1st order Euler scheme. Click to enlarge. Make sure you understand the slopes of the two curves, representing what we call a truncation error (in the sense of truncation of Taylor series. It's simply an error caused by limited accuracy numerical scheme and the necessary discretization of time.)

The reason the blue scaling breaks down and the red one levels off at very small dt's is the roundoff error in computer double-precision arithmetic. The more timesteps you take, the further you random-walk away from the true solution.

Finally, try this crazy idea and see if you can save some computational time by doing only 1 evaluation of r.h.s. per subinterval of time (dt) -- the evaluation is usually the most expensive part of the algorithm -- and still get 2nd order accuracy!
 
  x_mid = x + (dt/2)*f   ; use f from previous time step, low accuracy
  f = f(x_mid,t+dt/2)   ; quite high accuracy midpoint f
  x = x + dt * f   ; high-accuracy prediction  

No, you won't succeed... :-( This is still a 1st order scheme. If you could remember more old values of force (a multistep scheme), you could make a better prediction of midpoint x by effectively estimating the rate of change of the r.h.s., and have a higher order scheme while doing only 1 new evaluation of the r.h.s. per dt.

Amazingly, something like this unsuccessful trick does work for 2nd order ODEs like Newton's momentum equation dv/dt = f(x,t).

The simplest 1st order integration scheme is to turn one 2nd order ODE into two 1st order ODEs.
  v = v + dt f(x,t)
  x = x + dt v

In this form, if we take care to start (and later end) the computation in [0,T] time interval with de-sychronizing operation v = v + (dt/2)*f(x,t) ONLY in the first step & very last timestep (with no x-update), then our v's will be evaluated at midpoints of the time intervals [t,t+dt], while x's at the ends only. Positions and velocities will be leapfroging over each other, as we proceed through time. The scheme is called leapfrog integrator and is 2nd order accurate. It is even symplectic (preserves energy very accurately with no systematic drift).

If however you swap the order of the v and x updates in the above scheme, the accuracy is only 1st order. So, just in case, remember to always update velocities first. Even if you forget to desynchronize variables at the beginning, along the way your scheme will have 2nd order accuracy.

✈ Theory of flight. Phugoidal oscillations in times of N. Zhukovsky , F. Lanchester and today.

Equations of motion and fixed points (steady flight)

(v*)² = g/L, or Lv* = g/v*,
θ* = -D/L.

Notice that our D/L stands for the ratio of coefficients of drag and lift, but more generally it is the ratio of forces (proportionality of both forces to v² is not required.) A steady-state glide happens at θ* = -D/L ~ -1/30 rad (-2^o) for good gliders, and less than 1.5^o for the winning gliders. Polish glider called Diana 2 has a fantastic max glide ratio 50:1, meaning -1/50 rad = -1.15^o glide. For comparison, a typical older airliner (like TU-154 or Boeing 767) will glide some 17 times further than it sinks, at an angle 3.4^o to the ground. For that reason, standard final glide paths to airports tend to have a prescribed slope of at least 3^o, typically 3.5^o, so that a jetliner with engine malfunction has a chance to retract flaps, fly at optimium speed, and land on a runway faster than usual but in one piece. Flaps are needed for nice and slow landing. They increase drag by a larger factor than lift, which spoils the aircraft's aerodynamics (here, L/D); on a standard approach the jetliner must thus keep a noticeable fraction of engine thrust to maintain the prescribed shallow glideslope. If you live near the approach path to an airport, you may have noticed it.

More modern designs like Boeing 787 Dreamliner advertise glide ratio 21:1 (-2.73^o), which happens to be close to the performance of an excellent soaring bird albatross, who has L/D ≈ 22, so θ* ≈ -1/22 rad or -2.6^o. Albatross can soar 1000 km a day. Its flights take it around the globe. To reach such high glide ratios L/D, wings must have a very elongated shape. That's why both gliders and albatrosses have such wings. The wingspan of the latter reaches 3.4 m. Also, to keep wings fully extended albatrosses have special lock on tendons and do not use muscles to keep wings steady; they use ascending air currents and winds, and do not waste energy on flapping their wings either, unlike sparrows, who have a horribly bad L/D = 4 and sink 1 m per every 4 meters covered in gliding. OK, from birds, gliding frogs and flying fish, back to our equations, via a brief interlude on piloting and landing an airplane.

Paradoxes in pilotage

Look at the θ* equation: θ* = -D/L. You want to reduce |θ*|, so you need to reduce the D/L ratio. Lv²=W in stable flight (lift = weight). The only way to reduce D/L= (Dv²)/W is to add engine thrust force or rather, acceleration T, to reduce Dv² (acceleration of parasitic drag) by T. In no case should you re-trim the airplane or pull on the controls to raise the nose, while leaving the thrust unchanged. That will only INCREASE the total drag (by raising the so-called induced drag), lower the speed, and reduce your range. Even though the nose will be pointing higher (say, in the direction of the runway), fooling you to believe that you're now doing better, you will actually sink faster and crash short of the runway. Which will spoil your whole day.

Another issue is the speed. Let's say you glide and suddenly notice on the airspeed indicator that v=v* is dangerously low, too close to the stall speed of the aircraft, below which the smooth airflow is guaranteed to detach from most of wing's surface (L drops, D goes up dramatically, and the wing is more falling than flying). You need to correct v* up to the recommended, safe, value equal 1.3 times the stall speed. From driving a car you have an intuitive reaction to step on the gas, in order to speed up. Very wrong again! If you do increase engine power, the thing you change is θ* (which we just covered), but the speed v* won't go up.
Since v* = (g/L)^1/2, what you need to do is to lower the lift coefficient L, which depends on the angle of attack rougly linearly. You do that by pointing the nose down, which reduces the angle of attack, the lift coefficient L (and part of drag called induced drag). Congratulations, you saved the day avoiding a stall-and-spin accident!

There is a lot of nonlinear physics in flying. For an example of aerobatic maneuvers such as rolls, loops and an immelman (half-loop + half-roll) see a video with your prof. instructing a friend in an aerobatic-capable plane RV6A. To find out more, check the basic and advanced aerobatics textbooks.

Stability of the fixed point

v=v* (1+ε), and θ = θ* + φ,
where ε and φ are small, dimensionless deviations of speed and angle (angle is in radians).
After linearization, using g/v* = Lv*, we obtain \begin{gather} \frac{d}{dt} \begin{pmatrix} \epsilon \\ \theta \end{pmatrix} = \begin{pmatrix} -2Dv^* & -Lv^* \\ 2Lv^* & -D v^* \end{pmatrix} \; \begin{pmatrix} \epsilon \\ \theta \end{pmatrix} \end{gather} In shorter, matrix form, dx/dt = A x, where x = (ε,φ)^T (transposed or vertically written vector) and A is the 2x2 matrix shown above.

The trace of matrix A is negative,
τ = tr A = -3Dv* < 0,
whereas the determinant is always positive:
Δ = det A = 2(Lv*)²+2(Dv*)² > 0.
Considering that if D/L ≪ 1, then τ²-4Δ < 0, we have now established that the glider is in a stable equilibrium (cf. the classification of linear nodes on p. 138 in the Strogatz textbook).
If you don't believe that, convince yourself step by step, writing down the eigenproblem fully: A x = λ x, or (A-λ I)x = 0, where I is a unit 2x2 matrix. Non-trivial solutions exist only if det(A - λ I) = 0, which we write out as the characteristic equation
λ² - τ λ + Δ = 0
that yields two eigenvalues λ_1,2. Writing this quadratic function of λ as (λ-λ₁) (λ-λ₂)= λ² - (λ₁+λ₂) λ + λ₁λ₂, and comparing that from with characteristic equation we see that the sum of eigenvalues equals τ and is negative, while their product Δ = det A is positive. But λ₁ and λ₂ are not simply two negative real numbers. Inequality τ²-4Δ < 0 implies that they are a pair of complex-valued eigenvalues, because the explicit solutions of the quadratic equation for λ's are
$$ \lambda_{1,2} = \frac{\tau}{2} \pm \frac{i}{2} \sqrt{4\Delta - \tau^2}$$ The expression under square root is positive. The imaginary parts have opposite signs; we can construct sinusoidal oscillations from them. Most importantly, the real parts of eigenvalues are equal and negative: τ/2 < 0. This tells us that both independent solutions are damped in time (how fast, we will see below). The equilibrium point of the 2-D dynamical system is classified as stable spiral, and the flight of an airplane or bird is self-stabilizing.

As will become evident soon, we have found a fairly slow oscillation of an aircraft in flight, involving only pitch changes but no bank, both features distinguishing this mode from a faster, longitudinal+lateral oscillation called Dutch roll. It is known in aerodynamics as phugoid or phugoidal oscillation. We will further analyze phugoid oscillations below, after a historical interlude.

Historical notes

Lanchester's early aerodynamical paper "The soaring of birds and the possibilities of mechanical flight" was written in 1897, much before the first airplane was built. He was ahead of his time. [In 1927 he also designed a hybrid, "petrol-electric" passenger car; trucks used hybrid engines even earlier! One of such trucks pulled the telescopes up a dangerous path to Mt. Wilson observatory from Pasadena, CA, in 1908-1918. But it proved too weak, and was supplemented by an auxilliary donkey in front.] In the second part of a two-volume book on flight published in 1907 and 1908, both on microfilm in U of T library, Lancaster described the phugoid (minus the damping). The phenomenon is described in wikipedia, although in scant physical and mathematical detail.

There was a genius aerodynamicist who first explained the lift force as effect of circulation of air around airfoil also studied the theory of phugoid, earlier than Lanchester. He also predicted the possibility of making loops before an airplane was constructed, in a paper on gliding of birds. That man was Никола́й Его́рович Жуко́вский or Nikolai Zhukovsky (sometimes spelled Joukowsky).

I have recently dug up some works earlier and more comprehensive than Lanchester's by Zhukovsky (Жуковский Н.Е. О парении птиц. Мoscow, 1891), cf. this 1922 publication. It turns out that it was Zhukovsky who first formulated the differential equations that we set up at the beginning of my story of phugoid oscillation. He solved them by an approximate method, including the drag force, which Lanchester ignored (denoted F on the first picture):
  Here is more on airfoil theorem and the Joukowsky transform and airfoils.

By the way, Russian Empire's military pilot Пётр Николаевич Нестеров (Pyotr Nikolayevich Nesterov) after convincing himself by calculation that such a feat is possible to accomplish, performed in 1913 at an airport near Kiev the first loop (then called a dead loop) in a French-designed, Russian-built, 70 horse-power Nieuport IV airplane. Initially arrested for unauthorized action and endangerment of Army's property, Nesterov was soon released and given on order of merit.

Back to our eigenproblem

We already obtained two eigenvalues, which we can evaluate by substitution of τ and Δ: $$\lambda =-\frac{3Dv^*}{2} \pm i\sqrt{\frac{2g^2}{v^{*2}}-\frac{9D^2v^{*2}}{2}}.$$ In the square bracket, the second term is ~(L/D)² times smaller than 2(g/v*)² = 2(Lv*)² and can be neglected. We can write the eigenvalue as $\lambda = -t_d^{-1} \pm i \omega$, where $$\omega = \sqrt{2} \frac{g}{v^*}$$ is the oscillation frequency, and $$ t_d^{-1} = \frac{3Dv^*}{2} = \frac{3}{2^{3/2}}\,\frac{D}{L}\;\omega $$
is the inverse time of amplitude damping, both constants having the physical unit 1/s. The time scale t_d is the e-folding amplitude decrease time (decrease by factor e≈2.718). $$P = 2π/ω = \sqrt{2}\, \pi \frac{v^*}{g}$$ is the period of phugoidal oscillations.
Comparing t_d and P, we get $$ t_d = 2/(3 D v*) = (L/D) \sqrt{2} (3\pi)^{-1} P \simeq 0.150 (L/D) P.$$
Phugoid damping in good gliders with L/D > 30 takes t_d > 4.5 P, that is more than about 4.5 periods of oscillation. Badly designed aircraft are actually damping phugoids better, not worse! In them L/D is lower and damping takes fewer periods. But weak damping of phugoids is a small price to pay for small drag and high gliding ratio, Besides, they can be damped by dedicated damping systems deflecting control surfaces.

Very interestingly, phugoid period P does not depend on anything except the forward speed of flight: neither on how aerodynamic the aircraft's design is (what's its L/D ratio), nor how hot, dense, humid or viscous the air on a given day is, which depends on the atmospheric conditions and flight altitude, etc. The airship does not even have to be a glider, it could in fact be an airplane powered by engine(s). Could you redo the derivation with a constant thrust +f added along the path? This would surely change the equilibrium glide/ascent angle but not the formulae for oscillation and damping periods, in our simple theory. Would thrust change v* much? Do the calculation - you may be surprised and learn a bit about piloting. On approach, pilots most aircraft control glide angle but not the speed with engine thrust. They control speed by changing the AOA of wings with the elevator, a control surface attached to horizontal stabilizer. The stabilizer itself is responsible for short-term longitudinal stability (also damped oscillations), but the time scales of these are an order of magnitude shorter than what we discuss here.

In practical terms, our phugoid formula predicts $$ P \approx \frac{v}{100 {\rm m/s}} \; 45.3 \;{\rm s}.$$

For airplane flying at 306 km/h = 85 m/s, such as during the approach of flight number PLF 101 of a Polish presidential TU-154M aircraft on 10.04.10 that exhibited large vertical speed and glidepath variations resembling a phugoid (and subsequently collided with trees and crashed killing all 96 people on board), we obtain a phugoid period P ≈ 39 s, roughly in agreement with the black box recordings. The phugoid phase was unlucky in that the nose started dropping slightly at a moment when the trajectory should have been corrected the the opposite direction, but phugoid was not a cause of the accident, merely one of adverse circumstances. The crash was a consequence of erroneous pilotage in difficult meteorological conditions in Smolensk, Russia (low cloud base and dense fog). The crew was not supposed to descend below the MDH = minimum descent height, 120 m above the runway, because they had zero visibility there. They continued a too-steep descent until they emerged from the fog and saw the ground just 25 m below them. It was too late to avoid collisions with trees, resulting in loss of a 1/3 span of the left wing, an uncontrolled roll to the left and an inverted crash into a sparse forest. I was analyzing that catastrophe as an independent researcher, supporting the Polish governmental committee. Below is my reconstruction of vertical velocity and altitude deviations in the approach path, which should be smooth and much less steep (or else a go-around must be executed, which wasn't done). The vertical speed and approach cone indicated in green are the mandated ranges in a non-precision approach in TU-154M at Smolensk Severnyi airport (click the images)

 
Time in seconds is on the horizontal axis in the plot of v_z on the left, and the distance to runway threshold in km is on the horizontal axis in the figure on the right (altitude above the runway in meters is shown on the vertical axis).

Another example is a flight at cruise speed of 500 mph = 435 kt = 805 km/h = 224 m/s (see Etkin, book "Dynamics of Flight" 2005).
Etkin employs a much fuller set of dynamical equations, actually taking into account those things we said are not entering our equations, including all the so-called aerodynamic derivatives and aircraft's moments of inertia of an airliner. He gets the following answer: phugoid of period 115 s, e-folding damping over approximately 3 oscillation periods. In addition, the richer set of equations predicts more eigenvalues and a very quickly damped short-period oscillation, which we do not predict at all. We get the period P(500 mph) = 101.4 s. That is a 12% shorter P than from Etkin's calculation.

Incidentally, the 2005 edition is based on a 1972 edition, where the eigenproblem was solved on an computer in those long-gone times when computers were so rare that their use merited a mention in the book; it was an IBM 1130 mainframe computer at UTIAS (that's our UofT Institute for Aerospace Studies. Etkin was UofT Prof. and Dean. He died 10 years ago.)

OK, so we underpredicted the computed phugoidal period (Etkin 2005) shown here,

and also this actual flight test measurements at v* = 150 kt, where a small jet was oscillating in about P = 50 s, while we computed P = 35 s. Nevertheless, our result is still classical, as it gives the phugoid period exactly the same as in both the original derivation by Lanchester (1908), who used energy conservation principle and the assumption that engine thrust exactly cancels the drag, and in a lot of subsequent derivations making fewer assumptions than we made, as discussed in aerodynamics books such as Cook's "Flight Dynamics Principles" (2nd ed. 2007), from which we took the right-hand figure above.

If you want to learn about other types of airplane oscillations and instabilities, they are discussed e.g. here. Phugoids are sometimes shown in videos like this one, but they are so boringly slow that I guess they actually should be shown at 4x fast-forward speed.

Why are phugoids important? Because as already mentioned, unlike the short-period oscillations, the long climb and long descent in a phugoid causes large altitude deviation. A pilot who needs to maintain level flight has one more thing to actively correct, unless the plane is in cruise and is steered by multi-axis autopilot enforcing the longitudinal stability. However, on the final approach autopilots in most airplanes must by disengaged, and in the Smolensk crash it was not controlling the glidepath. Captain Sullenberger, who in January 2009 successfully ditched an Airbus 320 on the Hudson river in NYC, following two bird strikes and a complete loss of engine power, said that the phugoid-control circuit was interfering with his flare (increase of θ just before ditching the plane). This sounds a bit strange, but who knows. There are no other cases to compare.

🐕 A dog, 🦆 a duck, and the mystery of an invisible magnet

SOLUTION

There are different kinds of solutions. One of the previous years, half of students voted for a numerical solution first, and half for analytical solution first. Both are good choices for the first stab, as long as you attempt the other method as well.

As a matter of fact, if you start with numerical solution then you're immediately faced with the problem of what equation to ask the computer to solve. Let's talk about that first, then. I will use a compact kind of notation, using complex numbers in place of vectors, although at the cost of writing twice as many equations one can work directly with (x,y) or with (r,φ) coordinates of the dog in the program. Complex numbers just look.. less complex in a program. You may think of them as vectors, which they are. Let the position of the dog be r = x+iy, and D be another complex number giving duck's position.

The dog goes straight for the duck, so along the vector D-r, at speed s, i.e.
  dr/dt = s (D-r)/|D-r|
where the ducks position is known to uniformly move in a unit circle:
 D(t) = exp(it).
That's not the end, since both animals are moving in some circular or quasi-spiral way. Rather than talk about position r(t), we'd better switch to vector R,
 R = exp(-it) (D-r)
which is an instantaneous vector from the dog to its moving target, in a reference frame that rotates at unifom unit speed with the duck. In such a rotating frame, the duck is always at (0,0) position and R won't run wildly in quasi-circles, obscuring the plot. The required frame rotation is accomplished by turning any vector to the right by t radians, which is easily done in complex plane with multiplication by exp(-it).
We work out, by differentiation of the definition, the derivative of R to be
 dR/dt = -i exp(-it) (D-r) + exp(-it) [dD/dt - dr/dt].
Taking into account that dD/dt = i exp(it), and that dr/dt=sR/|R|, we obtain our evolution equation for a 2d dynamical system in the complex form

 dR/dt = s R/|R| + i(1-R)

The numerial solutions for R(t), t=0...500, were obtained with this second-order integrator (dt = 0.004 was sufficient for accuracy), for a series of 20 uniformly spaced values from 0 to 1 (s = 0.05, 0.1, 0.15, ...).

The higher the dog's speed s, the closer to the duck it eventually gets. For s=1, it reaches the duck in infinite time (and for s>1 in finite time). In every case, there is a fixed point to which the solution tends over the time of integration. Two green lines in the plot are circles drawn around the center of the pond (1,0), and have radii |R-1|=0.5 and 0.8 of the pond radius, for illustration of an important fact: the position of a (moving) dog with respect to (moving) duck tends to the radius equal s, as long as s is not larger than 1. You can see that two of the plotted trajectories (for s=0.5 and s=0.8) end at radius s.

We didn't have to do numerical solutions to find this fact, or even the precise asymptotic positions. Fixed positions can be found analytically from the evolution equations, as usual. Zeroing the r.h.s. we obtain R* (we'll skip the * for simplicity, to avoid confusion with dot product or multiplication):
 i(R-1) = s R/|R|.
Taking absolute value of this equation, we get |R-1| = s, precisely the point we have noticed empirically; R-1 is the final position of the dog in the rotating frame (dog enters a limit cycle that is a circle of radius s in nonrotating plane). Now we've proven it. Notice that we can also prove it ad absurdum: If |R-1| were not equal to s, the dog in its fixed position relative to the duck would travel at a linear speed differing from s (their angular speeds in fixed configuration being the same, the dog's linear speed is equal to its radius |R-1|, but that number would be different from s). This is a contradiction.
Furthermore, (R-1) as a vector is perpendicular to R, because the complex numbers representing them differ by a factor i/s, which represents rescaling by (1/s) and rotation by 90 degrees.

What is the locus of points such that vectors drawn from (0,0) and (1,0) to those points are at right angle? It's a circle whose diameter joins the two points, according to an ancient Greek geometry theorem. We can prove it by writing the usual equation for the sought-after curve in terms of coordinates x and y. We have R =(x,y), so the dot product of the two vectors R-1 and R must vanish: (x-1,y)*(x,y) = 0. This can be written as either x(x-1) +y² = 0 or more familiar form of the equation of circle (x - 1/2)² + y^2 = (1/2)². This is such a valuable result that we color half of this circle gold in our plots. So the fixed points are at intersection of two circles with radii 0.5 and s, trailing the motion of the duck by an angle θ = arccos(s). The family of solutions makes a beautiful shell:

A final comment should be made about a visible change of behavior of those solutions which have a chance to approach the duck closely, i.e. for 0.9 < s < 1. The dog rushes towards the duck (lagging behind its circular motion by about 30 degrees), then slows down its approach to the duck's unit circle and turns sharply near the point M=(0.0504,0.2188) to either approach the duck very nearly along the golden path or alternatively turns RIGHT and swims toward its fixed point described above. The origin of the invisible magnet at M is totally enigmatic to your instructor, but then he hasn't read any of the historical literature on this problem cited by Strogatz. If anybody knows or can find out why the trajectories have a bifurcation near M, please let us all know during tutorial. The complicated shape of trajectories certainly prevented 19. century mathematicians from finding solutions to the problem in terms of elementary functions (they do not exist). However, a mystery remains. The fixed points are quite smoothly distributed along the golden bow, with density decreasing towards its left end for obvious reasons (according to arccos function describing the angle betweeen the duck, the dog, and the fixed point). They do not have any bifurcations or concentration points. They're oblivious to M, while the trajectories, especially the ones with dog's speed s≈0.97447, clearly are not and are ruled by something else.

Perturbation equations for R3B problem (Hill's approximation) with radiation pressure.

x" = -3x/r³ +2 dy/dt +3x +b
(+b stands for a constant acceleration vector b along x-axis);

y" = -3y/r³ -2 dx/dt
where " is second derivative over time, b is rescaled β, the ratio of forces of external radiation pressure and gravity of the object located at the origin of coordinate system. Around the Earth, where radiation is from the sun, b = βa/r_L ~ 100 β = const. (nondimensional radiation pressure to sun's gravity ratio).

First order perturbation theory can be applied to an initially slightly elliptic orbit (e≪1) of a satellite circling around a gravity center at (0,0) (the Earth), subject to a small acceleration b pointing along the x-axis (representing radiation pressure from the sun). Satellite's orbit lies in (x,y) plane and has instantaneous orbital elements a,e,w (a=semi-major axis, e = eccentricity, w = pericenter longitude, here defined as the angle between the negative y-axis and the direction to pericenter; variable w(t) should really be ω(t), but w looks similar and is easier to type).

Carl Friedrich Gauss at the end of the 18th century created a perturbation theory, from which today we can obtain the following evolution equations (2-d dynamical system) after time averaging over one period of the unperturbed orbit (in 1st order perturbation theory we evaluate the perturbation along an unperturbed, known trajectory, here - ellipse with semi-diameter a):

da/dt = 0

dw/dt = -β n (1-e²)^1/2 sin(w)/e - n_p

de/dt = +β n (1-e²)^1/2 cos(w)

where n, n_p, b = const.

[Explanation: The perturbing force is given by nondimensional number β, the mean angular speed of the satellite is n, and n_p is a precession speed of the initial orbit due to frame rotation, if any.
A small rate of precession (n_p ≪ n) of satellite orbit is due to the Coriolis force in (x,y) frame if it rotates. Why would it rotate? Imagine that the (x,y) frame is attached to a planet that slowly orbits around a distant sun, while the x-axis always points away from the sun (so it rotates in inertial frame). The n_p is then the negative of angular speed of the planet.]

(i) If n_p is zero because the frame (x,y) is not rotating, solve exactly the evolution of w and e as a function of time t, starting from e₀ = e(0) ≪ 1 and w₀ = w(0) (arbitrary) at t=0.
Hint: Change variables to X = e cos(w), Y = e sin(w). Vector e=(X,Y) is called Laplace-Runge-Lenz vector in celestial mechanics. Draw its evolution in the orbital plane.

(ii) Additional question:

Explore the limit of e₀ → 0. The longitude of pericenter equation contains a term ~1/e, which is very large initially; if e ≪ 1 then the evolution of w is much faster than the evolution of e. Show that the system then becomes 1-D, and w adjusts to final value instantaneously. Write the expressions for w(t) and e(t)
The diameter of the orbit (equal to 2a) is constant, but eccentricity can grow to 1, which represents a needle-like orbit. Demonstrate that the orbit becomes needle-like (e=1) and must collide with the planet within a finite time; write an expression for the maximum survival time.

Try to solve the problem if n_p=const. is negative and small compared with n, but not zero. Is it possible to solve the Gauss perturbation equations analytically in that case?

Stability of Lagrange points in R3B problem

The origin of terms in the equations is given on the blackboard (gravitational forces of two massive bodies, centrifugal and Coriolis forces). Click on the pictures:

On the far right of the blackboard, you also see the typical way we write the gravity acceleration from a planet located at (0,0), using the usual dimensional distance r:   -GM_p r/r³. It is simply the absolute value of Newtonian acceleration GM_p/r² times the versor pointing toward the center of the planet, -r/r. After division by r_L, the planet's gravity looks like -3r/r³ in the new nondimensional coordinate, confusingly also denoted by r. Number 3 comes from the cube of r_L in the denominator (cf. the definition of r_L above).

A native Sardinian, Giuseppe Luigi Lagrangia (better known as Joseph-Luis Lagrange, though he moved from Berlin to Paris aged 51 in 1787, so rather late in life), has found that the R3B problem has 5 equilibria/fixed points, places where a massless test particle can rest forever: 3 collinear with star and planet (L1-L3) and two in exact equilateral configuration with the massive bodies (L4, L5). Hills equations represent two of them, L1 and L2, while the other are too far away to be captured by the approximate form of equations. The scientifically important question is the stability of Lagrange points. Let's find out the distance to L2 (the same as to L1) from the planet's center, and then analyze the stability of the L-points.

As you see in the snaphots, L-points are found at (x^*,y^*) = (±1,0). That is, L-points are indeed at the physical distance of r_L from the planet. Introducing small perturbations via:
(x,y) = (x^*,y^*) + (ξ,η),
we linearize the Hills equations, dropping terms of order O(ξ²,η²) and higher, since by assumption they are exceedingly small (even ξ,η ≪ 1). You can see that near the L-point the y-component of restoring force is -3η but x-component has a repelling, positive coefficient +9 in front of ξ. This is because L2, in common with other collinear points, are saddles of the force field potential.

We seek solutions in the form (ξ,η) = (A,B) exp(iωt), where the frequency of oscillations ω in general can be a complex. If Im(ω) < 0 then the initially small deviations from the L-point grow exponentially in time and the point is dynamically unstable. Conversly, Im(ω) > 0 means Lyapunov stability (slutions stay close to the L-point forever).

I don't have a snapshot of the standard matrix analysis of the linearized equations. You will finish this job in assignments set 3, problem 3. One ω² will be positive and the other negative. The existence of an imaginary frequency ω following from that means that L1 and L2 points are unconditionally unstable. The test particle will depart from these points upon the slightest perturbation. Despite that, it is precisely L2 where NASA and ESA have put their robotic observatories such as JWST and Planck, and sun-observing platform SOHO is near the L1 point! Why are these agencies not placing the spacecraft in the L4-L5 points, which are conditionally stable? (μ < 1/27 is neded for that, a requirement met by sun-Earth, and a lot of other binary systems).

The reason is that the triangular Lagrange points are maxima of effective potential. Too much fuel would have to be spent to climb the potential, in order to place objects there. (Incidentally, despite being maxima, these points are usually stable on account of strong Coriolis force, caused by frame rotation. This fictitious force bends the trajectories that try to roll down the effective potential and depart from the L4/L5 points, back toward them.) Meanwhile, L1 and L2 points are not only 100 times closer to Earth than L4 or L5, which is good for communication, their instability is so weak that it is easily corrected by occasional activation of thrusters on the spacecraft. It's one of the paradoxes of celestial mechanics: stable points are useless, while unstable points are useful.

Solve s" = s(1-s²) -v^q numerically

The second method we explored was perturbation theory. We noted that if we start very near the fixed point, the amplitude of oscillations is small, and the amplitude of v^q is small or even much smaller, if q > 1. Thus we can treat the power of v as a small addition to the periodic solution of s" = s(1-s²) equation linearized around the (s*,v*)=(1,0) fixed point, which reads x" = -2x, when we denote x=s-s*. The linear solution is sinusoidal function of time with angular frequency ω = √2. Linear equation has energy E = ½(ds/dt)² +½ωx = const. To find out the addition to this energy in time equal one period of the orbit, we integrate from 0 to full time period the product of the perturbation to acceleration (or force) multiplied by v, here -v^q+1. This is so because force component along the velocity vector times (dot product) velocity vector in physics is power, or the rate of change of energy E. Since v does small harmonic oscillations around zero and for even q the term v^q+1 is an odd function of time, the integral expressing ΔE is zero and the energy after a full circle in phase space is unchanged. Physically speaking, on half of a circle perturbation adds energy, and on the rest it removes an equal amount.

Not so if q is odd: the power of perturbation, (-v^q)v, stays negative and steadily saps energy out of the system's oscillations, making them die exponentially fast. We see this in the decay of amplitude in numerical results above.

Numerical investigation of chaos in Lorenz dynamical system

Scripts from tutorial that you can freely use:
You can for instance use the ideas from the scripts and re-write them in Python or Matlab or C. The IDL script is a template but not a final program that you need, so don't copy it mindlessly - it might not do what you want to do, have inefficiencies etc. The computational kernel in it was verified though, and that's what you want to use.

★   Chaos in Hills equations

In tutorial 3 we study the integration schemes of 2nd order Euler-like, and Leapfrog (Verlet) method. I used IDL (interactive data language) scripts that also run under GDL (free IDL), and you are encouraged to rewrite them in Python, Mathematica, Maple, Matlab or Octave (free Matlab replacement), Fortran 77 or Fortran 90/95/2003 (great language for number crunching), Julia, Haskell, or C/C++ . The re-writing is a good idea, as the skills you gain while doing it will be very useful for you outside this course.
inte1.pro
inte2.pro
inte2.png (results of inte2.pro)
Later, you can study this program: hill.pro
Examples of results are here:

Program hill.pro which was used in the tutorial. Solves Hill's problem for a range of starting positions outside the orbit of the planet (positive radial coordinate x). Vertical motions are neglected, but you can add the z-equation to the system, the only z-force is from the planet, so it will look quite similar to the x and y forces.
We start by identifying larger scale features x-y view of orbits (units are Hill sphere radii) , and the corresponding graph of
function x(x0) . Not clear if it's chaos. To establish chaos strictly, we'd need to evaluate and plot the Lyapunov exponents at different x0's.
Or one can start zooming in ... but there does not seem to be a self-similarity in the graphs, it's complicated system that may not have it, even though arbitrarily close starting conditions may diverge into different outcomes, and in that sense the outcome is likely a fractal object.