About numerics

This set requires numerical calculations. Write and submit the listing of your own codes in any language, e.g. Python scripts as .py text files or as PDF containing code and results, exported from Jupyter notebooks ( do not submit the notebooks themselves - the format is not readable in Quercus & we will not download and run your notebooks or codes!). Separate picture files are ok.

In a top comment section, sign your code, describe how you chose the value of timestep dt to assure an accurate solution. Put comments above all important lines of code, saying what they do & why. This way we I'll know it is your original code. If parts of the code are not yours, give credit to the author(s), cite the source; some points will be subtracted for that.

[15 p.] Problem 1. Precession in a funnel.

Two variables describing the particle are r(t) and φ(t). The first is the distance from the axis (not a 3-D distance from the bottom of the funnel). The third coordinate is then z = ξ r, where ξ = tan α, because the particle is always on the surface of the funnel.

Those of you who know Lagrangian mechanics from PHYB54, will be able to write the Lagrangian ℒ
[' and " are first and second derivatives over time and g = 9.81 m/s/s is the Earth's gravity constant]
ℒ (r,r',φ') = T - U = ½(1+ξ²) r'² + ½r²φ'² - g ξ r
and substituting it into Euler-Lagrange equations (click on the image, q := r or φ)

to derive two equations of motion (derive them if you want, but that is not required & there will be no extra points for it):

φ' = L/r²
r" = a^-1 (L²/r³ - gξ)

Here, L = const. is the specific angular momentum [per unit mass of the moving test particle, r · velocity component perpendicular to (r,z) plane]. L = φ' r² is conserved because the funnel is axisymmetric and so there is no azimuthal force and no torque on the particle (torque/mass = dL/dt = 0). Constant a = 1/cos²α = 1 + ξ² is a shorthand for the geometry of the funnel.

The dynamical system is conservative and frictionless. You should start the motion at t=0 with sufficiently large φ' and L but small (but non-zero) r', so that the particle circulates inside the funnel, doing small up-and-down oscillation around a purely horizontal circulation. The system has a neutrally stable circular orbit (r,φ) for every r.

Your first task is to show analytically what φ'(0) is needed to yield a circular orbit for any starting value R = r(0), and what is the value of L(R) for that circular orbit. That φ' is the azimuthal frequency (azimuthal speed in rad/s) of the strictly circular, as well as the average azimuthal speed of nearly-circular motion: ω_φ = φ'. (There is also an associated period of azimuthal circulation 2π/ω_φ, but let's just calculate ω_φ.)

Next, sketch (by hand) the acceleration r"(r) and the effective potential Φ(r), in which the autonomous r-motion with constant L takes place. Negative gradient -dΦ/dr is the r.h.s. of r-equation (acceleration).

What is the frequency ω_r of small radial oscillations? Compute ω_r² either as -1 times the slope at which acceleration r"(r) crosses the r-axis or as the curvature of the potential (∂_r)²Φ.
[Think of a harmonic oscillator r" = -ω²r, and its quadratic potential Φ=½ω²r², where ω²≡ω_r² equals both the second derivative of Φ and the negative slope of r"(r) at the fixed point r*=0 - but the relation is general! Incidentally, small oscillations around equilibrium of any system governed by Newtonian dynamics are approximated very well by harmonic oscillator, because all minima of potential locally look like a parabola]

The difference Δω = ω_φ - ω_r of the radial and azimuthal frequencies is the angular speed of precession of the orbit. If non-zero, the orbit looks like an ellipse turning in space, or rosette, i.e. it self-intersects in (r,φ) plane. If Δω < 0 then the orbit gradually rotates in the direction opposite to the sense of orbital motion, otherwise in the same direction.
Find analytically the funnel angle α for which there is no precession. In such a funnel, test particle exerting small deviations from circular motion will repeatedly pass through the starting point and its orbit is closed and periodic.

[15 p.] Problem 2. Numerical calculation of motion in a funnel

Compare that solution with the solution in case of either α = 60° or α = 45° (pick one value you like). Do you find agreement of the numerical precession speed with its analytical calculation in that case?

Write and submit your own code or Python script (see comments on numerics above). For clarity, show the trajectory seen from the top, i.e., the orbits should be shown in Cartesian (x,y) coordinates, not in polar coordinates, and not as separate functions r(t) and φ(t).

[30 p.] Problem 3. Divergence from L2 Lagrange point

To remind you: planet is at (0,0), the negative x axis points to toward the star, and we study the planar motion of a test particle in the Cartesian coordinates with rescaled length. Coordinates (x,y) are in units of Hill sphere radius (Roche lobe) radius, $r_L= (\mu/3)^{1/3} a$, where $a$ is the star-planet distance and $\mu =$ planet mass / (star+planet mass) $\ll 1$.
The beauty of Hill's equations is that they and the shapes of theor solutions apply to any small planet; the particular value of $\mu$ and $r_L$ is unimportant. It results in the same trajectory, in the rescaled coordinates (x,y). Here are the Hill's equations:
x" = -3x/r³ + 2y' + 3x
y" = -3y/r³ -2x'
where r² = x² + y².
As usual, '=d/dt is the first, while " is the second derivative over the nondimensional time t (unit of this time corresponds to time it takes the planet to turn 1 radian around the star). The meaning of the right-hand side terms has been described in the quasiblog. We already know that the positions of two fixed points are (x*,y*) = (∓1,0). We call them Lagrange points L1 and L2. One being a mirror reflection of the other, they have the same stability properties, and the same eigenvalues $\lambda$ and eigenvector directions.

A [5p]. Prove the constancy of Jacobi energy:
E_J = ½(x'²+ y'²) -3/r - (3/2)x² = const.
by explicitly evaluating the full time derivative dE_J/dt, and showing that it vanishes. Don't forget about the inner derivatives!

B [25p]. Numerical trajectory

Notice that we have four variables: x,x',y,y'. So we have a 4-D dynamical system. Program one of the many methods for solving dynamical equations. Preferably use the leapfrog scheme (looks quite like Euler scheme when the starting velocity is zero). If you use someone else's software (commercial/internet) to integrate the equations, state that and 25% of points for this section will be subtracted.

At the beginnig, set the constant timestep to $\Delta t \sim 0.05$. As you integrate the system in time, print out from time to time (printing in evey timestep creates too much output on the screen, but if you don't mind..) the value of Jacobi energy. Make sure the E_J you start with does not differ from the final value by more than $10^{-4}$, otherwise consider the trajectory inaccurate and repeat the whole calculation with a smaller timestep. State the satisfactory value of timestep and the relative E_J error you have obtained (this could be in the comments section of your code).

Compute the departure from point (1,0) toward the planet, the center of attraction at (0,0). If you place a fixed point at rest, it won't start moving. Therefore, begin from (x,y) = (0.9999,0.00005), and (x',y') = (0,0).

Continue integration of Hills equations until such time that the orbit goes around the planet about 7 times. If you extend the time, trajectory overlaps so much that it is hard to follow its course.

Plot the trajectory and overplot positions of L2 and the planet, as well as marker points of different colors, and size distinguishing them from the trajectory curve, every 0.1 time unit. Such visualization helps to see the speed in addition to direction of motion at every point of particle trajectory. Overplot the eigenvector direction (-1,0.5399) found analytically in Etudes, using a dotted line. Submit your well-commented program.

Graphics suggestions:
Try to plot the trajectory with 1:1 aspect ratio if you can, that is with equally sized units on x and y axes. In Python, one can precede plotting with command such as fig,ax = plt.subplots(), and then plot both the line of trajectory: plt.plot(trajx.trajy), and/or overlay individual points (x,y) points via plt.plt([x,x],[y,y],marker="."), and then do ax.set(aspect=1). Planet's position can be plotted like this: plt.plot([0,0],[0,0],"ro").

[25 p.] Problem 4 - Rössler system

A. Is this system a chaotic attractor?

Starting from starting point $x=y=z=1$, numerically evolve the system for $a=0.3$, $b=0.2$ and $c =6$.
Use numerical integrator of your choice (write your own script or program and enclose it as a file, cf. remarks on numerics at the beginning of assignment). Recommended method is Euler method with short time step $\Delta t = 10^{-3}$. Continue up to time $T > 50$, large enough to clarify the bahavior and illustrate the chaotic or non-chaotic behaviour clearly.
Plot the trajectory in the projection on $(x,z)$ plane (i.e. ignore $y$ while plotting). Then ignore $z$ and produce another plot with the ($x,y$) projection.

B. Bifurcation.

Vary parameter $c$ in the range $0 \lt c \lt 20$. Observe bifurcations (jumps in behaviour), describe your results and illustrate them with just a few plots. Include plots in one solutions PDF, or clearly describe and submit separately.

C. Lyapunov exponent.

Adopt $a=0.2$, $b=0.2$ and $c =5.7$. Create a variant of your code to follow two different trajectories ($x_1,y_1,z_1$) and ($x_2,y_2,z_2$) at the same time. Start the 1st trajectory from (1,1,1) point, and the second with $x = y = 1$, and $z = 1- 10^{-9}$, and plot the 3-d distance between points 1 and 2 as a function of time. Use log-linear scale, that is, plot the decimal logarithm of $D = ||(x_1,y_1,z_1)- (x_2,y_2,z_2)||$ vs. time. From the slope of the curve estimate the Lyapunov exponent $\lambda$. Desregarding little wiggles, we expect $D(t) \sim \exp(\lambda t)$, which is a straight line in a log-linear plot.