Homework 4, March 2025¶

1. Temperature profile of the Earth and Spherical Harmonics¶

Consider a model for temperature variation on the Earth Surface $$T(\theta,\phi)=−20+35 \sin(\theta)+5 \sin(\theta)\cos(\phi).$$ Here $\theta=0$ corresponds to the North Pole and $\theta=\pi/2$ to the equator. The first term, $−20+35\sin(\theta)$, represents the dominant latitudinal dependence. The additional term $\sin(\theta)\cos(\phi)$ produces a modest variation in $\phi$, which is zero at the poles and peaks at the equator.

Tasks:¶

1. 3D visualization: Plot the temperature profile of the Earth $T(\theta,\phi)$ in 3D using plot_surface or similar function.

2. Spherical Harmonics expansion: Expand $T(\theta,\phi)$ in terms of spherical harmonics up to $l<=20$ and and determine the expansion coefficients. Evaluate numerically coefficients $a_{lm}$ defined below for $l=0$, $l=1$, ...$l=20$, until you reproduce original function to a reasonable accuracy.

3. Numerical verification: Verify that the expansion accurately reproduces the original function $T(\theta,\phi)$. To do this, plot both the approximate and exact temperature profiles as a function of $\theta$ at $\phi=0$ and $\phi=\pi/2$ as a 2D plot.

Once you have values for $a_{lm}$, use the finite sum formula to reproduce the temperature profile $T'(\theta,\phi)$ and verify that the finite sum is close enough to the original function $T(\theta,\phi)$.

Useful formulas:

Since Spherical Harmonics form a complete orthogonal basis, any function of $\theta,\phi$ can be expanded as a series: \begin{eqnarray} T(\theta,\phi) = \sum_{l=0}^{l=\infty}\sum_{m=-l}^l a_{lm}Y_{lm}(\theta,\phi) \end{eqnarray} To compute coefficients $a_{lm}$ we multiply both sides by $Y^*_{lm}$ and integrate over the full angular domain. Using othogonality of $Y_{lm}$, we obtain \begin{eqnarray} a_{lm} = \int_0^{2\pi} d\phi \int_0^\pi d\theta \, \sin(\theta)\; T(\theta,\phi) Y^*_{lm}(\theta,\phi) \end{eqnarray}

2. Double pendulum¶

  1. Energy Conservation Analysis: Implement the total energy $E$ and analyze how well energy is conserved over time. To do this, plot $E(t)-E(0)$. The precision of energy conservation depends strongly on the choice of initial conditions. Select the following initial conditions $[\theta_1^0,\theta_2^0]$ to study their impact:

    • $[\pi/4,\pi/4]$
    • $[\pi/4,0]$
    • $[\pi/2,\pi/4]$
    • $[\pi/2,\pi/2]$
    • $[\pi/2,0]$
    • $[\pi,0]$

    Assume that $\dot{\theta}^0_1=0$ and $\dot{\theta}^0_2=0$.

    Use $N=1000$ time points, evenly spaced in the interval $[0,100]$ seconds.

    The total energy of the system is given by the Hamiltonian

\begin{eqnarray} H[\theta_1,\theta_2,p_{\theta_1},p_{\theta_2}] = \frac{2}{ml^2} \frac{3 p_{\theta_1}^2+12 p_{\theta_2}^2-9\cos(\theta_1-\theta_2) p_{\theta_1} p_{\theta_2}}{16-9 \cos^2(\theta_1-\theta_2)}-\frac{1}{2} m g l (3\cos\theta_1+\cos\theta_2) \end{eqnarray}

which, when evaluated with the computed solution, gives total energy $E(t)$, i.e., $E(t)=H[\theta_1(t),\theta_2(t),p_{\theta_1}(t),p_{\theta_2}(t)]$.

  1. Optional: Can you derive Hamiltonian from Lagrangian given below? The procedure follows the definition:

$H = \sum_i \dot{\theta}_i p_{\theta_i} - L$

  1. Sensitivity to Initial Conditions and Chaos: The double pendulum exhibits chaotic behavior for large amplitudes. To investigate this, analyze how small changes in initial conditions grow over time.

    • Plot how far are two solutions as a function of time using the distance function $$d(t) = \sqrt{(\theta_1-\theta_1')^2+(\theta_2-\theta_2')^2}.$$ Here $\theta_1$ and $\theta_2$ are solutions using initial conditions from problem (1), while $\theta_1'$ and $\theta_2'$ use very similar initial conditions, expect $\theta^0_1$ is slightly perturbed ($\theta^0_1\rightarrow \theta^0_1+10^{-6}$).

    • Compute the average distance for long times $$d_a = \frac{2}{T} \int_{T/2}^T dt d(t)$$ and estimate Lyapunov exponent $\lambda$, which is defined by $d_a = \exp(T \lambda)$.

Useful formulas:

Recall the equations on motion for double pendulum:

${\dot \theta_1} = \frac{6}{m\ell^2} \frac{ 2 p_{\theta_1} - 3 \cos(\theta_1-\theta_2) p_{\theta_2}}{16 - 9 \cos^2(\theta_1-\theta_2)}$

${\dot \theta_2} = \frac{6}{m\ell^2} \frac{ 8 p_{\theta_2} - 3 \cos(\theta_1-\theta_2) p_{\theta_1}}{16 - 9 \cos^2(\theta_1-\theta_2)}.$

${\dot p_{\theta_1}} = -\frac{1}{2} m \ell^2 \left [ {\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + 3 \frac{g}{\ell} \sin \theta_1 \right ]$

${\dot p_{\theta_2}} = -\frac{1}{2} m \ell^2 \left [ -{\dot \theta_1} {\dot \theta_2} \sin (\theta_1-\theta_2) + \frac{g}{\ell} \sin \theta_2 \right]$

where the generalized momenta are defined by

$$p_{\theta_1} \equiv \frac{\partial L}{\partial \dot{\theta}_1}=ml^2(\frac{4}{3}\dot{\theta}_1+\frac{1}{2}\cos(\theta_1-\theta_2)\dot{\theta}_2)$$$$p_{\theta_2} \equiv \frac{\partial L}{\partial \dot{\theta}_2}=ml^2(\frac{1}{3}\dot{\theta}_2+\frac{1}{2}\cos(\theta_1-\theta_2)\dot{\theta}_1)$$

and $L$ is \begin{eqnarray} L[\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2]=\frac{1}{2}ml^2\left(\frac{4}{3}\dot{\theta}_1^2+\frac{1}{3}\dot{\theta}_2^2+\cos(\theta_1-\theta_2)\dot{\theta}_1\dot{\theta}_2\right)+mgl\left(\frac{3}{2}\cos\theta_1+\frac{1}{2}\cos\theta_2\right) \end{eqnarray}