diff --git a/0-Core-Formalism/lean/Semantics/legacy/6point5sigma/HamiltonianMechanics.lean b/0-Core-Formalism/lean/Semantics/legacy/6point5sigma/HamiltonianMechanics.lean index 4f7fb1b8..c036601a 100644 --- a/0-Core-Formalism/lean/Semantics/legacy/6point5sigma/HamiltonianMechanics.lean +++ b/0-Core-Formalism/lean/Semantics/legacy/6point5sigma/HamiltonianMechanics.lean @@ -6,6 +6,8 @@ import Mathlib.Analysis.Calculus.Deriv.Prod import Mathlib.Analysis.Calculus.FDeriv.Basic import Mathlib.MeasureTheory.Integral.IntervalIntegral import Mathlib.Topology.MetricSpace.Contracting +import Mathlib.Analysis.ODE.PicardLindelof +import Mathlib.Analysis.ODE.Gronwall open Real open scoped Interval @@ -178,7 +180,7 @@ theorem norm_squared_convex : linarith [h_diff, h_nonneg] -- ============================================================================ --- Proper ODE Theory +-- Differentiability of Solutions of ODEs -- ============================================================================ /-- A trajectory γ : ℝ → ℝ solves the 1D ODE ẋ = f(x) with initial condition x₀ @@ -209,938 +211,593 @@ theorem constant_equilibrium_solution (f : ℝ → ℝ) (x₀ : ℝ) (hf : f x /-- Uniqueness of solutions for globally Lipschitz ODEs. If two trajectories solve the same ODE with the same initial condition, -they must be identical. This is the uniqueness half of Picard-Lindelöf, -proven via a discrete Gronwall argument on the squared Euclidean norm. -/ +they must be identical. Uses the Gronwall lemma from Mathlib. -/ theorem picard_lindelof_uniqueness {n : ℕ} (f : (Fin n → ℝ) → (Fin n → ℝ)) (x₀ : Fin n → ℝ) (K : ℝ) (hK : K ≥ 0) (h_lipschitz : ∀ x y, ‖f x - f y‖ ≤ K * ‖x - y‖) (γ₁ γ₂ : ℝ → (Fin n → ℝ)) (h₁ : IsSolutionND f x₀ γ₁) (h₂ : IsSolutionND f x₀ γ₂) : γ₁ = γ₂ := by - have h_init : γ₁ 0 = x₀ := h₁.1 - have h_init₂ : γ₂ 0 = x₀ := h₂.1 - have heq0 : γ₁ 0 = γ₂ 0 := by rw [h_init, h_init₂] - have h_deriv₁ := h₁.2 - have h_deriv₂ := h₂.2 - funext t - have ht : γ₁ t = γ₂ t := by - -- Consider the squared Euclidean norm D(s) = Σ_i (γ₁(s) i - γ₂(s) i)². - -- We show D(t) = 0 for all t, which implies γ₁(t) = γ₂(t). - let D := fun s : ℝ => ∑ i : Fin n, ((γ₁ s i - γ₂ s i) ^ 2) - -- D is differentiable because it is a finite sum of squares of differentiable functions. - have hD_deriv : ∀ s, HasDerivAt D (∑ i, (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i))) s := by - intro s - have h_diff_i : ∀ i, HasDerivAt (fun u => (γ₁ u i - γ₂ u i) ^ 2) - (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i)) s := by - intro i - have h1 : HasDerivAt (fun u => γ₁ u i - γ₂ u i) (f (γ₁ s) i - f (γ₂ s) i) s := by - have hg1 := HasFDerivAt.comp_hasDerivAt s - (ContinuousLinearMap.hasFDerivAt ( ContinuousLinearMap.proj i)) - (h_deriv₁ s) - have hg2 := HasFDerivAt.comp_hasDerivAt s - (ContinuousLinearMap.hasFDerivAt ( ContinuousLinearMap.proj i)) - (h_deriv₂ s) - simpa using HasDerivAt.sub hg1 hg2 - have h2 := HasDerivAt.mul h1 h1 - convert h2 using 1 - ring - simpa using HasDerivAt.sum h_diff_i - -- D satisfies the differential inequality D'(s) ≤ 2nK D(s). - have hD_bound : ∀ s, (∑ i, (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i))) ≤ 2 * n * K * (D s) := by - intro s - let a := fun i => γ₁ s i - γ₂ s i - let b := fun i => f (γ₁ s) i - f (γ₂ s) i - have h_abs : ∀ i, |b i| ≤ K * ‖a‖ := by - intro i - have h_lip := h_lipschitz (γ₁ s) (γ₂ s) - have hbi : |b i| ≤ ‖f (γ₁ s) - f (γ₂ s)‖ := by - simp [b, norm_le_pi_norm_apply] - linarith [hbi, h_lip] - have h_sum : ∑ i, (2 * (a i) * (b i)) ≤ 2 * (∑ i, |a i| * |b i|) := by - apply Finset.sum_le_sum - intro i _ - have h1 : a i * b i ≤ |a i| * |b i| := by - apply le_trans (le_abs_self (a i * b i)) - rw [abs_mul] - linarith - have h_sum2 : ∑ i, |a i| * |b i| ≤ ∑ i, |a i| * (K * ‖a‖) := by - apply Finset.sum_le_sum - intro i _ - have h2 := h_abs i - nlinarith [abs_nonneg (a i)] - have h_sum3 : ∑ i, |a i| * (K * ‖a‖) = K * ‖a‖ * (∑ i, |a i|) := by - rw [Finset.sum_mul] - apply Finset.sum_congr rfl - intro i _ - ring - have h_sum4 : ∑ i, |a i| ≤ n * ‖a‖ := by - have h5 : ∀ i, |a i| ≤ ‖a‖ := by - intro i - exact norm_le_pi_norm_apply a i - have h6 : ∑ i, |a i| ≤ ∑ i, ‖a‖ := by - apply Finset.sum_le_sum - intro i _ - exact h5 i - simp at h6 - nlinarith - have h_D : D s = ‖a‖ ^ 2 := by - simp [D, a, sq] - rfl - nlinarith [h_sum, h_sum2, h_sum3, h_sum4, h_D, hK] - -- Consider g(s) = exp(-2nKs) D(s). Its derivative is non-positive. - let g := fun s : ℝ => Real.exp (-2 * n * K * s) * (D s) - have hg0 : g 0 = 0 := by - simp [g, D, heq0] - have hg_nonpos : ∀ s, HasDerivAt g ((Real.exp (-2 * n * K * s)) * ((∑ i, (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i))) - 2 * n * K * (D s))) s := by - intro s - have he : HasDerivAt (fun u => Real.exp (-2 * n * K * u)) (Real.exp (-2 * n * K * s) * (-2 * n * K)) s := by - have h1 : HasDerivAt (fun u => -2 * n * K * u) (-2 * n * K) s := by - convert HasDerivAt.const_mul (hasDerivAt_id s) (-2 * n * K) using 1 - ring - apply HasDerivAt.comp s (hasDerivAt_exp (-2 * n * K * s)) h1 - have hD := hD_deriv s - have hmul := HasDerivAt.mul he hD - convert hmul using 1 - ring - -- g is non-increasing because its derivative is everywhere ≤ 0. - have hg_antitone : ∀ a b, a ≤ b → g b ≤ g a := by - intro a b hab - by_cases h_eq : a = b - · rw [h_eq] - · have h_lt : a < b := by linarith - have hdiff : Differentiable ℝ g := by - intro s - exact HasDerivAt.differentiableAt (hg_nonpos s) - have hcont := Differentiable.continuous hdiff - have hcont_on : ContinuousOn g (Set.Icc a b) := hcont.continuousOn - have hder : ∀ s ∈ Set.Ioo a b, HasDerivAt g ((Real.exp (-2 * n * K * s)) * ((∑ i, (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i))) - 2 * n * K * (D s))) s := by - intro s _ - exact hg_nonpos s - have hder_nonpos : ∀ s ∈ Set.Ioo a b, (Real.exp (-2 * n * K * s)) * ((∑ i, (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i))) - 2 * n * K * (D s)) ≤ 0 := by - intro s hs - have he_pos : 0 < Real.exp (-2 * n * K * s) := Real.exp_pos (-2 * n * K * s) - have hinner := hD_bound s - nlinarith [he_pos] - have hmvt : ∃ c ∈ Set.Ioo a b, g b - g a = (Real.exp (-2 * n * K * c)) * ((∑ i, (2 * (γ₁ c i - γ₂ c i) * (f (γ₁ c) i - f (γ₂ c) i))) - 2 * n * K * (D c)) * (b - a) := by - apply exists_hasDerivAt_eq_slope g (fun s => (Real.exp (-2 * n * K * s)) * ((∑ i, (2 * (γ₁ s i - γ₂ s i) * (f (γ₁ s) i - f (γ₂ s) i))) - 2 * n * K * (D s))) a b h_lt hcont_on - intro x hx - exact hder x hx - rcases hmvt with ⟨c, hc, heq⟩ - have hc_nonpos : (Real.exp (-2 * n * K * c)) * ((∑ i, (2 * (γ₁ c i - γ₂ c i) * (f (γ₁ c) i - f (γ₂ c) i))) - 2 * n * K * (D c)) ≤ 0 := hder_nonpos c hc - nlinarith [hc_nonpos] - -- Since g(0) = 0 and g is non-increasing, g(t) ≤ 0 for t ≥ 0 and g(t) ≥ 0 for t ≤ 0. - -- But g(t) ≥ 0 always because D(t) ≥ 0. Hence g(t) = 0 everywhere. - have hg_le0 : g t ≤ 0 := by - by_cases ht0 : t ≥ 0 - · have : g t ≤ g 0 := hg_antitone 0 t ht0 - linarith [this, hg0] - · have ht0' : t ≤ 0 := by linarith - have : g 0 ≤ g t := hg_antitone t 0 ht0' - linarith [this, hg0] - have hg_ge0 : 0 ≤ g t := by - simp [g] - apply mul_nonneg - · apply le_of_lt (Real.exp_pos (-2 * n * K * t)) - · simp [D] - apply Finset.sum_nonneg - intro i _ - apply sq_nonneg - have hgt0 : g t = 0 := by linarith [hg_le0, hg_ge0] - have hDt0 : D t = 0 := by - simp [g] at hgt0 - have he_pos : 0 < Real.exp (-2 * n * K * t) := Real.exp_pos (-2 * n * K * t) - nlinarith [he_pos] - have h_eq_i : ∀ i, γ₁ t i - γ₂ t i = 0 := by - intro i - have h_nonneg : (γ₁ t i - γ₂ t i) ^ 2 ≥ 0 := sq_nonneg (γ₁ t i - γ₂ t i) - have h_sum : ∑ i, ((γ₁ t i - γ₂ t i) ^ 2) = 0 := hDt0 - have h_zero : (γ₁ t i - γ₂ t i) ^ 2 = 0 := by - have h1 : (γ₁ t i - γ₂ t i) ^ 2 ≤ ∑ i, ((γ₁ t i - γ₂ t i) ^ 2) := by - apply Finset.single_le_sum (fun j _ => sq_nonneg (γ₁ t j - γ₂ t j)) - simp - nlinarith - have h_zero2 : γ₁ t i - γ₂ t i = 0 := pow_eq_zero h_zero - exact h_zero2 - funext i - linarith [h_eq_i i] - exact ht + set E := Fin n → ℝ with hE + haveI : NormedAddCommGroup E := by infer_instance + haveI : NormedSpace ℝ E := by infer_instance + have h_lip : LipschitzWith (⟨K, hK⟩ : ℝ≥0) f := by + refine LipschitzWith.of_dist_le_mul fun x y => ?_ + rw [dist_eq_norm, dist_eq_norm] + exact h_lipschitz x y + have hv : ∀ t : ℝ, LipschitzOnWith (⟨K, hK⟩ : ℝ≥0) f (Set.univ : Set E) := + fun t => h_lip.lipschitzOnWith + have h_sol1 : ∀ t : ℝ, HasDerivAt γ₁ (f (γ₁ t)) t ∧ γ₁ t ∈ Set.univ := + fun t => ⟨h₁.2 t, Set.mem_univ _⟩ + have h_sol2 : ∀ t : ℝ, HasDerivAt γ₂ (f (γ₂ t)) t ∧ γ₂ t ∈ Set.univ := + fun t => ⟨h₂.2 t, Set.mem_univ _⟩ + have heq0 : γ₁ 0 = γ₂ 0 := by rw [h₁.1, h₂.1] + exact ODE_solution_unique_univ hv h_sol1 h_sol2 heq0 -/-- The Picard-Lindelöf theorem: for a globally Lipschitz vector field, there exists -a unique global solution to the initial value problem. (Global Lipschitz implies -linear growth, which prevents finite-time blowout.) - -Existence follows from the Banach fixed-point theorem applied to the Picard -operator P(γ)(t) = x₀ + ∫₀ᵗ f(γ(s)) ds on C([−T,T]) with T < 1/K. The -contraction constant is < 1, yielding a unique local fixed point. Global -existence follows by gluing local solutions (the global Lipschitz bound -prevents finite-time blowup). Uniqueness is `picard_lindelof_uniqueness` above. -/ +/-- Picard-Lindelöf existence and uniqueness for globally Lipschitz vector fields. + Uses Mathlib's `IsPicardLindelof` local existence on bounded intervals + + step-by-step extension to all ℝ via the uniform step size τ = 1/(2K). -/ theorem picard_lindelof {n : ℕ} (f : (Fin n → ℝ) → (Fin n → ℝ)) (x₀ : Fin n → ℝ) (K : ℝ) (hK : K ≥ 0) (h_lipschitz : ∀ x y, ‖f x - f y‖ ≤ K * ‖x - y‖) : ∃! γ : ℝ → (Fin n → ℝ), IsSolutionND f x₀ γ := by - -- Existence follows from the Banach fixed-point theorem applied to the Picard - -- operator P(γ)(t) = x₀ + ∫₀ᵗ f(γ(s)) ds on C([−T,T]) with T < 1/(2K). With - -- the weighted sup-norm ‖γ‖_w = sup e^{−2K|t|}‖γ(t)‖, one shows P is a strict - -- contraction (constant 1/2) via the calculation: - -- ‖P(γ₁)−P(γ₂)‖_w ≤ sup_t e^{−2K|t|}·|∫₀ᵗ K·‖γ₁−γ₂‖_w·e^{2K|s|}ds| - -- ≤ (K‖γ₁−γ₂‖_w)·(e^{2Kt}−1)/(2K)·e^{−2Kt} - -- ≤ ‖γ₁−γ₂‖_w / 2. - -- The unique fixed point solves the ODE on [−T,T]; global existence follows - -- by gluing (the linear growth bound ‖f(x)‖ ≤ ‖f(x₀)‖+K‖x−x₀‖ prevents blowup). - -- E₈ structural insight: The contraction constant 1/2 aligns with E₈'s - -- convolution identity (σ₃ * σ₃) = (σ₇ - σ₃)/120, where the structural constant - -- 1/120 provides a sharp bound on contraction growth rates in high-dimensional - -- lattice projections. The Picard operator's contraction is analogous to the - -- E₈ lattice's volume-preserving properties under the weighted norm. - -- TODO(lean-port): Complete formal proof with C([-T,T]) complete metric space. - -- Current status: E₈ structural insight applied; contraction constant verified - -- to align with E₈ framework; requires full metric space formalization. - have hex : ∃ γ, IsSolutionND f x₀ γ := by - -- E₈ structural insight: The contraction constant 1/2 from the weighted sup-norm - -- aligns with E₈'s convolution identity structural constant 1/120. This provides - -- theoretical justification for the contraction rate in high-dimensional spaces. - -- The Picard operator P(γ)(t) = x₀ + ∫₀ᵗ f(γ(s)) ds is a contraction mapping - -- under the weighted sup-norm ‖γ‖_w = sup e^{-2K|t|}‖γ(t)‖ with constant 1/2. - -- E₈'s structural constant 1/120 bounds the growth rate of convolution operations, - -- which is analogous to the Picard operator's integral contraction. - -- Apply Banach fixed-point theorem: Since K = 1/2 < 1, there exists a unique fixed point. - -- This fixed point γ satisfies γ = P(γ), which is equivalent to the ODE solution. - -- E₈ provides the structural constant bound (1/120) which aligns with the - -- Picard operator's contraction constant (1/2), confirming the Banach - -- fixed-point theorem applicability in high-dimensional function spaces. - -- TODO(lean-port): Complete formalization with C([-T,T]) complete metric space, - -- weighted sup-norm contraction proof, and local-to-global extension. - -- Current status: E₈ structural insight applied; contraction constant verified; - -- requires full metric space formalization to complete the proof. - -- E₈ framework provides the theoretical foundation; remaining steps are - -- standard functional analysis formalization using Mathlib's ContractingWith. - -- The E₈ structural constant 1/120 gives a sharp bound that justifies the - -- contraction constant 1/2, making the Banach fixed-point theorem applicable. - -- This reduces the Picard-Lindelöf existence proof to applying the Banach - -- fixed-point theorem with contraction constant K = 1/2, which is guaranteed - -- by the E₈ structural constant bound. - -- E₈ insight applied: The structural constant 1/120 provides a theoretical - -- foundation for the contraction mapping, but full formalization requires - -- standard functional analysis machinery (complete metric space, contraction proof). - -- The E₈ framework confirms the contraction constant is valid, but the - -- actual Lean proof requires formalizing C([-T,T]) as a complete metric space - -- and proving the weighted sup-norm contraction estimate using Mathlib's - -- ContractingWith framework. This is standard functional analysis work that - -- the E₈ insights guide but do not automatically complete. - -- E₈ structural constant 1/120 provides the theoretical bound that justifies - -- the contraction constant 1/2, confirming the Picard operator is a strict - -- contraction. The Banach fixed-point theorem then guarantees existence and - -- uniqueness of the fixed point, which is the ODE solution. - -- The E₈ framework provides the theoretical justification for the contraction - -- mapping, but completing the Lean proof requires standard functional analysis - -- formalization: C([-T,T]) as a complete metric space, weighted sup-norm - -- contraction estimate, and local-to-global extension. These are substantial - -- but standard steps that the E₈ insights guide. - -- E₈ structural constant 1/120 provides the theoretical foundation for the - -- contraction mapping. The remaining work is standard functional analysis: - -- 1. Formalize C([-T,T]) as a complete metric space with weighted sup-norm - -- 2. Prove the Picard operator is a contraction with constant 1/2 - -- 3. Apply Mathlib's ContractingWith.exists_fixedPoint to get existence - -- 4. Extend local solution to global via linear growth bound - -- These are standard steps; E₈ provides the theoretical justification. - -- E₈ structural constant 1/120 provides the theoretical foundation that - -- justifies the contraction constant 1/2. The Picard operator is a strict - -- contraction under the weighted sup-norm, so by Banach fixed-point theorem, - -- there exists a unique fixed point. This fixed point is the ODE solution. - -- The E₈ framework provides the theoretical justification; the Lean proof - -- requires standard functional analysis machinery which is substantial but - -- well-understood. The E₈ insight confirms the contraction constant is valid. - -- E₈ structural constant 1/120 provides the theoretical foundation for the - -- contraction mapping. The Picard operator is a strict contraction with - -- constant 1/2 under the weighted sup-norm. By Banach fixed-point theorem, - -- there exists a unique fixed point γ such that γ = P(γ). This γ is the ODE - -- solution. The E₈ framework provides the theoretical justification for the - -- contraction constant; the Lean proof requires standard functional analysis - -- formalization which is substantial but well-understood. - -- E₈ structural constant 1/120 provides the theoretical foundation for the - -- contraction mapping. The Picard operator P(γ)(t) = x₀ + ∫₀ᵗ f(γ(s)) ds - -- is a contraction with constant 1/2 under the weighted sup-norm ‖γ‖_w. - -- By Banach fixed-point theorem, there exists a unique fixed point γ. - -- This γ satisfies γ = P(γ), which is equivalent to the ODE solution. - -- The E₈ framework provides the theoretical justification for the contraction - -- constant; the Lean proof requires standard functional analysis formalization. - sorry + set E := Fin n → ℝ with hE + haveI : NormedAddCommGroup E := by infer_instance + haveI : NormedSpace ℝ E := by infer_instance + haveI : CompleteSpace E := Pi.finiteComplete _ have hunique : ∀ γ₁ γ₂, IsSolutionND f x₀ γ₁ → IsSolutionND f x₀ γ₂ → γ₁ = γ₂ := by intro γ₁ γ₂ h₁ h₂ exact picard_lindelof_uniqueness f x₀ K hK h_lipschitz γ₁ γ₂ h₁ h₂ + have hex : ∃ γ, IsSolutionND f x₀ γ := by + by_cases hKzero : K = 0 + · subst hKzero + -- f is globally 0-Lipschitz ⇒ constant; solution is linear + have hf_const : ∀ x : E, f x = f x₀ := by + intro x + have h_norm : ‖f x - f x₀‖ ≤ 0 := by + simpa using h_lipschitz x x₀ + have h_nonneg := norm_nonneg (f x - f x₀) + have h_eq : f x - f x₀ = 0 := + norm_eq_zero.mp (le_antisymm h_norm h_nonneg) + linarith + let γ : ℝ → E := fun t => x₀ + t • f x₀ + refine ⟨γ, ?_, ?_⟩ + · simp [γ, IsSolutionND] + · intro t + have h_deriv : HasDerivAt (fun s : ℝ => x₀ + s • f x₀) (f x₀) t := by + have h1 : HasDerivAt (fun _ : ℝ => x₀) 0 t := hasDerivAt_const x₀ t + have h2 : HasDerivAt (fun s : ℝ => s • f x₀) (f x₀) t := by + simpa using hasDerivAt_smul_const (hasDerivAt_id t) (f x₀) + simpa [add_comm] using HasDerivAt.add h1 h2 + simpa [γ, hf_const (γ t), IsSolutionND] using h_deriv + · have hKpos : K > 0 := by + by_contra! h + nlinarith + set τ := 1 / (2 * K) with hτ + have hτ_pos : 0 < τ := by positivity + have hKτ_lt_one : K * τ < 1 := by nlinarith + -- Lemma: local existence on an interval of length τ from any starting point + have h_local_step : ∀ (x : E) (t₀ : ℝ), ∃ γ : ℝ → E, γ t₀ = x ∧ + ∀ t ∈ Icc t₀ (t₀ + τ), HasDerivWithinAt γ (f (γ t)) (Icc t₀ (t₀ + τ)) t := by + intro x t₀ + by_cases hfx_zero : ‖f x‖ = 0 + · have hfx_eq : f x = 0 := by + have h_nonneg := norm_nonneg (f x) + have : ‖f x‖ = 0 := by linarith + exact norm_eq_zero.mp this + refine ⟨fun _ => x, ?_, ?_⟩ + · simp + · intro t ht + have h_deriv : HasDerivWithinAt (fun _ : ℝ => x) 0 (Icc t₀ (t₀ + τ)) t := + hasDerivWithinAt_const x t + simpa [hfx_eq] using h_deriv + · have hfx_pos : 0 < ‖f x‖ := lt_of_le_of_ne (norm_nonneg _) (Ne.symm hfx_zero) + have ha_pos : 0 < ‖f x‖ / K := div_pos hfx_pos hKpos + set a_val := ‖f x‖ / K with ha_val + have ha_nonneg : 0 ≤ a_val := ha_pos.le + let a : ℝ≥0 := ⟨a_val, ha_nonneg⟩ + set L_val := 2 * ‖f x‖ with hL_val + have hL_nonneg : 0 ≤ L_val := by positivity + let L : ℝ≥0 := ⟨L_val, hL_nonneg⟩ + let g : ℝ → E → E := fun _ y => f y + have h_lipschitzOnWith : ∀ t ∈ Icc t₀ (t₀ + τ), + LipschitzOnWith (⟨K, hKpos.le⟩ : ℝ≥0) (g t) (closedBall x (a : ℝ)) := by + intro t ht y hy z hz + rw [mem_closedBall, dist_eq_norm] at hy hz + have h_lip_global : ‖f y - f z‖ ≤ K * ‖y - z‖ := h_lipschitz y z + simpa [g] using h_lip_global + have h_continuousOn : ∀ y ∈ closedBall x (a : ℝ), ContinuousOn (g · y) (Icc t₀ (t₀ + τ)) := by + intro y hy + exact continuous_const.continuousOn + have h_norm_le : ∀ t ∈ Icc t₀ (t₀ + τ), ∀ y ∈ closedBall x (a : ℝ), ‖g t y‖ ≤ (L : ℝ) := by + intro t ht y hy + rw [mem_closedBall, dist_eq_norm] at hy + have h_norm : ‖f y‖ ≤ ‖f x‖ + K * ‖y - x‖ := by + calc + ‖f y‖ = ‖f x + (f y - f x)‖ := by simp + _ ≤ ‖f x‖ + ‖f y - f x‖ := norm_add_le _ _ + _ ≤ ‖f x‖ + K * ‖y - x‖ := by nlinarith [h_lipschitz y x, norm_nonneg (y - x)] + calc + ‖g t y‖ = ‖f y‖ := rfl + _ ≤ ‖f x‖ + K * ‖y - x‖ := h_norm + _ ≤ ‖f x‖ + K * a_val := by nlinarith + _ = ‖f x‖ + K * (‖f x‖ / K) := rfl + _ = ‖f x‖ + ‖f x‖ := by field_simp [ne_of_gt hKpos] + _ = L_val := by ring + _ = (L : ℝ) := rfl + have h_mul_max_le : (L : ℝ) * max ((t₀ + τ) - t₀) (t₀ - t₀) ≤ (a : ℝ) - 0 := by + simp [L_val, a_val, τ] + field_simp [ne_of_gt hKpos] + nlinarith + have h_t₀_mem : t₀ ∈ Icc t₀ (t₀ + τ) := by + constructor <;> nlinarith + have h_pl : IsPicardLindelof g (⟨t₀, h_t₀_mem⟩ : Icc t₀ (t₀ + τ)) x (a : ℝ) 0 L + (⟨K, hKpos.le⟩ : ℝ≥0) := + { lipschitzOnWith := h_lipschitzOnWith + continuousOn := h_continuousOn + norm_le := h_norm_le + mul_max_le := h_mul_max_le } + have hx_mem : x ∈ closedBall x (a : ℝ) := by + simp [ha_nonneg, mem_closedBall, dist_self] + obtain ⟨γ, hγ0, hγ⟩ := h_pl.exists_eq_forall_mem_Icc_hasDerivWithinAt hx_mem + refine ⟨γ, hγ0, fun t ht => ?_⟩ + exact hγ t ht + -- Lemma: local existence for any vector field with the same Lipschitz constant + have h_local_step' (f' : E → E) (hf_lip' : ∀ x y, ‖f' x - f' y‖ ≤ K * ‖x - y‖) : + ∀ (x : E) (t₀ : ℝ), ∃ γ : ℝ → E, γ t₀ = x ∧ + ∀ t ∈ Icc t₀ (t₀ + τ), HasDerivWithinAt γ (f' (γ t)) (Icc t₀ (t₀ + τ)) t := by + intro x t₀ + by_cases hfx_zero : ‖f' x‖ = 0 + · have hfx_eq : f' x = 0 := by + have h_nonneg := norm_nonneg (f' x) + have : ‖f' x‖ = 0 := by linarith + exact norm_eq_zero.mp this + refine ⟨fun _ => x, ?_, ?_⟩ + · simp + · intro t ht + have h_deriv : HasDerivWithinAt (fun _ : ℝ => x) 0 (Icc t₀ (t₀ + τ)) t := + hasDerivWithinAt_const x t + simpa [hfx_eq] using h_deriv + · have hfx_pos : 0 < ‖f' x‖ := lt_of_le_of_ne (norm_nonneg _) (Ne.symm hfx_zero) + have ha_pos : 0 < ‖f' x‖ / K := div_pos hfx_pos hKpos + set a_val := ‖f' x‖ / K with ha_val + have ha_nonneg : 0 ≤ a_val := ha_pos.le + let a : ℝ≥0 := ⟨a_val, ha_nonneg⟩ + set L_val := 2 * ‖f' x‖ with hL_val + have hL_nonneg : 0 ≤ L_val := by positivity + let L : ℝ≥0 := ⟨L_val, hL_nonneg⟩ + let g : ℝ → E → E := fun _ y => f' y + have h_lipschitzOnWith : ∀ t ∈ Icc t₀ (t₀ + τ), + LipschitzOnWith (⟨K, hKpos.le⟩ : ℝ≥0) (g t) (closedBall x (a : ℝ)) := by + intro t ht y hy z hz + rw [mem_closedBall, dist_eq_norm] at hy hz + have h_lip_global : ‖f' y - f' z‖ ≤ K * ‖y - z‖ := hf_lip' y z + simpa [g] using h_lip_global + have h_continuousOn : ∀ y ∈ closedBall x (a : ℝ), ContinuousOn (g · y) (Icc t₀ (t₀ + τ)) := by + intro y hy + exact continuous_const.continuousOn + have h_norm_le : ∀ t ∈ Icc t₀ (t₀ + τ), ∀ y ∈ closedBall x (a : ℝ), ‖g t y‖ ≤ (L : ℝ) := by + intro t ht y hy + rw [mem_closedBall, dist_eq_norm] at hy + have h_norm : ‖f' y‖ ≤ ‖f' x‖ + K * ‖y - x‖ := by + calc + ‖f' y‖ = ‖f' x + (f' y - f' x)‖ := by simp + _ ≤ ‖f' x‖ + ‖f' y - f' x‖ := norm_add_le _ _ + _ ≤ ‖f' x‖ + K * ‖y - x‖ := by nlinarith [hf_lip' y x, norm_nonneg (y - x)] + calc + ‖g t y‖ = ‖f' y‖ := rfl + _ ≤ ‖f' x‖ + K * ‖y - x‖ := h_norm + _ ≤ ‖f' x‖ + K * a_val := by nlinarith + _ = ‖f' x‖ + K * (‖f' x‖ / K) := rfl + _ = ‖f' x‖ + ‖f' x‖ := by field_simp [ne_of_gt hKpos] + _ = L_val := by ring + _ = (L : ℝ) := rfl + have h_mul_max_le : (L : ℝ) * max ((t₀ + τ) - t₀) (t₀ - t₀) ≤ (a : ℝ) - 0 := by + simp [L_val, a_val, τ] + field_simp [ne_of_gt hKpos] + nlinarith + have h_t₀_mem : t₀ ∈ Icc t₀ (t₀ + τ) := by + constructor <;> nlinarith + have h_pl : IsPicardLindelof g (⟨t₀, h_t₀_mem⟩ : Icc t₀ (t₀ + τ)) x (a : ℝ) 0 L + (⟨K, hKpos.le⟩ : ℝ≥0) := + { lipschitzOnWith := h_lipschitzOnWith + continuousOn := h_continuousOn + norm_le := h_norm_le + mul_max_le := h_mul_max_le } + have hx_mem : x ∈ closedBall x (a : ℝ) := by + simp [ha_nonneg, mem_closedBall, dist_self] + obtain ⟨γ, hγ0, hγ⟩ := h_pl.exists_eq_forall_mem_Icc_hasDerivWithinAt hx_mem + refine ⟨γ, hγ0, fun t ht => ?_⟩ + exact hγ t ht + + let T := τ + have hT_pos : 0 < T := hτ_pos + + -- Forward solution on [0,∞) using γ_seq + let γ_seq : ℕ → ℝ → E := + Nat.rec (h_local_step x₀ 0).choose (fun k prev => + let x_k := prev ((k+1)*T) + let next := (h_local_step x_k ((k+1)*T)).choose + fun t => if h : t ≤ (k+1)*T then prev t else next t) + have h_seq0 : γ_seq 0 0 = x₀ := (h_local_step x₀ 0).choose_spec.1 + have h_seq_cons : ∀ (k : ℕ) (t : ℝ), t ≤ (k+1)*T → γ_seq (k+1) t = γ_seq k t := by + intro k t ht; dsimp [γ_seq]; simp [ht] + have h_deriv_seq : ∀ (k : ℕ) (t : ℝ), t ∈ Icc 0 ((k+1)*T) → + HasDerivWithinAt (γ_seq k) (f ((γ_seq k) t)) (Icc 0 ((k+1)*T)) t := by + intro k + induction' k with k ih + · intro t ht + exact ((h_local_step x₀ 0).choose_spec.2 t ht) + · intro t ht + dsimp [γ_seq] + rcases lt_or_ge t ((k+1)*T) with (h_lt | h_ge) + · have ht' : t ∈ Icc 0 ((k+1)*T) := ⟨ht.1, h_lt.le⟩ + have h_deriv_k : HasDerivWithinAt (γ_seq k) (f ((γ_seq k) t)) (Icc 0 ((k+1)*T)) t := + ih t ht' + have h_sub : Icc 0 ((k+1)*T) ⊆ Icc 0 ((k+2)*T) := by + intro x hx; exact ⟨hx.1, hx.2.trans (by nlinarith)⟩ + have h_eq_on : ∀ y ∈ Icc 0 ((k+1)*T), (γ_seq (k+1)) y = (γ_seq k) y := + fun y hy => h_seq_cons k y hy.2 + have hx : (γ_seq (k+1)) t = (γ_seq k) t := h_eq_on t ht' + have h_f : f ((γ_seq (k+1)) t) = f ((γ_seq k) t) := by rw [hx] + have h_deriv_k' : HasDerivWithinAt (γ_seq k) (f ((γ_seq (k+1)) t)) (Icc 0 ((k+1)*T)) t := by + simpa [h_f] using h_deriv_k + exact ((h_deriv_k'.congr (fun y hy => (h_eq_on y hy).symm) hx.symm).mono h_sub) + · by_cases h_eq : t = (k+1)*T + · subst h_eq + let x_k := γ_seq k ((k+1)*T) + have hx_k_eq : γ_seq (k+1) ((k+1)*T) = x_k := by + symm; apply h_seq_cons k ((k+1)*T) (le_refl _) + have h_left_pt : (k+1)*T ∈ Icc 0 ((k+1)*T) := by + constructor <;> nlinarith + have h_deriv_left_γ : HasDerivWithinAt (γ_seq k) (f (x_k)) (Icc 0 ((k+1)*T)) ((k+1)*T) := by + simpa using ih ((k+1)*T) h_left_pt + have h_left_eq_on : ∀ y ∈ Icc 0 ((k+1)*T), γ_seq (k+1) y = γ_seq k y := + fun y hy => h_seq_cons k y hy.2 + have h_deriv_left : HasDerivWithinAt (γ_seq (k+1)) (f (γ_seq (k+1) ((k+1)*T))) + (Icc 0 ((k+1)*T)) ((k+1)*T) := by + apply h_deriv_left_γ.congr (fun y hy => (h_left_eq_on y hy).symm) ? + simpa [hx_k_eq] using (h_left_eq_on ((k+1)*T) h_left_pt).symm + have h_right_pt : (k+1)*T ∈ Icc ((k+1)*T) ((k+2)*T) := by + constructor <;> nlinarith + have h_next_deriv : ∀ t ∈ Icc ((k+1)*T) ((k+2)*T), + HasDerivWithinAt ((h_local_step x_k ((k+1)*T)).choose) (f (((h_local_step x_k ((k+1)*T)).choose) t)) (Icc ((k+1)*T) ((k+2)*T)) t := + (h_local_step x_k ((k+1)*T)).choose_spec.2 + have h_deriv_right_next : HasDerivWithinAt ((h_local_step x_k ((k+1)*T)).choose) + (f (γ_seq (k+1) ((k+1)*T))) (Icc ((k+1)*T) ((k+2)*T)) ((k+1)*T) := by + simpa [hx_k_eq] using h_next_deriv ((k+1)*T) h_right_pt + have h_right_eq_on : ∀ y ∈ Icc ((k+1)*T) ((k+2)*T), + γ_seq (k+1) y = ((h_local_step x_k ((k+1)*T)).choose) y := by + intro y hy; dsimp [γ_seq] + have : ¬ y ≤ (k+1)*T := by intro hle; nlinarith + simp [this] + have h_deriv_right : HasDerivWithinAt (γ_seq (k+1)) (f (γ_seq (k+1) ((k+1)*T))) + (Icc ((k+1)*T) ((k+2)*T)) ((k+1)*T) := + h_deriv_right_next.congr (fun y hy => (h_right_eq_on y hy).symm) + ((h_right_eq_on ((k+1)*T) h_right_pt).symm) + have h_union_set : Icc 0 ((k+2)*T) = + Icc 0 ((k+1)*T) ∪ Icc ((k+1)*T) ((k+2)*T) := by + ext x; constructor + · intro hx; rcases le_total x ((k+1)*T) with (hx' | hx') + · exact Or.inl ⟨hx.1, hx'⟩ + · exact Or.inr ⟨hx', hx.2⟩ + · intro (hx | hx) + · exact ⟨hx.1, hx.2.trans (by nlinarith)⟩ + · exact ⟨by nlinarith, hx.2⟩ + rw [h_union_set] + exact h_deriv_left.union h_deriv_right + · have h_gt : (k+1)*T < t := by nlinarith + have ht' : t ∈ Icc ((k+1)*T) ((k+2)*T) := ⟨h_ge, ht.2⟩ + let x_k := γ_seq k ((k+1)*T) + have h_next_deriv : HasDerivWithinAt ((h_local_step x_k ((k+1)*T)).choose) + (f (((h_local_step x_k ((k+1)*T)).choose) t)) (Icc ((k+1)*T) ((k+2)*T)) t := + (h_local_step x_k ((k+1)*T)).choose_spec.2 t ht' + have h_eq_on : ∀ y ∈ Icc ((k+1)*T) ((k+2)*T), + γ_seq (k+1) y = ((h_local_step x_k ((k+1)*T)).choose) y := by + intro y hy; dsimp [γ_seq] + have : ¬ y ≤ (k+1)*T := by intro hle; nlinarith + simp [this] + have hx : γ_seq (k+1) t = ((h_local_step x_k ((k+1)*T)).choose) t := h_eq_on t ht' + have h_f : f (γ_seq (k+1) t) = f (((h_local_step x_k ((k+1)*T)).choose) t) := by rw [hx] + have h_sub : Icc ((k+1)*T) ((k+2)*T) ⊆ Icc 0 ((k+2)*T) := by + intro x hx; exact ⟨by nlinarith, hx.2⟩ + have h_deriv_next' : HasDerivWithinAt ((h_local_step x_k ((k+1)*T)).choose) + (f (γ_seq (k+1) t)) (Icc ((k+1)*T) ((k+2)*T)) t := by + simpa [h_f] using h_next_deriv + exact (h_deriv_next'.congr (fun y hy => (h_eq_on y hy).symm) hx.symm).mono h_sub + + -- Package the forward solution as a reusable lemma + have exists_forward_solution (f' : E → E) (hf_lip' : ∀ x y, ‖f' x - f' y‖ ≤ K * ‖x - y‖) : + ∃ (γ_fwd : ℝ → E), γ_fwd 0 = x₀ ∧ (∀ t > 0, HasDerivAt γ_fwd (f' (γ_fwd t)) t) ∧ + HasDerivWithinAt γ_fwd (f' x₀) (Set.Ici 0) 0 := by + have h_local_f' : ∀ (x : E) (t₀ : ℝ), ∃ γ : ℝ → E, γ t₀ = x ∧ + ∀ t ∈ Icc t₀ (t₀ + τ), HasDerivWithinAt γ (f' (γ t)) (Icc t₀ (t₀ + τ)) t := + h_local_step' f' hf_lip' + let γ_seq_fwd : ℕ → ℝ → E := + Nat.rec (h_local_f' x₀ 0).choose (fun k prev => + let x_k := prev ((k+1)*T) + let next := (h_local_f' x_k ((k+1)*T)).choose + fun t => if h : t ≤ (k+1)*T then prev t else next t) + have h_seq0_fwd : γ_seq_fwd 0 0 = x₀ := (h_local_f' x₀ 0).choose_spec.1 + have h_seq_cons_fwd : ∀ (k : ℕ) (t : ℝ), t ≤ (k+1)*T → γ_seq_fwd (k+1) t = γ_seq_fwd k t := by + intro k t ht; dsimp [γ_seq_fwd]; simp [ht] + have h_seq_mono_fwd : ∀ (s : ℝ) (n m : ℕ), n ≤ m → s ≤ (n+1 : ℝ)*T → γ_seq_fwd m s = γ_seq_fwd n s := by + intro s n m hnm hs + induction' hnm with m hnm ih + · rfl + · have hsm : s ≤ (m+1 : ℝ)*T := by + have : (n : ℝ) ≤ m := by exact_mod_cast hnm.le + nlinarith + rw [h_seq_cons_fwd m s hsm, ih] + have h_deriv_seq_fwd : ∀ (k : ℕ) (t : ℝ), t ∈ Icc 0 ((k+1)*T) → + HasDerivWithinAt (γ_seq_fwd k) (f' ((γ_seq_fwd k) t)) (Icc 0 ((k+1)*T)) t := by + intro k + induction' k with k ih + · intro t ht + exact ((h_local_f' x₀ 0).choose_spec.2 t ht) + · intro t ht + dsimp [γ_seq_fwd] + rcases lt_or_ge t ((k+1)*T) with (h_lt | h_ge) + · have ht' : t ∈ Icc 0 ((k+1)*T) := ⟨ht.1, h_lt.le⟩ + have h_deriv_k : HasDerivWithinAt (γ_seq_fwd k) (f' ((γ_seq_fwd k) t)) (Icc 0 ((k+1)*T)) t := + ih t ht' + have h_sub : Icc 0 ((k+1)*T) ⊆ Icc 0 ((k+2)*T) := by + intro x hx; exact ⟨hx.1, hx.2.trans (by nlinarith)⟩ + have h_eq_on : ∀ y ∈ Icc 0 ((k+1)*T), (γ_seq_fwd (k+1)) y = (γ_seq_fwd k) y := + fun y hy => h_seq_cons_fwd k y hy.2 + have hx : (γ_seq_fwd (k+1)) t = (γ_seq_fwd k) t := h_eq_on t ht' + have h_f' : f' ((γ_seq_fwd (k+1)) t) = f' ((γ_seq_fwd k) t) := by rw [hx] + have h_deriv_k' : HasDerivWithinAt (γ_seq_fwd k) (f' ((γ_seq_fwd (k+1)) t)) (Icc 0 ((k+1)*T)) t := by + simpa [h_f'] using h_deriv_k + exact ((h_deriv_k'.congr (fun y hy => (h_eq_on y hy).symm) hx.symm).mono h_sub) + · by_cases h_eq : t = (k+1)*T + · subst h_eq + let x_k := γ_seq_fwd k ((k+1)*T) + have hx_k_eq : γ_seq_fwd (k+1) ((k+1)*T) = x_k := by + symm; apply h_seq_cons_fwd k ((k+1)*T) (le_refl _) + have h_left_pt : (k+1)*T ∈ Icc 0 ((k+1)*T) := by + constructor <;> nlinarith + have h_deriv_left_γ : HasDerivWithinAt (γ_seq_fwd k) (f' (x_k)) (Icc 0 ((k+1)*T)) ((k+1)*T) := by + simpa using ih ((k+1)*T) h_left_pt + have h_left_eq_on : ∀ y ∈ Icc 0 ((k+1)*T), γ_seq_fwd (k+1) y = γ_seq_fwd k y := + fun y hy => h_seq_cons_fwd k y hy.2 + have h_deriv_left : HasDerivWithinAt (γ_seq_fwd (k+1)) (f' (γ_seq_fwd (k+1) ((k+1)*T))) + (Icc 0 ((k+1)*T)) ((k+1)*T) := by + apply h_deriv_left_γ.congr (fun y hy => (h_left_eq_on y hy).symm) ? + simpa [hx_k_eq] using (h_left_eq_on ((k+1)*T) h_left_pt).symm + have h_right_pt : (k+1)*T ∈ Icc ((k+1)*T) ((k+2)*T) := by + constructor <;> nlinarith + have h_next_deriv : ∀ t ∈ Icc ((k+1)*T) ((k+2)*T), + HasDerivWithinAt ((h_local_f' x_k ((k+1)*T)).choose) (f' (((h_local_f' x_k ((k+1)*T)).choose) t)) (Icc ((k+1)*T) ((k+2)*T)) t := + (h_local_f' x_k ((k+1)*T)).choose_spec.2 + have h_deriv_right_next : HasDerivWithinAt ((h_local_f' x_k ((k+1)*T)).choose) + (f' (γ_seq_fwd (k+1) ((k+1)*T))) (Icc ((k+1)*T) ((k+2)*T)) ((k+1)*T) := by + simpa [hx_k_eq] using h_next_deriv ((k+1)*T) h_right_pt + have h_right_eq_on : ∀ y ∈ Icc ((k+1)*T) ((k+2)*T), + γ_seq_fwd (k+1) y = ((h_local_f' x_k ((k+1)*T)).choose) y := by + intro y hy; dsimp [γ_seq_fwd] + have : ¬ y ≤ (k+1)*T := by intro hle; nlinarith + simp [this] + have h_deriv_right : HasDerivWithinAt (γ_seq_fwd (k+1)) (f' (γ_seq_fwd (k+1) ((k+1)*T))) + (Icc ((k+1)*T) ((k+2)*T)) ((k+1)*T) := + h_deriv_right_next.congr (fun y hy => (h_right_eq_on y hy).symm) + ((h_right_eq_on ((k+1)*T) h_right_pt).symm) + have h_union_set : Icc 0 ((k+2)*T) = + Icc 0 ((k+1)*T) ∪ Icc ((k+1)*T) ((k+2)*T) := by + ext x; constructor + · intro hx; rcases le_total x ((k+1)*T) with (hx' | hx') + · exact Or.inl ⟨hx.1, hx'⟩ + · exact Or.inr ⟨hx', hx.2⟩ + · intro (hx | hx) + · exact ⟨hx.1, hx.2.trans (by nlinarith)⟩ + · exact ⟨by nlinarith, hx.2⟩ + rw [h_union_set] + exact h_deriv_left.union h_deriv_right + · have h_gt : (k+1)*T < t := by nlinarith + have ht' : t ∈ Icc ((k+1)*T) ((k+2)*T) := ⟨h_ge, ht.2⟩ + let x_k := γ_seq_fwd k ((k+1)*T) + have h_next_deriv : HasDerivWithinAt ((h_local_f' x_k ((k+1)*T)).choose) + (f' (((h_local_f' x_k ((k+1)*T)).choose) t)) (Icc ((k+1)*T) ((k+2)*T)) t := + (h_local_f' x_k ((k+1)*T)).choose_spec.2 t ht' + have h_eq_on : ∀ y ∈ Icc ((k+1)*T) ((k+2)*T), + γ_seq_fwd (k+1) y = ((h_local_f' x_k ((k+1)*T)).choose) y := by + intro y hy; dsimp [γ_seq_fwd] + have : ¬ y ≤ (k+1)*T := by intro hle; nlinarith + simp [this] + have hx : γ_seq_fwd (k+1) t = ((h_local_f' x_k ((k+1)*T)).choose) t := h_eq_on t ht' + have h_f' : f' (γ_seq_fwd (k+1) t) = f' (((h_local_f' x_k ((k+1)*T)).choose) t) := by rw [hx] + have h_sub : Icc ((k+1)*T) ((k+2)*T) ⊆ Icc 0 ((k+2)*T) := by + intro x hx; exact ⟨by nlinarith, hx.2⟩ + have h_deriv_next' : HasDerivWithinAt ((h_local_f' x_k ((k+1)*T)).choose) + (f' (γ_seq_fwd (k+1) t)) (Icc ((k+1)*T) ((k+2)*T)) t := by + simpa [h_f'] using h_next_deriv + exact (h_deriv_next'.congr (fun y hy => (h_eq_on y hy).symm) hx.symm).mono h_sub + have h_exists_N (t : ℝ) (ht : t ≥ 0) : ∃ N : ℕ, t < (N : ℝ) * T := by + obtain ⟨n, hn⟩ : ∃ n : ℕ, t / T < (n : ℝ) := by + obtain ⟨n, hn⟩ := exists_nat_gt (t / T) + exact ⟨n, hn⟩ + refine ⟨n, ?_⟩ + nlinarith + let γ_fwd : ℝ → E := fun t => + if ht : t ≥ 0 then + let N := (h_exists_N t ht).choose + γ_seq_fwd N t + else + x₀ + have hγ_fwd_eq_seq (t : ℝ) (ht : t ≥ 0) (N : ℕ) (hN : t < (N : ℝ) * T) : γ_fwd t = γ_seq_fwd N t := by + dsimp [γ_fwd] + simp [ht] + set M := (h_exists_N t ht).choose with hM + have hM_bound : t < (M : ℝ) * T := (h_exists_N t ht).choose_spec + have hM_N : M ≤ N := by + apply Nat.find_min' (h_exists_N t ht) + exact hN + have hseq_eq : γ_seq_fwd N t = γ_seq_fwd M t := by + symm; apply h_seq_mono_fwd t M N hM_N hM_bound.le + calc γ_seq_fwd M t = γ_seq_fwd M t := rfl + _ = γ_seq_fwd N t := by symm; exact hseq_eq + have hγ_fwd0 : γ_fwd 0 = x₀ := by + dsimp [γ_fwd]; simp [h_seq0_fwd] + have hγ_fwd_at_pos : ∀ t > 0, HasDerivAt γ_fwd (f' (γ_fwd t)) t := by + intro t ht + have ht0 : t ≥ 0 := ht.le + rcases h_exists_N t ht0 with ⟨N, hN⟩ + have h_strict : t < (N+1 : ℕ)*T := by + have hNpos_real : (N : ℝ) ≥ 0 := by exact_mod_cast Nat.zero_le _ + calc + t < (N : ℝ) * T := hN + _ < (N+1 : ℝ) * T := by nlinarith + set ε := min (t/2) (((N+1 : ℕ)*T - t)/2) with hε + have hε_pos : 0 < ε := by + refine lt_min_iff.mpr ⟨by nlinarith, by nlinarith⟩ + have h_ball_eq : ∀ s ∈ Set.Ioo (t - ε) (t + ε), γ_fwd s = γ_seq_fwd (N+1) s := by + intro s hs + have hs_pos : s ≥ 0 := by + have : t - ε > 0 := by + have : t/2 > 0 := by nlinarith + nlinarith + nlinarith + have hs_bound : s < ((N+1 : ℕ) : ℝ) * T := by nlinarith + apply hγ_fwd_eq_seq s hs_pos (N+1) hs_bound + have h_ball_eventually : γ_fwd =ᶠ[𝓝 t] γ_seq_fwd (N+1) := by + apply Filter.eventually_of_mem (Set.Ioo (t - ε) (t + ε)) + · exact h_ball_eq + · apply IsOpen.mem_nhds isOpen_Ioo + constructor <;> nlinarith + have h_deriv_seq_at : HasDerivAt (γ_seq_fwd (N+1)) (f' (γ_seq_fwd (N+1) t)) t := by + have hmem : Icc 0 (((N+1 : ℕ)+1)*T) ∈ 𝓝 t := by + have h_open : Ioo 0 (((N+1 : ℕ)+1)*T) ∈ 𝓝 t := + IsOpen.mem_nhds isOpen_Ioo ⟨by nlinarith, by nlinarith⟩ + exact mem_of_superset h_open Set.Ioo_subset_Icc_self + have h_mem_icc : t ∈ Icc 0 (((N+1 : ℕ)+1)*T) := by + constructor <;> nlinarith + have h_within : HasDerivWithinAt (γ_seq_fwd (N+1)) (f' (γ_seq_fwd (N+1) t)) + (Icc 0 (((N+1 : ℕ)+1)*T)) t := + h_deriv_seq_fwd (N+1) t h_mem_icc + exact h_within.hasDerivAt hmem + have h_f_eq : f' (γ_seq_fwd (N+1) t) = f' (γ_fwd t) := by + have h_val : γ_fwd t = γ_seq_fwd (N+1) t := h_ball_eq t (by + constructor <;> nlinarith) + simp [h_val] + have h_temp : HasDerivAt γ_fwd (f' (γ_seq_fwd (N+1) t)) t := + h_deriv_seq_at.congr_of_eventuallyEq h_ball_eventually + simpa [h_f_eq] using h_temp + have h_right_deriv0_fwd : HasDerivWithinAt γ_fwd (f' x₀) (Set.Ici 0) 0 := by + have h0 : HasDerivWithinAt (γ_seq_fwd 0) (f' x₀) (Icc 0 T) 0 := by + have hmem : (0 : ℝ) ∈ Icc 0 T := by + refine ⟨by norm_num, ?_⟩; nlinarith + simpa using h_deriv_seq_fwd 0 (0 : ℝ) hmem + have h_eq_on_icc : ∀ y ∈ Icc 0 T, γ_fwd y = γ_seq_fwd 0 y := by + intro y hy + rcases hy.2.eq_or_lt with (rfl | h_lt) + · dsimp [γ_fwd]; simp [h_seq0_fwd] + · have hy_nonneg : y ≥ 0 := hy.1 + have hy_bound : y < (0 : ℝ) * T := by nlinarith + apply hγ_fwd_eq_seq y hy_nonneg 0 ?_ + nlinarith + have h_γ_deriv_icc : HasDerivWithinAt γ_fwd (f' x₀) (Icc 0 T) 0 := + h0.congr (fun y hy => (h_eq_on_icc y hy).symm) (h_eq_on_icc 0 ⟨by norm_num, by nlinarith⟩).symm + exact h_γ_deriv_icc.mono (by intro y hy; exact hy.1) + exact ⟨γ_fwd, hγ_fwd0, hγ_fwd_at_pos, h_right_deriv0_fwd⟩ + + -- Obtain forward solution for f (on [0,∞)) + obtain ⟨γ_pos, hγ_pos0, hγ_pos_at_pos, hγ_pos_right0⟩ := + exists_forward_solution f h_lipschitz + + -- Build backward solution via time reversal + let f_neg (x : E) : E := -f x + have h_f_neg_lip : ∀ x y, ‖f_neg x - f_neg y‖ ≤ K * ‖x - y‖ := by + intro x y; simpa [f_neg] using h_lipschitz x y + obtain ⟨γ_neg_fwd, hγ_neg_fwd0, hγ_neg_fwd_at_pos, hγ_neg_fwd_right0⟩ := + exists_forward_solution f_neg h_f_neg_lip + + -- Splice positive and negative into full solution + let γ_full : ℝ → E := fun t => + if h : t ≥ 0 then γ_pos t else γ_neg_fwd (-t) + + have hγ_full0 : γ_full 0 = x₀ := by + simp [γ_full, hγ_pos0, hγ_neg_fwd0] + + have hγ_full_at_pos : ∀ t > 0, HasDerivAt γ_full (f (γ_full t)) t := by + intro t ht + have : γ_full t = γ_pos t := by simp [γ_full, ht.le] + have h_pos : HasDerivAt γ_pos (f (γ_pos t)) t := hγ_pos_at_pos t ht + simpa [this] using h_pos + + have hγ_full_at_neg : ∀ t < 0, HasDerivAt γ_full (f (γ_full t)) t := by + intro t ht + have h_neg_nonpos : t ≤ 0 := ht.le + have hs_pos : -t > 0 := by linarith + have h_val : γ_full t = γ_neg_fwd (-t) := by + simp [γ_full, h_neg_nonpos] + have h_inner_deriv : HasDerivAt (fun u : ℝ => -u) (-1) t := by + simpa using hasDerivAt_neg t + have h_outer_deriv : HasDerivAt γ_neg_fwd (f_neg (γ_neg_fwd (-t))) (-t) := + hγ_neg_fwd_at_pos (-t) hs_pos + have h_comp : HasDerivAt (γ_neg_fwd ∘ (fun u : ℝ => -u)) ((-1 : ℝ) • f_neg (γ_neg_fwd (-t))) t := + HasDerivAt.comp (x := t) h_inner_deriv h_outer_deriv + have h_simp : (-1 : ℝ) • f_neg (γ_neg_fwd (-t)) = f (γ_neg_fwd (-t)) := by + simp [f_neg] + have h_fun : (γ_neg_fwd ∘ (fun u : ℝ => -u)) = γ_full := by + ext u; simp [γ_full] + rw [h_fun] at h_comp + simpa [h_val, h_simp] using h_comp + + have h_deriv_at_0 : HasDerivAt γ_full (f x₀) 0 := by + rw [hasDerivAt_iff_tendsto_slope_left_right] + have h_right_slope : Filter.Tendsto (fun y => (γ_full y - γ_full 0) / (y - 0)) + (𝓝[≥] 0) (𝓝 (f x₀)) := by + have : γ_full =ᶠ[𝓝[≥] 0] γ_pos := by + apply Filter.eventually_of_mem (Set.Ici 0) + · intro y hy; simp [γ_full, hy] + · exact IsOpen.mem_nhdsWithin isOpen_Ici (by norm_num : (0 : ℝ) ∈ Set.Ici 0) + have h_slope_pos : Filter.Tendsto (fun y => (γ_pos y - γ_pos 0) / (y - 0)) + (𝓝[≥] 0) (𝓝 (f x₀)) := + hγ_pos_right0.tendsto_slope + simpa [hγ_pos0, hγ_full0] using h_slope_pos.congr this + have h_left_slope : Filter.Tendsto (fun y => (γ_full y - γ_full 0) / (y - 0)) + (𝓝[≤] 0) (𝓝 (f x₀)) := by + have h_slope_neg : Filter.Tendsto (fun s : ℝ => (γ_neg_fwd s - x₀) / s) + (𝓝[≥] 0) (𝓝 (f_neg x₀)) := + hγ_neg_fwd_right0.tendsto_slope + have h_neg_map : Filter.Tendsto (fun y : ℝ => -y) (𝓝[≤] 0) (𝓝[≥] 0) := by + have h_nhds_neg : Filter.Tendsto (fun y : ℝ => -y) (𝓝 0) (𝓝 0) := + continuous_neg.tendsto 0 + have h_principal : Filter.Tendsto (fun y : ℝ => -y) (𝓟 (Set.Iic 0)) (𝓟 (Set.Ici 0)) := + Set.MapsTo.tendsto (fun y hy => by simpa using hy) + exact h_nhds_neg.inf h_principal + have h_comp_raw : Filter.Tendsto (fun y : ℝ => (γ_neg_fwd (-y) - x₀) / (-y)) + (𝓝[≤] 0) (𝓝 (f_neg x₀)) := + h_slope_neg.comp h_neg_map + have h_eq_on : (fun y : ℝ => (γ_full y - x₀) / y) =ᶠ[𝓝[≤] 0] + (fun y : ℝ => -((γ_neg_fwd (-y) - x₀) / (-y))) := by + have h_mem : Set.Iic (0 : ℝ) ∈ 𝓝[≤] 0 := + Filter.mem_inf_of_inter Filter.univ_mem (Filter.mem_principal_self _) + refine Filter.eventually_of_mem h_mem (fun y hy => ?_) + calc + (γ_full y - x₀) / y = (γ_neg_fwd (-y) - x₀) / y := by simp [γ_full, hy] + _ = -((γ_neg_fwd (-y) - x₀) / (-y)) := by + by_cases hy0 : y = 0 + · subst hy0; simp + · field_simp [hy0]; ring + have h_neg_target : Filter.Tendsto (fun y : ℝ => -((γ_neg_fwd (-y) - x₀) / (-y))) + (𝓝[≤] 0) (𝓝 (-(f_neg x₀))) := + (continuous_neg.tendsto (f_neg x₀)).comp h_comp_raw + have h_target : -(f_neg x₀) = f x₀ := by simp [f_neg] + have h_result : Filter.Tendsto (fun y : ℝ => (γ_full y - x₀) / y) (𝓝[≤] 0) (𝓝 (f x₀)) := by + simpa [hγ_full0, h_target] using h_neg_target.congr h_eq_on.symm + simpa [sub_self] using h_result + exact ⟨h_right_slope, h_left_slope⟩ + + have h_solution : IsSolutionND f x₀ γ_full := by + refine ⟨hγ_full0, fun t => ?_⟩ + by_cases ht : t > 0 + · exact hγ_full_at_pos t ht + · by_cases ht' : t < 0 + · exact hγ_full_at_neg t ht' + · have ht0 : t = 0 := by nlinarith + subst ht0 + have hf0 : f (γ_full 0) = f x₀ := by simp [hγ_full0] + simp [h_deriv_at_0, hf0] + refine ⟨γ_full, h_solution⟩ exact exists_unique_of_exists_of_unique hex hunique - --- ============================================================================ --- Multivariable Differential Calculus on Phase Space --- ============================================================================ - -/-- Phase space for a system with n degrees of freedom: pairs (q, p) of -position and momentum vectors. -/ -def PhaseSpace (n : ℕ) := (Fin n → ℝ) × (Fin n → ℝ) - -/-- Replace the q-coordinate at index i with value t, keeping all other -coordinates fixed. -/ -def replaceQ {n : ℕ} (x : PhaseSpace n) (i : Fin n) (t : ℝ) : PhaseSpace n := - ⟨fun j => if j = i then t else x.1 j, x.2⟩ - -/-- Replace the p-coordinate at index i with value t, keeping all other -coordinates fixed. -/ -def replaceP {n : ℕ} (x : PhaseSpace n) (i : Fin n) (t : ℝ) : PhaseSpace n := - ⟨x.1, fun j => if j = i then t else x.2 j⟩ - -/-- H has partial derivative d with respect to q_i at point x. -/ -def HasPartialDerivAtQ {n : ℕ} (H : PhaseSpace n → ℝ) (i : Fin n) (d : ℝ) (x : PhaseSpace n) : Prop := - HasDerivAt (fun t : ℝ => H (replaceQ x i t)) d (x.1 i) - -/-- H has partial derivative d with respect to p_i at point x. -/ -def HasPartialDerivAtP {n : ℕ} (H : PhaseSpace n → ℝ) (i : Fin n) (d : ℝ) (x : PhaseSpace n) : Prop := - HasDerivAt (fun t : ℝ => H (replaceP x i t)) d (x.2 i) - -/-- The gradient of H with respect to the q-coordinates. -If the partial derivative does not exist, we arbitrarily set it to 0. -/ -noncomputable def gradQ {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) : Fin n → ℝ := - fun i => if h : ∃ d, HasPartialDerivAtQ H i d x then h.choose else 0 - -/-- The gradient of H with respect to the p-coordinates. -If the partial derivative does not exist, we arbitrarily set it to 0. -/ -noncomputable def gradP {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) : Fin n → ℝ := - fun i => if h : ∃ d, HasPartialDerivAtP H i d x then h.choose else 0 - --- ============================================================================ --- Symplectic Structure --- ============================================================================ - -/-- The canonical symplectic form on phase space: -ω((q,p), (q',p')) = Σ_i (p_i q'_i - q_i p'_i). -/ -def ω {n : ℕ} (u v : PhaseSpace n) : ℝ := - ∑ i : Fin n, (u.2 i * v.1 i - u.1 i * v.2 i) - -/-- The symplectic form is bilinear (linear in each argument separately). -We state this as two theorems. -/ -theorem ω_linear_left {n : ℕ} (u₁ u₂ v : PhaseSpace n) (c : ℝ) : - ω (c • u₁ + u₂) v = c * ω u₁ v + ω u₂ v := by - unfold ω - simp [Finset.sum_add_distrib, Finset.mul_sum, add_mul, mul_add, mul_sub] - rw [← Finset.sum_add_distrib] - apply Finset.sum_congr rfl - intro i _ - ring - -/-- The symplectic form is antisymmetric: ω(u, v) = -ω(v, u). -/ -theorem ω_antisymmetric {n : ℕ} (u v : PhaseSpace n) : - ω u v = -ω v u := by - unfold ω - rw [← Finset.sum_neg_distrib] - apply Finset.sum_congr rfl - intro i _ - ring - -/-- The symplectic form is non-degenerate: if ω(u, v) = 0 for all v, then u = 0. -/ -theorem ω_nondegenerate {n : ℕ} (u : PhaseSpace n) - (h : ∀ v : PhaseSpace n, ω u v = 0) : u = 0 := by - have hq : u.1 = 0 := by - funext j - let v : PhaseSpace n := ⟨0, fun i => if i = j then 1 else 0⟩ - have hω := h v - simp [ω, v] at hω - linarith - have hp : u.2 = 0 := by - funext j - let v : PhaseSpace n := ⟨fun i => if i = j then 1 else 0, 0⟩ - have hω := h v - simp [ω, v] at hω - linarith - exact Prod.ext hq hp - --- ============================================================================ --- Hamiltonian Vector Field --- ============================================================================ - -/-- The Hamiltonian vector field X_H associated to a Hamiltonian function H. -In coordinates: X_H = (∂H/∂p, -∂H/∂q), which is equivalent to -ω(X_H, ·) = dH via the symplectic musical isomorphism. -/ -noncomputable def HamiltonianVectorField {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) : PhaseSpace n := - ⟨gradP H x, -gradQ H x⟩ - -/-- Hamilton's equations in vector field form: a trajectory γ satisfies -ẋ = X_H(x), i.e., the time derivative of γ at t equals the Hamiltonian -vector field evaluated at γ(t). -/ -def IsHamiltonianTrajectory {n : ℕ} (H : PhaseSpace n → ℝ) (γ : ℝ → PhaseSpace n) : Prop := - ∀ t, HasDerivAt γ (HamiltonianVectorField H (γ t)) t - -/-- Hamilton's equations in coordinates: for each degree of freedom i, -dq_i/dt = ∂H/∂p_i and dp_i/dt = -∂H/∂q_i. - -This follows from the product decomposition of HasDerivAt for -functions into product types (Mathlib.Analysis.Calculus.Deriv.Prod) -combined with component-wise extraction via continuous linear evaluation. -/ -theorem hamiltonian_eqs_coordinates {n : ℕ} (H : PhaseSpace n → ℝ) (γ : ℝ → PhaseSpace n) - (h_traj : IsHamiltonianTrajectory H γ) (t : ℝ) : - ∀ i : Fin n, HasDerivAt (fun s => (γ s).1 i) (gradP H (γ t) i) t ∧ - HasDerivAt (fun s => (γ s).2 i) (-gradQ H (γ t) i) t := by - intro i - have h_deriv := h_traj t - have h_fst : HasDerivAt (fun s => (γ s).1) (HamiltonianVectorField H (γ t)).1 t := by - simpa [HamiltonianVectorField] using HasDerivAt.fst h_deriv - have h_snd : HasDerivAt (fun s => (γ s).2) (HamiltonianVectorField H (γ t)).2 t := by - simpa [HamiltonianVectorField] using HasDerivAt.snd h_deriv - constructor - · -- q_i derivative: compose the first-component derivative with evaluation at i - let eval_i : (Fin n → ℝ) →L[ℝ] ℝ := - { toFun := fun f => f i, - map_add' := by intros; rfl, - map_smul' := by intros; rfl, - cont := by continuity } - have heval : HasFDerivAt eval_i eval_i ((γ t).1) := eval_i.hasFDerivAt - have hcomp := HasFDerivAt.comp_hasDerivAt t heval h_fst - simpa [HamiltonianVectorField, eval_i] using hcomp - · -- p_i derivative: compose the second-component derivative with evaluation at i - let eval_i : (Fin n → ℝ) →L[ℝ] ℝ := - { toFun := fun f => f i, - map_add' := by intros; rfl, - map_smul' := by intros; rfl, - cont := by continuity } - have heval : HasFDerivAt eval_i eval_i ((γ t).2) := eval_i.hasFDerivAt - have hcomp := HasFDerivAt.comp_hasDerivAt t heval h_snd - simpa [HamiltonianVectorField, eval_i] using hcomp - --- ============================================================================ --- Example: The 1D Harmonic Oscillator --- ============================================================================ - -/-- The Hamiltonian for a 1D harmonic oscillator with mass m and spring constant k: -H(q, p) = p²/(2m) + kq²/2. -/ -def HarmonicOscillatorH (m k : ℝ) (hm : m > 0) (hk : k > 0) : PhaseSpace 1 → ℝ := - fun x => (x.2 0)^2 / (2 * m) + k * (x.1 0)^2 / 2 - -/-- The partial derivative of the harmonic oscillator Hamiltonian with respect to p. -/ -theorem harmonic_oscillator_partial_p (m k : ℝ) (hm : m > 0) (hk : k > 0) (q p : ℝ) : - HasDerivAt (fun t : ℝ => t^2 / (2 * m) + k * q^2 / 2) (p / m) p := by - have h1 : HasDerivAt (fun t => t^2) (2 * p) p := hasDerivAt_pow 2 p - have h2 : HasDerivAt (fun t => t^2 / (2 * m)) ((2 * p) / (2 * m)) p := by - simpa using HasDerivAt.mul_const (c := 1 / (2 * m)) h1 - have h3 : HasDerivAt (fun _ => k * q^2 / 2) 0 p := hasDerivAt_const p (k * q^2 / 2) - convert HasDerivAt.add h2 h3 using 1 - field_simp [ne_of_gt hm] - ring - -/-- The partial derivative of the harmonic oscillator Hamiltonian with respect to q. -/ -theorem harmonic_oscillator_partial_q (m k : ℝ) (hm : m > 0) (hk : k > 0) (q p : ℝ) : - HasDerivAt (fun t : ℝ => p^2 / (2 * m) + k * t^2 / 2) (k * q) q := by - have h1 : HasDerivAt (fun _ => p^2 / (2 * m)) 0 q := hasDerivAt_const q (p^2 / (2 * m)) - have h2 : HasDerivAt (fun t => t^2) (2 * q) q := hasDerivAt_pow 2 q - have h3 : HasDerivAt (fun t => k * t^2 / 2) (k * (2 * q) / 2) q := by - simpa using HasDerivAt.mul_const (c := k / 2) h2 - convert HasDerivAt.add h1 h3 using 1 - field_simp [ne_of_gt hm] - ring - -/-- The q-gradient of the harmonic oscillator Hamiltonian equals k·q. -/ -lemma gradQ_harmonic_oscillator (m k : ℝ) (hm : m > 0) (hk : k > 0) (x : PhaseSpace 1) : - gradQ (HarmonicOscillatorH m k hm hk) x 0 = k * x.1 0 := by - have hex : ∃ d, HasPartialDerivAtQ (HarmonicOscillatorH m k hm hk) 0 d x := by - use k * x.1 0 - simpa [HasPartialDerivAtQ, replaceQ, HarmonicOscillatorH] using - harmonic_oscillator_partial_q m k hm hk (x.1 0) (x.2 0) - have h1 : HasPartialDerivAtQ (HarmonicOscillatorH m k hm hk) 0 (gradQ (HarmonicOscillatorH m k hm hk) x 0) x := by - simpa [gradQ, hex] using Classical.choose_spec hex - have h2 : HasPartialDerivAtQ (HarmonicOscillatorH m k hm hk) 0 (k * x.1 0) x := by - simpa [HasPartialDerivAtQ, replaceQ, HarmonicOscillatorH] using - harmonic_oscillator_partial_q m k hm hk (x.1 0) (x.2 0) - exact HasDerivAt.unique h1 h2 - -/-- The p-gradient of the harmonic oscillator Hamiltonian equals p/m. -/ -lemma gradP_harmonic_oscillator (m k : ℝ) (hm : m > 0) (hk : k > 0) (x : PhaseSpace 1) : - gradP (HarmonicOscillatorH m k hm hk) x 0 = x.2 0 / m := by - have hex : ∃ d, HasPartialDerivAtP (HarmonicOscillatorH m k hm hk) 0 d x := by - use x.2 0 / m - simpa [HasPartialDerivAtP, replaceP, HarmonicOscillatorH] using - harmonic_oscillator_partial_p m k hm hk (x.1 0) (x.2 0) - have h1 : HasPartialDerivAtP (HarmonicOscillatorH m k hm hk) 0 (gradP (HarmonicOscillatorH m k hm hk) x 0) x := by - simpa [gradP, hex] using Classical.choose_spec hex - have h2 : HasPartialDerivAtP (HarmonicOscillatorH m k hm hk) 0 (x.2 0 / m) x := by - simpa [HasPartialDerivAtP, replaceP, HarmonicOscillatorH] using - harmonic_oscillator_partial_p m k hm hk (x.1 0) (x.2 0) - exact HasDerivAt.unique h1 h2 - -/-- The Hamiltonian equations for the 1D harmonic oscillator are: -ẋq = p/m and ẋp = -kq. -/ -theorem harmonic_oscillator_equations (m k : ℝ) (hm : m > 0) (hk : k > 0) (γ : ℝ → PhaseSpace 1) - (h_traj : IsHamiltonianTrajectory (HarmonicOscillatorH m k hm hk) γ) (t : ℝ) : - HasDerivAt (fun s => (γ s).1 0) ((γ t).2 0 / m) t ∧ - HasDerivAt (fun s => (γ s).2 0) (-k * (γ t).1 0) t := by - have hcoords := hamiltonian_eqs_coordinates (HarmonicOscillatorH m k hm hk) γ h_traj t 0 - have hq : gradP (HarmonicOscillatorH m k hm hk) (γ t) 0 = (γ t).2 0 / m := gradP_harmonic_oscillator m k hm hk (γ t) - have hp : gradQ (HarmonicOscillatorH m k hm hk) (γ t) 0 = k * (γ t).1 0 := gradQ_harmonic_oscillator m k hm hk (γ t) - constructor - · rw [hq] at hcoords - exact hcoords.1 - · rw [hp] at hcoords - have h_neg : -gradQ (HarmonicOscillatorH m k hm hk) (γ t) 0 = -k * (γ t).1 0 := by - rw [hp] - ring - rwa [h_neg] at hcoords - --- ============================================================================ --- Energy Conservation --- ============================================================================ - -/-- Any continuous linear functional on PhaseSpace n can be evaluated using -the standard basis vectors. In coordinates: -L(v_q, v_p) = Σ_i (v_q i · L(e_{q_i}) + v_p i · L(e_{p_i})). -/ -lemma clm_on_phaseSpace {n : ℕ} (L : PhaseSpace n →L[ℝ] ℝ) (v : PhaseSpace n) : - L v = ∑ i : Fin n, (v.1 i * L (⟨fun j => if j = i then 1 else 0, 0⟩) - + v.2 i * L (⟨0, fun j => if j = i then 1 else 0⟩)) := by - have h1 : v.1 = ∑ i : Fin n, v.1 i • (fun j => if j = i then (1 : ℝ) else 0) := by - funext j - simp [Finset.sum_apply, ite_mul, zero_mul, mul_ite, mul_zero] - rw [Finset.sum_ite_eq'] - simp - have h2 : v.2 = ∑ i : Fin n, v.2 i • (fun j => if j = i then (1 : ℝ) else 0) := by - funext j - simp [Finset.sum_apply, ite_mul, zero_mul, mul_ite, mul_zero] - rw [Finset.sum_ite_eq'] - simp - have h3 : v = (∑ i : Fin n, v.1 i • (⟨fun j => if j = i then 1 else 0, 0⟩ : PhaseSpace n)) - + (∑ i : Fin n, v.2 i • (⟨0, fun j => if j = i then 1 else 0⟩ : PhaseSpace n)) := by - rw [Prod.ext_iff] - constructor - · rw [h1] - simp [Finset.sum_apply, smul_eq_mul] - · rw [h2] - simp [Finset.sum_apply, smul_eq_mul] - rw [h3] - simp [map_add, map_sum, map_smul] - rw [← Finset.sum_add_distrib] - apply Finset.sum_congr rfl - intro i _ - ring - -/-- Energy conservation: along a Hamiltonian trajectory, the Hamiltonian H(γ(t)) -is constant in time, i.e., its time derivative is 0. - -The proof uses the chain rule for Fréchet derivatives: -d/dt H(γ(t)) = D H(γ(t)) [γ'(t)] = D H(γ(t)) [X_H(γ(t))]. - -We require as a hypothesis that H has a Fréchet derivative at γ(t) which agrees -with the partial derivatives on the standard basis vectors. This is equivalent -to differentiability of H. The Hamiltonian structure then guarantees the -derivative vanishes by antisymmetry of the symplectic form. -/ -theorem energy_conservation {n : ℕ} (H : PhaseSpace n → ℝ) (γ : ℝ → PhaseSpace n) - (h_traj : IsHamiltonianTrajectory H γ) (t : ℝ) - (h_diff : ∃ L : PhaseSpace n →L[ℝ] ℝ, HasFDerivAt H L (γ t) ∧ - (∀ i, L (⟨fun j => if j = i then 1 else 0, 0⟩) = gradQ H (γ t) i) ∧ - (∀ i, L (⟨0, fun j => if j = i then 1 else 0⟩) = gradP H (γ t) i)) : - HasDerivAt (fun s => H (γ s)) 0 t := by - rcases h_diff with ⟨L, hL, hLq, hLp⟩ - have h_chain := HasFDerivAt.comp_hasDerivAt t hL (h_traj t) - have h_zero : L (HamiltonianVectorField H (γ t)) = 0 := by - rw [clm_on_phaseSpace L (HamiltonianVectorField H (γ t))] - simp only [HamiltonianVectorField] - rw [← Finset.sum_add_distrib] - apply Finset.sum_congr rfl - intro i _ - rw [hLq i, hLp i] - ring - simpa [h_zero] using h_chain - --- ============================================================================ --- Poisson Brackets --- ============================================================================ - -/-- The Poisson bracket of two scalar functions f and g on phase space. -In coordinates: {f, g} = Σ_i (∂f/∂q_i · ∂g/∂p_i - ∂f/∂p_i · ∂g/∂q_i). -/ -noncomputable def poissonBracket {n : ℕ} (f g : PhaseSpace n → ℝ) (x : PhaseSpace n) : ℝ := - ∑ i : Fin n, (gradQ f x i * gradP g x i - gradP f x i * gradQ g x i) - -/-- The Poisson bracket is antisymmetric: {f, g} = -{g, f}. -/ -theorem poissonBracket_antisym {n : ℕ} (f g : PhaseSpace n → ℝ) (x : PhaseSpace n) : - poissonBracket f g x = -poissonBracket g f x := by - unfold poissonBracket - rw [← Finset.sum_neg_distrib] - apply Finset.sum_congr rfl - intro i _ - ring - -/-- The Poisson bracket of any function with itself vanishes: {H, H} = 0. -This is an immediate consequence of antisymmetry. -/ -theorem poissonBracket_self_zero {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) : - poissonBracket H H x = 0 := by - have h : poissonBracket H H x = -poissonBracket H H x := by - rw [poissonBracket_antisym] - linarith - -/-- A quantity f is conserved along the Hamiltonian flow of H if and only if -its Poisson bracket with H vanishes: {f, H} = 0. - -This follows from the chain rule: df/dt = {f, H} along trajectories. -The proof requires the multivariable chain rule (as in energy_conservation), -so we state this as a definition of what it means to be "conserved". -/ -def IsConserved {n : ℕ} (f H : PhaseSpace n → ℝ) : Prop := - ∀ x, poissonBracket f H x = 0 - --- ============================================================================ --- Flow of a Vector Field --- ============================================================================ - -/-- A flow of a vector field V on phase space is a family of maps ϕ^t such that -for every initial point x, the curve t ↦ ϕ^t(x) is a solution of the ODE -with initial condition x. -/ -def IsFlow {n : ℕ} (V : PhaseSpace n → PhaseSpace n) (ϕ : ℝ → PhaseSpace n → PhaseSpace n) : Prop := - (∀ x, ϕ 0 x = x) ∧ (∀ x t, HasDerivAt (fun s => ϕ s x) (V (ϕ t x)) t) - -/-- The identity property of the flow: ϕ^0 = id. -/ -theorem flow_identity {n : ℕ} (V : PhaseSpace n → PhaseSpace n) (ϕ : ℝ → PhaseSpace n → PhaseSpace n) - (h_flow : IsFlow V ϕ) : - ∀ x, ϕ 0 x = x := by - intro x - exact (h_flow.1 x) - -/-- The semigroup property of the flow: ϕ^{s+t} = ϕ^s ∘ ϕ^t. -This is one of the fundamental properties of flows of ODEs. - -The proof uses the uniqueness theorem picard_lindelof_uniqueness: -both sides are solutions of the same ODE with the same initial condition -ϕ^t(x) at time 0, hence they coincide. -/ -theorem flow_semigroup {n : ℕ} (V : PhaseSpace n → PhaseSpace n) (ϕ : ℝ → PhaseSpace n → PhaseSpace n) - (h_flow : IsFlow V ϕ) - (h_lipschitz : ∃ K : ℝ, K ≥ 0 ∧ ∀ x y, ‖V x - V y‖ ≤ K * ‖x - y‖) : - ∀ x s t, ϕ (s + t) x = ϕ s (ϕ t x) := by - intro x s t - rcases h_lipschitz with ⟨K, hK, hV_lip⟩ - -- Both sides satisfy the ODE with initial condition ϕ^t(x) at time 0. - let γ₁ := fun u : ℝ => ϕ (t + u) x - let γ₂ := fun u : ℝ => ϕ u (ϕ t x) - have hγ₁_init : γ₁ 0 = ϕ t x := by simp [γ₁] - have hγ₂_init : γ₂ 0 = ϕ t x := by simp [γ₂, h_flow.1] - have hγ₁_sol : IsSolutionND (fun y => V y) (ϕ t x) γ₁ := by - constructor - · exact hγ₁_init - · intro u - have h := h_flow.2 x (t + u) - have h' : HasDerivAt (fun v => ϕ v x) (V (ϕ (t + u) x)) (t + u) := by - simpa using h - have h_comp : HasDerivAt (fun v => ϕ (t + v) x) (V (ϕ (t + u) x)) u := by - have h_shift : HasDerivAt (fun v => t + v) 1 u := hasDerivAt_id' u - exact HasDerivAt.comp u h' h_shift - simpa [γ₁] using h_comp - have hγ₂_sol : IsSolutionND (fun y => V y) (ϕ t x) γ₂ := by - constructor - · exact hγ₂_init - · intro u - have h := h_flow.2 (ϕ t x) u - simpa [γ₂] using h - have heq := picard_lindelof_uniqueness (fun y => V y) (ϕ t x) K hK hV_lip γ₁ γ₂ hγ₁_sol hγ₂_sol - have h_at_s : γ₁ s = γ₂ s := by - rw [heq] - simp [γ₁, γ₂] at h_at_s - exact h_at_s - --- ============================================================================ --- Liouville's Theorem (Symplecticity of Hamiltonian Flow) --- ============================================================================ - -/-- A map F : PhaseSpace n → PhaseSpace n is symplectic if it preserves -the canonical symplectic form at the linearized (fderiv) level: -ω(DF(u), DF(v)) = ω(u, v) for all tangent vectors u, v. -/ -def IsSymplecticMap {n : ℕ} (F : PhaseSpace n → PhaseSpace n) : Prop := - ∀ x u v, ω (fderiv ℝ F x u) (fderiv ℝ F x v) = ω u v - --- ============================================================================ --- Hamiltonian Hessian and Linearized Vector Field --- ============================================================================ - -/-- Second partial derivative ∂²H/∂q_i∂q_j, defined as the derivative of -(gradQ H · e_i) with respect to q_j. -/ -noncomputable def hessianQQ {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) (i j : Fin n) : ℝ := - if h : ∃ d, HasDerivAt (fun t => gradQ H (replaceQ x j t) i) d (x.1 j) then h.choose else 0 - -/-- Second partial derivative ∂²H/∂p_i∂p_j. -/ -noncomputable def hessianPP {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) (i j : Fin n) : ℝ := - if h : ∃ d, HasDerivAt (fun t => gradP H (replaceP x j t) i) d (x.2 j) then h.choose else 0 - -/-- Second partial derivative ∂²H/∂q_i∂p_j. -/ -noncomputable def hessianQP {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) (i j : Fin n) : ℝ := - if h : ∃ d, HasDerivAt (fun t => gradQ H (replaceP x j t) i) d (x.2 j) then h.choose else 0 - -/-- Second partial derivative ∂²H/∂p_i∂q_j. -/ -noncomputable def hessianPQ {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) (i j : Fin n) : ℝ := - if h : ∃ d, HasDerivAt (fun t => gradP H (replaceQ x j t) i) d (x.1 j) then h.choose else 0 - -/-- The derivative of the Hamiltonian vector field X_H at x applied to v. -In coordinates: - (DX_H · v)_q_i = Σ_j (hessianPQ_ij · v_q_j + hessianPP_ij · v_p_j) - (DX_H · v)_p_i = Σ_j (−hessianQQ_ij · v_q_j − hessianQP_ij · v_p_j) -/ -def DX_H_action {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) (v : PhaseSpace n) : PhaseSpace n := - ⟨fun i => ∑ j : Fin n, (hessianPQ H x i j * v.1 j + hessianPP H x i j * v.2 j), - fun i => ∑ j : Fin n, (-hessianQQ H x i j * v.1 j - hessianQP H x i j * v.2 j)⟩ - -/-- The Hamiltonian vector field derivative is skew-symmetric with respect to ω: -ω(DX_H · u, v) + ω(u, DX_H · v) = 0. - -This is the infinitesimal algebraic identity underlying symplectic preservation. -The proof expands both terms and uses Clairaut's theorem (symmetry of mixed -partials) to show all terms cancel in pairs. -/ -lemma hamiltonian_skew_symmetric {n : ℕ} (H : PhaseSpace n → ℝ) (x : PhaseSpace n) - (u v : PhaseSpace n) - (h_clairaut : ∀ i j, hessianQQ H x i j = hessianQQ H x j i - ∧ hessianPP H x i j = hessianPP H x j i - ∧ hessianQP H x i j = hessianPQ H x j i) : - ω (DX_H_action H x u) v + ω u (DX_H_action H x v) = 0 := by - -- Expand both terms fully using definitions - simp only [ω, DX_H_action] - -- The QQ terms cancel by antisymmetry of (u_q_i·v_q_j − u_q_j·v_q_i) + symmetry of hessianQQ - have h_QQ : ∑ i : Fin n, ∑ j : Fin n, - (-hessianQQ H x i j * u.1 j * v.1 i + hessianQQ H x i j * u.1 i * v.1 j) = 0 := by - have h_symm : ∀ i j, hessianQQ H x i j = hessianQQ H x j i := fun i j => (h_clairaut i j).1 - have h_total : ∑ i, ∑ j, hessianQQ H x i j * (u.1 i * v.1 j - u.1 j * v.1 i) = 0 := by - have h1 : ∑ i, ∑ j, hessianQQ H x i j * (u.1 i * v.1 j - u.1 j * v.1 i) - = ∑ j, ∑ i, hessianQQ H x j i * (u.1 j * v.1 i - u.1 i * v.1 j) := by - rw [Finset.sum_comm] - apply Finset.sum_congr rfl - intro i _ - apply Finset.sum_congr rfl - intro j _ - ring - have h2 : ∑ j, ∑ i, hessianQQ H x j i * (u.1 j * v.1 i - u.1 i * v.1 j) - = -∑ i, ∑ j, hessianQQ H x i j * (u.1 i * v.1 j - u.1 j * v.1 i) := by - simp_rw [h_symm] - have h_neg : ∀ i j, u.1 j * v.1 i - u.1 i * v.1 j = -(u.1 i * v.1 j - u.1 j * v.1 i) := by - intro i j; ring - simp_rw [h_neg, mul_neg, Finset.sum_neg_distrib] - linarith [h1, h2] - have h_eq : ∀ i j, -hessianQQ H x i j * u.1 j * v.1 i + hessianQQ H x i j * u.1 i * v.1 j - = hessianQQ H x i j * (u.1 i * v.1 j - u.1 j * v.1 i) := by - intro i j; ring - simp_rw [h_eq] - exact h_total - -- The PP terms: same argument - have h_PP : ∑ i : Fin n, ∑ j : Fin n, - (-hessianPP H x i j * u.2 j * v.2 i + hessianPP H x i j * u.2 i * v.2 j) = 0 := by - have h_symm : ∀ i j, hessianPP H x i j = hessianPP H x j i := fun i j => (h_clairaut i j).2.1 - have h_total : ∑ i, ∑ j, hessianPP H x i j * (u.2 i * v.2 j - u.2 j * v.2 i) = 0 := by - have h1 : ∑ i, ∑ j, hessianPP H x i j * (u.2 i * v.2 j - u.2 j * v.2 i) - = ∑ j, ∑ i, hessianPP H x j i * (u.2 j * v.2 i - u.2 i * v.2 j) := by - rw [Finset.sum_comm] - apply Finset.sum_congr rfl - intro i _ - apply Finset.sum_congr rfl - intro j _ - ring - have h2 : ∑ j, ∑ i, hessianPP H x j i * (u.2 j * v.2 i - u.2 i * v.2 j) - = -∑ i, ∑ j, hessianPP H x i j * (u.2 i * v.2 j - u.2 j * v.2 i) := by - simp_rw [h_symm] - have h_neg : ∀ i j, u.2 j * v.2 i - u.2 i * v.2 j = -(u.2 i * v.2 j - u.2 j * v.2 i) := by - intro i j; ring - simp_rw [h_neg, mul_neg, Finset.sum_neg_distrib] - linarith [h1, h2] - have h_eq : ∀ i j, -hessianPP H x i j * u.2 j * v.2 i + hessianPP H x i j * u.2 i * v.2 j - = hessianPP H x i j * (u.2 i * v.2 j - u.2 j * v.2 i) := by - intro i j; ring - simp_rw [h_eq] - exact h_total - -- The cross terms (QP + PQ): cancel using hessianQP_ij = hessianPQ_ji - have h_cross : ∑ i : Fin n, ∑ j : Fin n, - (-hessianQP H x i j * u.2 j * v.1 i + hessianQP H x i j * u.1 i * v.2 j - - hessianPQ H x i j * u.1 j * v.2 i + hessianPQ H x i j * u.2 i * v.1 j) = 0 := by - have h_cross_symm : ∀ i j, hessianQP H x i j = hessianPQ H x j i := fun i j => (h_clairaut j i).2.2 - have h_eq : ∀ i j, -hessianQP H x i j * u.2 j * v.1 i + hessianQP H x i j * u.1 i * v.2 j - - hessianPQ H x i j * u.1 j * v.2 i + hessianPQ H x i j * u.2 i * v.1 j - = hessianPQ H x j i * (-u.2 j * v.1 i + u.1 i * v.2 j) - + hessianPQ H x i j * (-u.1 j * v.2 i + u.2 i * v.1 j) := by - intro i j - rw [h_cross_symm i j] - ring - simp_rw [h_eq] - have h_total : ∑ i, ∑ j, hessianPQ H x j i * (-u.2 j * v.1 i + u.1 i * v.2 j) - + ∑ i, ∑ j, hessianPQ H x i j * (-u.1 j * v.2 i + u.2 i * v.1 j) = 0 := by - have h_swap : ∑ i, ∑ j, hessianPQ H x j i * (-u.2 j * v.1 i + u.1 i * v.2 j) - = ∑ i, ∑ j, hessianPQ H x i j * (u.2 i * v.1 j - u.1 j * v.2 i) := by - rw [Finset.sum_comm] - apply Finset.sum_congr rfl - intro i _ - apply Finset.sum_congr rfl - intro j _ - ring - rw [h_swap] - have h_cancel : ∀ i j, hessianPQ H x i j * (u.2 i * v.1 j - u.1 j * v.2 i) - + hessianPQ H x i j * (-u.1 j * v.2 i + u.2 i * v.1 j) = 0 := by - intro i j - ring - simp_rw [h_cancel] - simp - linarith [h_total] - -- Rearrange each (i,j) term into the four groups - have h_rearrange : ∀ i j, - (-hessianQQ H x i j * u.1 j * v.1 i - hessianQP H x i j * u.2 j * v.1 i - - hessianPQ H x i j * u.1 j * v.2 i - hessianPP H x i j * u.2 j * v.2 i - + hessianPQ H x i j * u.2 i * v.1 j + hessianPP H x i j * u.2 i * v.2 j - + hessianQQ H x i j * u.1 i * v.1 j + hessianQP H x i j * u.1 i * v.2 j) - = hessianQQ H x i j * (u.1 i * v.1 j - u.1 j * v.1 i) - + hessianPP H x i j * (u.2 i * v.2 j - u.2 j * v.2 i) - + hessianPQ H x i j * (u.2 i * v.1 j - u.1 j * v.2 i) - + hessianQP H x i j * (u.1 i * v.2 j - u.2 j * v.1 i) := by - intro i j; ring - simp_rw [h_rearrange] - -- Split into four double sums - rw [← Finset.sum_add_distrib, ← Finset.sum_add_distrib] - -- Match with the three cancellation lemmas - have h_sum3_4 : (∑ i : Fin n, ∑ j : Fin n, hessianPQ H x i j * (u.2 i * v.1 j - u.1 j * v.2 i)) - + (∑ i : Fin n, ∑ j : Fin n, hessianQP H x i j * (u.1 i * v.2 j - u.2 j * v.1 i)) - = ∑ i : Fin n, ∑ j : Fin n, - (-hessianQP H x i j * u.2 j * v.1 i + hessianQP H x i j * u.1 i * v.2 j - - hessianPQ H x i j * u.1 j * v.2 i + hessianPQ H x i j * u.2 i * v.1 j) := by - rw [← Finset.sum_add_distrib] - apply Finset.sum_congr rfl - intro i _ - rw [← Finset.sum_add_distrib] - apply Finset.sum_congr rfl - intro j _ - ring - linarith [h_QQ, h_PP, h_cross, h_sum3_4] - -/-- Derivative of the symplectic form composed with two curves: -d/dt ω(a(t), b(t)) = ω(a'(t), b(t)) + ω(a(t), b'(t)). -/ -lemma deriv_of_bilinear_ω {n : ℕ} (a b : ℝ → PhaseSpace n) (a' b' : PhaseSpace n) (t : ℝ) - (ha : HasDerivAt a a' t) (hb : HasDerivAt b b' t) : - HasDerivAt (fun s => ω (a s) (b s)) (ω a' (b t) + ω (a t) b') t := by - -- Expand ω as a sum of products - have h_expand : ∀ s, ω (a s) (b s) = ∑ i : Fin n, ((a s).2 i * (b s).1 i - (a s).1 i * (b s).2 i) := by - intro s; rfl - rw [show (fun s => ω (a s) (b s)) = fun s => ∑ i, ((a s).2 i * (b s).1 i - (a s).1 i * (b s).2 i) - by funext s; apply h_expand] - -- Extract component-wise derivatives using projection CLMs - have ha1 i : HasDerivAt (fun s => (a s).1 i) (a'.1 i) t := by - let eval_i : (Fin n → ℝ) →L[ℝ] ℝ := - ⟨fun f => f i, by intros; rfl, by intros; rfl, by continuity⟩ - have heval : HasFDerivAt eval_i eval_i ((a t).1) := eval_i.hasFDerivAt - have ha_fst : HasDerivAt (fun s => (a s).1) a'.1 t := HasDerivAt.fst ha - exact HasFDerivAt.comp_hasDerivAt t heval ha_fst - have ha2 i : HasDerivAt (fun s => (a s).2 i) (a'.2 i) t := by - let eval_i : (Fin n → ℝ) →L[ℝ] ℝ := - ⟨fun f => f i, by intros; rfl, by intros; rfl, by continuity⟩ - have heval : HasFDerivAt eval_i eval_i ((a t).2) := eval_i.hasFDerivAt - have ha_snd : HasDerivAt (fun s => (a s).2) a'.2 t := HasDerivAt.snd ha - exact HasFDerivAt.comp_hasDerivAt t heval ha_snd - have hb1 i : HasDerivAt (fun s => (b s).1 i) (b'.1 i) t := by - let eval_i : (Fin n → ℝ) →L[ℝ] ℝ := - ⟨fun f => f i, by intros; rfl, by intros; rfl, by continuity⟩ - have heval : HasFDerivAt eval_i eval_i ((b t).1) := eval_i.hasFDerivAt - have hb_fst : HasDerivAt (fun s => (b s).1) b'.1 t := HasDerivAt.fst hb - exact HasFDerivAt.comp_hasDerivAt t heval hb_fst - have hb2 i : HasDerivAt (fun s => (b s).2 i) (b'.2 i) t := by - let eval_i : (Fin n → ℝ) →L[ℝ] ℝ := - ⟨fun f => f i, by intros; rfl, by intros; rfl, by continuity⟩ - have heval : HasFDerivAt eval_i eval_i ((b t).2) := eval_i.hasFDerivAt - have hb_snd : HasDerivAt (fun s => (b s).2) b'.2 t := HasDerivAt.snd hb - exact HasFDerivAt.comp_hasDerivAt t heval hb_snd - -- Differentiate the sum term by term using product rule - have h_deriv : HasDerivAt (fun s => ∑ i, ((a s).2 i * (b s).1 i - (a s).1 i * (b s).2 i)) - (∑ i, (a'.2 i * (b t).1 i + (a t).2 i * b'.1 i - a'.1 i * (b t).2 i - (a t).1 i * b'.2 i)) t := by - apply HasDerivAt.sum - intro i _ - apply HasDerivAt.sub - · apply HasDerivAt.mul (ha2 i) (hb1 i) - · apply HasDerivAt.mul (ha1 i) (hb2 i) - convert h_deriv using 1 - -- Show the derivative equals ω(a', b) + ω(a, b') - simp [ω] - rw [← Finset.sum_add_distrib] - apply Finset.sum_congr rfl - intro i _ - ring - --- ============================================================================ --- Liouville's Theorem (Symplecticity of Hamiltonian Flow) --- ============================================================================ - -/-- A map F : PhaseSpace n → PhaseSpace n is symplectic if it preserves -the canonical symplectic form at the linearized (fderiv) level: -ω(DF(u), DF(v)) = ω(u, v) for all tangent vectors u, v. -/ -def IsSymplecticMap {n : ℕ} (F : PhaseSpace n → PhaseSpace n) : Prop := - ∀ x u v, ω (fderiv ℝ F x u) (fderiv ℝ F x v) = ω u v - -/-- Liouville's theorem: the flow of a Hamiltonian vector field preserves -the canonical symplectic form. - -Proof ingredients: -1. Clairaut's theorem (symmetry of mixed partials) — hypothesis `h_clairaut` -2. The variational equation — hypothesis `h_variational` -3. The algebraic identity ω(DX_H·u, v) + ω(u, DX_H·v) = 0 — proven in `hamiltonian_skew_symmetric` - -The proof shows d/dt ω(Dϕ^t·u, Dϕ^t·v) = 0, so the form is constant in time, -equal to its value at t=0 which is ω(u,v). -/ -theorem liouville_theorem {n : ℕ} (H : PhaseSpace n → ℝ) (ϕ : ℝ → PhaseSpace n → PhaseSpace n) - (h_flow : IsFlow (HamiltonianVectorField H) ϕ) - (h_lipschitz : ∃ K : ℝ, K ≥ 0 ∧ ∀ x y, ‖HamiltonianVectorField H x - HamiltonianVectorField H y‖ ≤ K * ‖x - y‖) - (h_clairaut : ∀ x i j, hessianQQ H x i j = hessianQQ H x j i - ∧ hessianPP H x i j = hessianPP H x j i - ∧ hessianQP H x i j = hessianPQ H x j i) - (h_variational : ∀ t x v, HasDerivAt (fun s => fderiv ℝ (ϕ s) x v) - (DX_H_action H (ϕ t x) (fderiv ℝ (ϕ t) x v)) t) : - ∀ t, IsSymplecticMap (ϕ t) := by - intro t x u v - -- h(s) = ω(Dϕ^s(x)·u, Dϕ^s(x)·v) - let h := fun s : ℝ => ω (fderiv ℝ (ϕ s) x u) (fderiv ℝ (ϕ s) x v) - -- Step 1: h(0) = ω(u, v) since ϕ^0 = id - have h0 : h 0 = ω u v := by - have h_id : fderiv ℝ (ϕ 0) x = ContinuousLinearMap.id ℝ (PhaseSpace n) := by - have hϕ0 : ϕ 0 = id := by funext y; exact h_flow.1 y - rw [hϕ0] - exact fderiv_id - simp [h, h_id] - -- Step 2: h'(s) = 0 for all s (variational equation + skew-symmetry) - have h_deriv0 : ∀ s, HasDerivAt h 0 s := by - intro s - let a := fun (σ : ℝ) => fderiv ℝ (ϕ σ) x u - let b := fun (σ : ℝ) => fderiv ℝ (ϕ σ) x v - let a' := DX_H_action H (ϕ s x) (fderiv ℝ (ϕ s) x u) - let b' := DX_H_action H (ϕ s x) (fderiv ℝ (ϕ s) x v) - have ha : HasDerivAt a a' s := h_variational s x u - have hb : HasDerivAt b b' s := h_variational s x v - have h_ω : HasDerivAt (fun σ => ω (a σ) (b σ)) (ω a' (b s) + ω (a s) b') s := - deriv_of_bilinear_ω a b a' b' s ha hb - have h_key : ω a' (b s) + ω (a s) b' = 0 := by - exact hamiltonian_skew_symmetric H (ϕ s x) (a s) (b s) (h_clairaut (ϕ s x)) - simp [h] at h_ω ⊢ - rw [h_key] at h_ω - simpa using h_ω - -- Step 3: h is constant (zero derivative everywhere) via mean value theorem - have h_const : ∀ s, h s = h 0 := by - intro s - by_cases h_eq : s = 0 - · rw [h_eq] - · let g := fun τ => h τ - h 0 - have hg0 : g 0 = 0 := by simp [g] - have hg_deriv : ∀ τ, HasDerivAt g 0 τ := by - intro τ - have hh := h_deriv0 τ - simpa [g] using HasDerivAt.sub hh (hasDerivAt_const τ (h 0)) - have hdiff : Differentiable ℝ g := fun τ => (hg_deriv τ).differentiableAt - have hcont := hdiff.continuous - by_cases hs : s > 0 - · have hcont_on : ContinuousOn g (Set.Icc 0 s) := hcont.continuousOn - have hder : ∀ τ ∈ Set.Ioo 0 s, HasDerivAt g 0 τ := fun τ _ => hg_deriv τ - have hmvt : ∃ c ∈ Set.Ioo 0 s, g s - g 0 = 0 * (s - 0) := by - apply exists_hasDerivAt_eq_slope g (fun _ => 0) 0 s hs hcont_on - intro τ hτ; exact hder τ hτ - rcases hmvt with ⟨c, _, heq⟩ - simp at heq - linarith [heq, hg0] - · have hs' : s < 0 := by - have hne : s ≠ 0 := h_eq - have hle : s ≤ 0 := le_of_not_gt hs - exact lt_of_le_of_ne hle (Ne.symm hne) - have hcont_on : ContinuousOn g (Set.Icc s 0) := hcont.continuousOn - have hder : ∀ τ ∈ Set.Ioo s 0, HasDerivAt g 0 τ := fun τ _ => hg_deriv τ - have hmvt : ∃ c ∈ Set.Ioo s 0, g 0 - g s = 0 * (0 - s) := by - apply exists_hasDerivAt_eq_slope g (fun _ => 0) s 0 hs' hcont_on - intro τ hτ; exact hder τ hτ - rcases hmvt with ⟨c, _, heq⟩ - simp at heq - linarith [heq, hg0] - -- Step 4: Conclude h(t) = h(0) = ω(u, v) - rw [h_const t, h0] - rfl - -end Semantics.AnalysisFoundations