At latitude $\ell$, the Earth's angular velocity in east–north–up coordinates is $\mathbf{\Omega }=(0,\Omega \cos \ell ,\Omega \sin \ell )$. Assume a small horizontal patch, constant $\Omega$, negligible vertical velocity, and a vertical constraint balancing effective gravity. Absorb the steady centrifugal force into that effective gravity. With $f=2\Omega \sin \ell$ and physical horizontal force per unit mass $F_x^{\prime },F_y^{\prime }$, the Coriolis force givesLarge-scale atmospheric and oceanic flows deflect to their right in the northern hemisphere and left in the southern hemisphere; their balance with a horizontal pressure gradient produces geostrophic flow. A Foucault pendulum has an oscillation plane that precesses relative to the ground at angular rate $-\Omega \sin \ell$. Both effects vanish at the equator within this horizontal approximation.

Suppress primes on the angular velocity while retaining rotating components. The canonical momentum is $\partial L^{\prime }/\partial {\dot{\mathbf{r}}}^{\prime }=m({\dot{\mathbf{r}}}^{\prime }+\mathbf{\Omega }\times {\mathbf{r}}^{\prime })$. Varying $L^{\prime }$ with respect to position and applying the Euler-Lagrange equations givesThe terms are the physical force, the velocity-dependent Coriolis force, the outward centrifugal force, and the Euler force caused by changing angular velocity. The common origin is inertial; an accelerating origin would introduce another translational apparent force.

The inertial Lagrangian is $L=m|\dot{\mathbf{r}}{|}^2/2-U(\mathbf{r},t)$. Put $\mathbf{r}=R(t){\mathbf{r}}^{\prime }$ and define $U^{\prime }({\mathbf{r}}^{\prime },t)=U(R(t){\mathbf{r}}^{\prime },t)$. The transport identity and invariance of the Euclidean norm under a rotation giveIts expansion contains the ordinary kinetic energy, $m{\dot{\mathbf{r}}}^{\prime }\cdot ({\mathbf{\Omega }}^{\prime }\times {\mathbf{r}}^{\prime })$, and $m|{\mathbf{\Omega }}^{\prime }\times {\mathbf{r}}^{\prime }{|}^2/2$. These velocity and position terms generate the apparent forces in the rotating frame.

Write the geometric vector as $\mathbf{V}={\sum }_i{V}_i^{\prime }{\mathbf{e}}_i^{\prime }(t)$. A rotating orthonormal basis obeys ${\dot{\mathbf{e}}}_i^{\prime }=\mathbf{\Omega }\times {\mathbf{e}}_i^{\prime }$. The product rule gives the transport identityBoth terms must be expressed in the same basis before adding them. If $R(t)$ maps rotating components into inertial components, the precise component equation is $d(R{\mathbf{V}}^{\prime })/dt=R({\dot{\mathbf{V}}}^{\prime }+{\mathbf{\Omega }}^{\prime }\times {\mathbf{V}}^{\prime })$. The printed expression suppresses this basis conversion.

The orbiting planet produces a potential $V(r,\vartheta ,\phi -\Omega t)$. It has explicit time dependence and exerts an azimuthal torque, so generally $\dot{E}=V_t$ and ${\dot{L}}_z=-V_{\phi }$ are nonzero. Nevertheless a simultaneous time-translation symmetry $\delta t=\lambda$ and rotation $\delta \phi =\Omega \lambda$ leaves the Lagrangian and $dt$ unchanged. Noether's theorem gives $Q=\Omega p_{\phi }-E$ and henceEquivalently $V_t=-\Omega V_{\phi }$ implies $d(E-\Omega L_z)/dt=0$. This Jacobi integral assumes a rigidly rotating potential with constant $\Omega$ and negligible back-reaction from the test particle.

For a central potential use $L=m|\dot{\mathbf{r}}{|}^2/2-V(|\mathbf{r}|)$. An infinitesimal rotation about any constant unit vector $\mathbf{n}$ has $\varphi =0$ and $\mathbf{\psi }=\mathbf{n}\times \mathbf{r}$. Noether's theorem gives $Q=m\dot{\mathbf{r}}\cdot (\mathbf{n}\times \mathbf{r})=\mathbf{n}\cdot \mathbf{L}$, where $\mathbf{L}=\mathbf{r}\times m\dot{\mathbf{r}}$. Since $\mathbf{n}$ is arbitrary, every component of angular momentum is conserved. A time-translation symmetry has $\varphi =1$, $\mathbf{\psi }=0$, giving $Q=-E$ andas conserved quantities.

The appropriate invariance is invariance of the action element $Ldt$, including the change of the time measure. Under $\delta t=\lambda \varphi$ and $\delta q_i=\lambda {\psi }_i$, differentiation of $q_i^{\prime }(t^{\prime })$ gives $\delta {\dot{q}}_i=\lambda (D_t{\psi }_i-{\dot{q}}_i{D}_t\varphi )$. ThereforeUsing the Euler-Lagrange equations and $D_{t}E=-L_t$, this becomesAn action symmetry therefore yields the Noether chargeFor $\varphi =0$ this reduces to the usual spatial form of Noether's theorem. The same on-shell calculation permits velocity-dependent generators when they define a legitimate variational symmetry. Merely preserving the numerical value of $L$ is insufficient for transformations of time: the Jacobian matrix and determinant $d{t}^{\prime }/dt$ matters. If the action changes by a boundary term $\lambda D_{t}B$, the conserved quantity is $Q-B$.

Define the canonical momentum $p_i=\partial L/\partial {\dot{q}}_i$ and the energy function $E={\sum }_i{p}_i{\dot{q}}_i-L$. Along an Euler-Lagrange equation solution,Thus an autonomous Lagrangian has conserved $E$. For a simple harmonic oscillator, $L=m{\dot{x}}^2/2-k{x}^2/2$ and $\boxed{E=m{\dot{x}}^2/2+k{x}^2/2}$ is constant. For ordinary time-independent mechanical coordinates this is the total mechanical energy.

An explicit numerical method normally has a simple, inexpensive step and needs no simultaneous solve, but stability can force very small steps, especially for stiff differential equations. An implicit numerical method requires a linear or nonlinear solve, storage and suitable solver tolerances; its larger stability region can permit larger steps when stability is the main restriction. Accuracy, phase fidelity and preservation of physical structure still constrain the step. For this oscillator, backward Euler method is stable but dissipative, while forward Euler method is cheap but unstable without damping. Neither statement makes all implicit methods superior to all explicit ones; method order, problem structure and total cost determine the useful choice.

Backward Euler method has update ${\mathbf{y}}^{n+1}=(I-hM{)}^{-1}{\mathbf{y}}^n$ and amplification eigenvalues $g_{\pm }=(1-h{\lambda }_{\pm }{)}^{-1}$. Hencefor every $h\ge 0$. It is therefore unconditionally stable for this linear system. Even with $\gamma =0$, it damps every nonzero numerical mode: stability is obtained at the cost of artificial damping and phase error. Large $h$ does not accurately resolve the oscillation period.

The characteristic polynomial is ${\lambda }^2+2\gamma \lambda +k^2=0$. Since $0<\gamma <|k|$,Forward Euler method has amplification eigenvalues $g_{\pm }=1+h{\lambda }_{\pm }$. Their squared modulus isThusStrict inequality gives numerical damping; at the upper boundary the amplitude is constant despite physical decay, so that boundary is a poor accuracy choice. In the undamped limit $\gamma =0$, $|g{|}^2=1+k^2{h}^2>1$ for every $h>0$: forward Euler method is unstable for any nonzero step on an undamped harmonic oscillator.

A numerical evolution has numerical stability if small perturbations in initial data, forcing or arithmetic are not amplified without control over the time interval of interest. For this constant-coefficient decaying linear system, an absolute-stability test requires the homogeneous numerical amplification factors to have magnitude at most one, with no growing Jordan-block contribution on the unit circle. Stability is distinct from accuracy: a strongly overdamped numerical solution can be stable while representing the physical oscillation poorly.

Put $y^{(0)}=y$ and $y^{(1)}=\dot{y}$. ThenThe second row is the original damped harmonic oscillator equation and the first row defines the auxiliary velocity. This converts a second-order ordinary differential equation into two first-order equations without changing the required two initial values.

For an error Fourier mode $e_j^n={\hat{e}}^n{e}^{i\theta j}$, the periodic discrete Laplacian has eigenvalue $-4{\sin }^2(\theta /2)/h^2$. The Backward Euler method therefore has amplification factorvon Neumann stability analysis thus imposes no upper limit on positive $\Delta \tau$. The constant mode has $g=1$ and is not damped; compatible forcing and a gauge condition remain necessary as discussed in part (e). Large steps can accelerate relaxation of nonconstant modes, but require a suitably accurate linear solver and can worsen its conditioning. Unconditional linear stability is not a statement about time accuracy or cost.

Let $r=\Delta \tau /h^2$ and impose periodic indices modulo $J$. With ${\mathbf{\Phi }}^n=({\Phi }_0^n,\dots ,{\Phi }_{J-1}^n{)}^T$, the linear system iswhere $S$ shifts one grid cell cyclically. For $J\ge 3$, $A$ has diagonal $1+2r$, adjacent and corner entries $-r$, and all other entries zero. It is a cyclic tridiagonal matrix, rather than a nonperiodic tridiagonal one.

The periodic Laplacian annihilates the constant vector, so its steady equation requires $\boxed{\sum _j{\rho }_j=0}$. A strictly positive gravitational density violates this condition: summing the update gives ${\bar{\Phi }}^{n+1}={\bar{\Phi }}^n-4\pi G\Delta \tau \bar{\rho }$, an indefinitely drifting mean. Periodic gravitational calculations normally solve for $\rho -\bar{\rho }$ and fix the arbitrary mean gravitational potential. The per-step matrix $A$ itself remains invertible; the issue is existence of a periodic steady solution of Poisson's equation.

A backward Euler method with spatial second-order central differences givesThe source is fixed in pseudo-time; if it were evolving, its evaluation would also need specifying. Each pseudo-time step therefore requires a linear system solve rather than a purely local explicit update.

At a steady state in pseudo-time, ${\partial }_{\tau }\Phi =0$ and the evolution becomes ${\Phi }_{xx}=4\pi G\rho$, the desired Poisson equation. Deviations from a steady solution obey a diffusion equation, so nonconstant error modes decay. This turns an elliptic boundary-value solve into repeated relaxation steps using familiar time-evolution methods. Boundary data and source compatibility still matter; pseudo-time does not remove a nonexistent steady solution.

For a second-order scalar partial differential equation, classification refers to its principal second-derivative form. An elliptic partial differential equation has a definite principal form, exemplified by Poisson's equation; it typically determines a spatial field from boundary data. A parabolic partial differential equation has a distinguished first-order evolution direction and diffusive second spatial derivatives, exemplified by $u_t=D{u}_{xx}$. Its initial data relax or spread with time. Thus the pseudo-time equation is parabolic, while its steady-state spatial equation is elliptic.

Let the step be $h$. Forward Euler method gives $x_1=x_0+h{v}_0$, $v_1=v_0+hF(x_0)$. Applying the same map with $-h$ gives $x_2=x_1-h{v}_1=x_0-h^{2}F(x_0)$, which generally differs from $x_0$. It is therefore not its own inverse under time reversal.

A kick–drift–kick leapfrog integrator instead hasThe step $-h$ first recovers $v_{1/2}$, then $x_0$, then $v_0$, so $\boxed{{\Phi }_{-h}={\Phi }_h^{-1}}$. This assumes a position-dependent force and the same fixed step. Newtonian mechanics is reversible; respecting that structure helps avoid artificial directional bias and often improves long-term qualitative behaviour. The usual leapfrog also has a structure preserving symplectic geometry, which supports bounded oscillatory energy error, but reversibility by itself does not imply exact conservation of energy. Adaptive stepping requires special care to preserve these properties.

Use a smooth exact solution such as $u(x,t)=e^{-D{q}^{2}t}\sin (qx)$ on a domain with matching boundary conditions, and compare numerical and exact values at a fixed physical time. To measure spatial order, halve $h$ while making temporal error negligible; the error norm should fall by a factor of four. To measure temporal order, halve $k$ on a sufficiently fine fixed spatial grid and again expect a factor of four. Report observed order $p={\log }_2(E_h/E_{h/2})$ and use several refinements to identify the asymptotic regime.

Check more than one smooth mode, boundary implementation and expected conservation or decay properties. A manufactured solution with an added known source can test nontrivial boundaries and forcing. Repeat in maximum and root-mean-square norms, stop before finite-precision or solver tolerances dominate, and separate stability failures from errors in convergence order. A single small error or a single visually smooth solution does not establish accuracy.