IndisputableMonolith.Relativity.Calculus.Derivatives
IndisputableMonolith/Relativity/Calculus/Derivatives.lean · 839 lines · 51 declarations
show as:
view math explainer →
1import Mathlib
2import IndisputableMonolith.Relativity.Geometry.Tensor
3
4namespace IndisputableMonolith
5namespace Relativity
6namespace Calculus
7
8open scoped Topology
9open Filter Real
10
11/-- Standard basis vector `e_μ`. -/
12def basisVec (μ : Fin 4) : Fin 4 → ℝ := fun ν => if ν = μ then 1 else 0
13
14@[simp] lemma basisVec_self (μ : Fin 4) : basisVec μ μ = 1 := by simp [basisVec]
15
16@[simp] lemma basisVec_ne {μ ν : Fin 4} (h : ν ≠ μ) : basisVec μ ν = 0 := by
17 simp [basisVec, h]
18
19/-- Coordinate ray `x + t e_μ`. -/
20def coordRay (x : Fin 4 → ℝ) (μ : Fin 4) (t : ℝ) : Fin 4 → ℝ :=
21 fun ν => x ν + t * basisVec μ ν
22
23@[simp] lemma coordRay_apply (x : Fin 4 → ℝ) (μ : Fin 4) (t : ℝ) (ν : Fin 4) :
24 coordRay x μ t ν = x ν + t * basisVec μ ν := rfl
25
26@[simp] lemma coordRay_zero (x : Fin 4 → ℝ) (μ : Fin 4) : coordRay x μ 0 = x := by
27 funext ν; simp [coordRay]
28
29@[simp] lemma coordRay_coordRay (x : Fin 4 → ℝ) (μ : Fin 4) (s t : ℝ) :
30 coordRay (coordRay x μ s) μ t = coordRay x μ (s + t) := by
31 funext ν; simp [coordRay]; ring
32
33/-- Directional derivative `∂_μ f(x)` via real derivative along the coordinate ray. -/
34noncomputable def partialDeriv_v2 (f : (Fin 4 → ℝ) → ℝ) (μ : Fin 4)
35 (x : Fin 4 → ℝ) : ℝ :=
36 deriv (fun t => f (coordRay x μ t)) 0
37
38/-- The derivative of a constant function is zero. -/
39lemma partialDeriv_v2_const {f : (Fin 4 → ℝ) → ℝ} {c : ℝ} (h : ∀ y, f y = c) (μ : Fin 4) (x : Fin 4 → ℝ) :
40 partialDeriv_v2 f μ x = 0 := by
41 unfold partialDeriv_v2
42 have h_const : (fun t => f (coordRay x μ t)) = (fun _ => c) := by
43 funext t
44 rw [h]
45 rw [h_const]
46 exact deriv_const (0 : ℝ) c
47
48/-- Second derivative `∂_μ∂_ν f(x)` as iterated directional derivatives. -/
49noncomputable def secondDeriv (f : (Fin 4 → ℝ) → ℝ) (μ ν : Fin 4)
50 (x : Fin 4 → ℝ) : ℝ :=
51 deriv (fun s => partialDeriv_v2 f μ (coordRay x ν s)) 0
52
53/-- Laplacian `∇² = Σ_{i=1}^3 ∂²/∂xᵢ²`. -/
54noncomputable def laplacian (f : (Fin 4 → ℝ) → ℝ) (x : Fin 4 → ℝ) : ℝ :=
55 secondDeriv f ⟨1, by decide⟩ ⟨1, by decide⟩ x +
56 secondDeriv f ⟨2, by decide⟩ ⟨2, by decide⟩ x +
57 secondDeriv f ⟨3, by decide⟩ ⟨3, by decide⟩ x
58
59/-- Linearity of the directional derivative. -/
60lemma deriv_add_lin (f g : (Fin 4 → ℝ) → ℝ) (μ : Fin 4)
61 (x : Fin 4 → ℝ) (hf : DifferentiableAt ℝ (fun t => f (coordRay x μ t)) 0)
62 (hg : DifferentiableAt ℝ (fun t => g (coordRay x μ t)) 0) :
63 partialDeriv_v2 (fun y => f y + g y) μ x =
64 partialDeriv_v2 f μ x + partialDeriv_v2 g μ x := by
65 unfold partialDeriv_v2
66 exact deriv_add hf hg
67
68/-- Linearity of directional derivative (scalar multiplication). -/
69lemma partialDeriv_v2_smul (f : (Fin 4 → ℝ) → ℝ) (c : ℝ) (μ : Fin 4)
70 (x : Fin 4 → ℝ) (hf : DifferentiableAt ℝ (fun t => f (coordRay x μ t)) 0) :
71 partialDeriv_v2 (fun y => c * f y) μ x = c * partialDeriv_v2 f μ x := by
72 unfold partialDeriv_v2
73 exact deriv_const_mul c hf
74
75/-- Localized version of second derivative linearity (scalar multiplication).
76 This only requires differentiability in a neighborhood of the point x. -/
77lemma secondDeriv_smul_local (f : (Fin 4 → ℝ) → ℝ) (c : ℝ) (μ ν : Fin 4)
78 (x : Fin 4 → ℝ)
79 (h1 : ∀ᶠ s in 𝓝 0, DifferentiableAt ℝ (fun t => f (coordRay (coordRay x ν s) μ t)) 0)
80 (h2 : DifferentiableAt ℝ (fun s => partialDeriv_v2 f μ (coordRay x ν s)) 0) :
81 secondDeriv (fun y => c * f y) μ ν x = c * secondDeriv f μ ν x := by
82 unfold secondDeriv
83 have h_ev : ∀ᶠ s in 𝓝 0, partialDeriv_v2 (fun z => c * f z) μ (coordRay x ν s) =
84 c * partialDeriv_v2 f μ (coordRay x ν s) := by
85 apply h1.mono
86 intro s hs
87 exact partialDeriv_v2_smul f c μ (coordRay x ν s) hs
88 rw [Filter.EventuallyEq.deriv_eq h_ev]
89 exact deriv_const_mul c h2
90
91/-- Second derivative linearity (scalar multiplication). -/
92lemma secondDeriv_smul (f : (Fin 4 → ℝ) → ℝ) (c : ℝ) (μ ν : Fin 4)
93 (x : Fin 4 → ℝ)
94 (h1 : ∀ y, DifferentiableAt ℝ (fun t => f (coordRay y μ t)) 0)
95 (h2 : DifferentiableAt ℝ (fun s => partialDeriv_v2 f μ (coordRay x ν s)) 0) :
96 secondDeriv (fun y => c * f y) μ ν x = c * secondDeriv f μ ν x := by
97 unfold secondDeriv
98 have h_partial : ∀ y, partialDeriv_v2 (fun z => c * f z) μ y = c * partialDeriv_v2 f μ y := by
99 intro y
100 exact partialDeriv_v2_smul f c μ y (h1 y)
101 simp only [h_partial]
102 exact deriv_const_mul c h2
103
104/-- Laplacian linearity (scalar multiplication). -/
105lemma laplacian_smul (f : (Fin 4 → ℝ) → ℝ) (c : ℝ) (x : Fin 4 → ℝ)
106 (h1 : ∀ μ y, DifferentiableAt ℝ (fun t => f (coordRay y μ t)) 0)
107 (h2 : ∀ μ ν, DifferentiableAt ℝ (fun s => partialDeriv_v2 f μ (coordRay x ν s)) 0) :
108 laplacian (fun y => c * f y) x = c * laplacian f x := by
109 unfold laplacian
110 simp only [secondDeriv_smul f c _ _ x (h1 _) (h2 _ _)]
111 ring
112
113/-- Product rule for directional derivative. -/
114lemma partialDeriv_v2_mul (f g : (Fin 4 → ℝ) → ℝ) (μ : Fin 4)
115 (x : Fin 4 → ℝ) (hf : DifferentiableAt ℝ (fun t => f (coordRay x μ t)) 0)
116 (hg : DifferentiableAt ℝ (fun t => g (coordRay x μ t)) 0) :
117 partialDeriv_v2 (fun y => f y * g y) μ x =
118 f x * partialDeriv_v2 g μ x + g x * partialDeriv_v2 f μ x := by
119 unfold partialDeriv_v2
120 have h_mul : deriv (fun ε => f (coordRay x μ ε) * g (coordRay x μ ε)) 0 =
121 deriv (fun ε => f (coordRay x μ ε)) 0 * g (coordRay x μ 0) +
122 f (coordRay x μ 0) * deriv (fun ε => g (coordRay x μ ε)) 0 :=
123 deriv_mul hf hg
124 rw [h_mul]
125 simp only [coordRay_zero]
126 ring
127
128/-- Spatial norm squared `x₁² + x₂² + x₃²`. -/
129def spatialNormSq (x : Fin 4 → ℝ) : ℝ := x 1 ^ 2 + x 2 ^ 2 + x 3 ^ 2
130
131theorem spatialNormSq_nonneg (x : Fin 4 → ℝ) : 0 ≤ spatialNormSq x := by
132 unfold spatialNormSq
133 positivity
134
135theorem spatialNormSq_eq_zero_iff (x : Fin 4 → ℝ) : spatialNormSq x = 0 ↔ x 1 = 0 ∧ x 2 = 0 ∧ x 3 = 0 := by
136 unfold spatialNormSq
137 constructor
138 · intro h
139 have h1 := sq_nonneg (x 1)
140 have h2 := sq_nonneg (x 2)
141 have h3 := sq_nonneg (x 3)
142 have h1_zero : x 1 ^ 2 = 0 := by linarith
143 have h2_zero : x 2 ^ 2 = 0 := by linarith
144 have h3_zero : x 3 ^ 2 = 0 := by linarith
145 simp only [sq_eq_zero_iff] at h1_zero h2_zero h3_zero
146 exact ⟨h1_zero, h2_zero, h3_zero⟩
147 · intro h
148 simp [h]
149
150/-- Spatial radius `r = √(x₁² + x₂² + x₃²)`. -/
151noncomputable def spatialRadius (x : Fin 4 → ℝ) : ℝ := Real.sqrt (spatialNormSq x)
152
153theorem spatialRadius_pos_iff (x : Fin 4 → ℝ) : 0 < spatialRadius x ↔ 0 < spatialNormSq x := by
154 unfold spatialRadius
155 rw [Real.sqrt_pos]
156
157theorem spatialRadius_ne_zero_iff (x : Fin 4 → ℝ) : spatialRadius x ≠ 0 ↔ spatialNormSq x ≠ 0 := by
158 unfold spatialRadius
159 rw [Real.sqrt_ne_zero (spatialNormSq_nonneg x)]
160
161/-- Nonzero spatial radius is automatically positive. -/
162theorem spatialRadius_pos_of_ne_zero (x : Fin 4 → ℝ) (hx : spatialRadius x ≠ 0) :
163 0 < spatialRadius x := by
164 have h_ne : spatialNormSq x ≠ 0 := (spatialRadius_ne_zero_iff x).mp hx
165 exact (spatialRadius_pos_iff x).2 (lt_of_le_of_ne (spatialNormSq_nonneg x) h_ne.symm)
166
167/-- Temporal coordinate ray doesn't change spatial components. -/
168lemma coordRay_temporal_spatial (x : Fin 4 → ℝ) (s : ℝ) (i : Fin 4) (hi : i ≠ 0) :
169 (coordRay x 0 s) i = x i := by
170 simp [coordRay, basisVec, hi]
171
172/-- spatialNormSq is invariant under temporal coordinate ray. -/
173lemma spatialNormSq_coordRay_temporal (x : Fin 4 → ℝ) (s : ℝ) :
174 spatialNormSq (coordRay x 0 s) = spatialNormSq x := by
175 unfold spatialNormSq
176 have h1 : (coordRay x 0 s) 1 = x 1 := coordRay_temporal_spatial x s 1 (by decide)
177 have h2 : (coordRay x 0 s) 2 = x 2 := coordRay_temporal_spatial x s 2 (by decide)
178 have h3 : (coordRay x 0 s) 3 = x 3 := coordRay_temporal_spatial x s 3 (by decide)
179 rw [h1, h2, h3]
180
181/-- spatialRadius is invariant under temporal coordinate ray. -/
182lemma spatialRadius_coordRay_temporal (x : Fin 4 → ℝ) (s : ℝ) :
183 spatialRadius (coordRay x 0 s) = spatialRadius x := by
184 unfold spatialRadius
185 rw [spatialNormSq_coordRay_temporal]
186
187/-- For any spatial index `i ∈ {1,2,3}`, `x_i² ≤ spatialNormSq x`. -/
188private lemma sq_le_spatialNormSq_1 (x : Fin 4 → ℝ) :
189 x 1 ^ 2 ≤ spatialNormSq x := by
190 unfold spatialNormSq; nlinarith [sq_nonneg (x 2), sq_nonneg (x 3)]
191
192private lemma sq_le_spatialNormSq_2 (x : Fin 4 → ℝ) :
193 x 2 ^ 2 ≤ spatialNormSq x := by
194 unfold spatialNormSq; nlinarith [sq_nonneg (x 1), sq_nonneg (x 3)]
195
196private lemma sq_le_spatialNormSq_3 (x : Fin 4 → ℝ) :
197 x 3 ^ 2 ≤ spatialNormSq x := by
198 unfold spatialNormSq; nlinarith [sq_nonneg (x 1), sq_nonneg (x 2)]
199
200/-- Helper: closed form of `spatialNormSq (coordRay x j s)` when `j` is a fixed
201 spatial index. The sum changes only at index `j`, where `x j` becomes `x j + s`. -/
202private lemma spatialNormSq_coordRay_spatial_1 (x : Fin 4 → ℝ) (s : ℝ) :
203 spatialNormSq (coordRay x 1 s) = (x 1 + s) ^ 2 + x 2 ^ 2 + x 3 ^ 2 := by
204 unfold spatialNormSq coordRay basisVec
205 rw [if_pos (rfl : (1 : Fin 4) = 1),
206 if_neg (by decide : (2 : Fin 4) ≠ 1),
207 if_neg (by decide : (3 : Fin 4) ≠ 1)]
208 ring
209
210private lemma spatialNormSq_coordRay_spatial_2 (x : Fin 4 → ℝ) (s : ℝ) :
211 spatialNormSq (coordRay x 2 s) = x 1 ^ 2 + (x 2 + s) ^ 2 + x 3 ^ 2 := by
212 unfold spatialNormSq coordRay basisVec
213 rw [if_neg (by decide : (1 : Fin 4) ≠ 2),
214 if_pos (rfl : (2 : Fin 4) = 2),
215 if_neg (by decide : (3 : Fin 4) ≠ 2)]
216 ring
217
218private lemma spatialNormSq_coordRay_spatial_3 (x : Fin 4 → ℝ) (s : ℝ) :
219 spatialNormSq (coordRay x 3 s) = x 1 ^ 2 + x 2 ^ 2 + (x 3 + s) ^ 2 := by
220 unfold spatialNormSq coordRay basisVec
221 rw [if_neg (by decide : (1 : Fin 4) ≠ 3),
222 if_neg (by decide : (2 : Fin 4) ≠ 3),
223 if_pos (rfl : (3 : Fin 4) = 3)]
224 ring
225
226/-- `spatialRadius` stays nonzero under sufficiently small coordinate perturbations.
227
228 Quantitative version: if `r = spatialRadius x ≠ 0` and `|s| < r/2`, then the
229 perturbed point `coordRay x ν s = x + s · e_ν` still has nonzero spatial radius.
230
231 Proof: case-split on `ν ∈ {0,1,2,3}`.
232 - `ν = 0`: temporal direction, `spatialRadius (coordRay x 0 s) = spatialRadius x` (proved).
233 - `ν ∈ {1,2,3}`: only the `ν`-th spatial component changes by `s`, so
234 `spatialNormSq (coordRay x ν s) = ‖x‖² + 2 s · x_ν + s²`. Using `|x_ν| ≤ r`
235 and `|s| < r/2`, the polynomial lower bound `(r - |s|)² ≤ ‖x‖² + 2 s x_ν + s²`
236 gives `spatialNormSq > 0` and hence `spatialRadius ≠ 0`.
237
238 Closes one of the §XXIII.B′ Mathlib calculus axioms. -/
239theorem spatialRadius_coordRay_ne_zero (x : Fin 4 → ℝ) (ν : Fin 4) (s : ℝ)
240 (hx : spatialRadius x ≠ 0) (hs : |s| < spatialRadius x / 2) :
241 spatialRadius (coordRay x ν s) ≠ 0 := by
242 rw [spatialRadius_ne_zero_iff]
243 have hr_pos : 0 < spatialRadius x := spatialRadius_pos_of_ne_zero x hx
244 have hr_sq : spatialRadius x ^ 2 = spatialNormSq x := by
245 unfold spatialRadius; rw [Real.sq_sqrt (spatialNormSq_nonneg x)]
246 have h_s_lo : -(spatialRadius x / 2) < s := (abs_lt.mp hs).1
247 have h_s_hi : s < spatialRadius x / 2 := (abs_lt.mp hs).2
248 have h_x1_le : x 1 ^ 2 ≤ spatialRadius x ^ 2 := hr_sq ▸ sq_le_spatialNormSq_1 x
249 have h_x2_le : x 2 ^ 2 ≤ spatialRadius x ^ 2 := hr_sq ▸ sq_le_spatialNormSq_2 x
250 have h_x3_le : x 3 ^ 2 ≤ spatialRadius x ^ 2 := hr_sq ▸ sq_le_spatialNormSq_3 x
251 -- Case split on ν using its underlying ℕ.
252 obtain ⟨k, hk⟩ := ν
253 match k, hk with
254 | 0, _ =>
255 rw [show ((⟨0, by decide⟩ : Fin 4) = (0 : Fin 4)) from rfl,
256 spatialNormSq_coordRay_temporal]
257 exact (spatialRadius_ne_zero_iff x).mp hx
258 | 1, _ =>
259 rw [show ((⟨1, by decide⟩ : Fin 4) = (1 : Fin 4)) from rfl,
260 spatialNormSq_coordRay_spatial_1]
261 intro h_zero
262 have h_sq1 : (x 1 + s) ^ 2 = 0 := by
263 nlinarith [sq_nonneg (x 1 + s), sq_nonneg (x 2), sq_nonneg (x 3)]
264 have h_sq2 : x 2 ^ 2 = 0 := by
265 nlinarith [sq_nonneg (x 1 + s), sq_nonneg (x 2), sq_nonneg (x 3)]
266 have h_sq3 : x 3 ^ 2 = 0 := by
267 nlinarith [sq_nonneg (x 1 + s), sq_nonneg (x 2), sq_nonneg (x 3)]
268 have h_eq1 : x 1 + s = 0 := sq_eq_zero_iff.mp h_sq1
269 have h_eq2 : x 2 = 0 := sq_eq_zero_iff.mp h_sq2
270 have h_eq3 : x 3 = 0 := sq_eq_zero_iff.mp h_sq3
271 have h_sn_eq : spatialNormSq x = s ^ 2 := by
272 show x 1 ^ 2 + x 2 ^ 2 + x 3 ^ 2 = s ^ 2
273 have h_x1 : x 1 = -s := by linarith
274 rw [h_x1, h_eq2, h_eq3]; ring
275 have h_r_sq_eq_s_sq : spatialRadius x ^ 2 = s ^ 2 := by rw [hr_sq, h_sn_eq]
276 -- r > 0 and r² = s², so r = √(s²) = |s|. Then |s| < r/2 = |s|/2 contradicts |s| ≥ 0 ∧ r > 0.
277 have h_r_eq_abs : spatialRadius x = |s| := by
278 have h_pos : 0 ≤ spatialRadius x := le_of_lt hr_pos
279 have h_abs_pos : 0 ≤ |s| := abs_nonneg s
280 have h_eq : Real.sqrt (spatialRadius x ^ 2) = Real.sqrt (s ^ 2) := by
281 rw [h_r_sq_eq_s_sq]
282 rw [Real.sqrt_sq h_pos] at h_eq
283 rw [show s ^ 2 = |s| ^ 2 from (sq_abs s).symm, Real.sqrt_sq h_abs_pos] at h_eq
284 exact h_eq
285 rw [h_r_eq_abs] at hs
286 linarith [abs_nonneg s]
287 | 2, _ =>
288 rw [show ((⟨2, by decide⟩ : Fin 4) = (2 : Fin 4)) from rfl,
289 spatialNormSq_coordRay_spatial_2]
290 intro h_zero
291 have h_sq1 : x 1 ^ 2 = 0 := by
292 nlinarith [sq_nonneg (x 1), sq_nonneg (x 2 + s), sq_nonneg (x 3)]
293 have h_sq2 : (x 2 + s) ^ 2 = 0 := by
294 nlinarith [sq_nonneg (x 1), sq_nonneg (x 2 + s), sq_nonneg (x 3)]
295 have h_sq3 : x 3 ^ 2 = 0 := by
296 nlinarith [sq_nonneg (x 1), sq_nonneg (x 2 + s), sq_nonneg (x 3)]
297 have h_eq1 : x 1 = 0 := sq_eq_zero_iff.mp h_sq1
298 have h_eq2 : x 2 + s = 0 := sq_eq_zero_iff.mp h_sq2
299 have h_eq3 : x 3 = 0 := sq_eq_zero_iff.mp h_sq3
300 have h_sn_eq : spatialNormSq x = s ^ 2 := by
301 show x 1 ^ 2 + x 2 ^ 2 + x 3 ^ 2 = s ^ 2
302 have h_x2 : x 2 = -s := by linarith
303 rw [h_eq1, h_x2, h_eq3]; ring
304 have h_r_sq_eq_s_sq : spatialRadius x ^ 2 = s ^ 2 := by rw [hr_sq, h_sn_eq]
305 have h_r_eq_abs : spatialRadius x = |s| := by
306 have h_pos : 0 ≤ spatialRadius x := le_of_lt hr_pos
307 have h_abs_pos : 0 ≤ |s| := abs_nonneg s
308 have h_eq : Real.sqrt (spatialRadius x ^ 2) = Real.sqrt (s ^ 2) := by
309 rw [h_r_sq_eq_s_sq]
310 rw [Real.sqrt_sq h_pos] at h_eq
311 rw [show s ^ 2 = |s| ^ 2 from (sq_abs s).symm, Real.sqrt_sq h_abs_pos] at h_eq
312 exact h_eq
313 rw [h_r_eq_abs] at hs
314 linarith [abs_nonneg s]
315 | 3, _ =>
316 rw [show ((⟨3, by decide⟩ : Fin 4) = (3 : Fin 4)) from rfl,
317 spatialNormSq_coordRay_spatial_3]
318 intro h_zero
319 have h_sq1 : x 1 ^ 2 = 0 := by
320 nlinarith [sq_nonneg (x 1), sq_nonneg (x 2), sq_nonneg (x 3 + s)]
321 have h_sq2 : x 2 ^ 2 = 0 := by
322 nlinarith [sq_nonneg (x 1), sq_nonneg (x 2), sq_nonneg (x 3 + s)]
323 have h_sq3 : (x 3 + s) ^ 2 = 0 := by
324 nlinarith [sq_nonneg (x 1), sq_nonneg (x 2), sq_nonneg (x 3 + s)]
325 have h_eq1 : x 1 = 0 := sq_eq_zero_iff.mp h_sq1
326 have h_eq2 : x 2 = 0 := sq_eq_zero_iff.mp h_sq2
327 have h_eq3 : x 3 + s = 0 := sq_eq_zero_iff.mp h_sq3
328 have h_sn_eq : spatialNormSq x = s ^ 2 := by
329 show x 1 ^ 2 + x 2 ^ 2 + x 3 ^ 2 = s ^ 2
330 have h_x3 : x 3 = -s := by linarith
331 rw [h_eq1, h_eq2, h_x3]; ring
332 have h_r_sq_eq_s_sq : spatialRadius x ^ 2 = s ^ 2 := by rw [hr_sq, h_sn_eq]
333 have h_r_eq_abs : spatialRadius x = |s| := by
334 have h_pos : 0 ≤ spatialRadius x := le_of_lt hr_pos
335 have h_abs_pos : 0 ≤ |s| := abs_nonneg s
336 have h_eq : Real.sqrt (spatialRadius x ^ 2) = Real.sqrt (s ^ 2) := by
337 rw [h_r_sq_eq_s_sq]
338 rw [Real.sqrt_sq h_pos] at h_eq
339 rw [show s ^ 2 = |s| ^ 2 from (sq_abs s).symm, Real.sqrt_sq h_abs_pos] at h_eq
340 exact h_eq
341 rw [h_r_eq_abs] at hs
342 linarith [abs_nonneg s]
343
344/-- Radial inverse function `1/r^n` where r is the spatial radius.
345 Used for gravitational potentials. -/
346noncomputable def radialInv (n : ℕ) (x : Fin 4 → ℝ) : ℝ :=
347 1 / (spatialRadius x) ^ n
348
349/-- Differentiability of a coordinate ray component. -/
350theorem differentiableAt_coordRay_i (x : Fin 4 → ℝ) (μ i : Fin 4) :
351 DifferentiableAt ℝ (fun t => (coordRay x μ t) i) 0 := by
352 simp only [coordRay_apply]
353 apply DifferentiableAt.add
354 · apply differentiableAt_const
355 · apply DifferentiableAt.mul
356 · apply differentiableAt_id
357 · apply differentiableAt_const
358
359/-- Differentiability of a squared coordinate ray component. -/
360theorem differentiableAt_coordRay_i_sq (x : Fin 4 → ℝ) (μ i : Fin 4) :
361 DifferentiableAt ℝ (fun t => (coordRay x μ t) i ^ 2) 0 := by
362 apply DifferentiableAt.pow (differentiableAt_coordRay_i x μ i) 2
363
364/-- Closed form for ∂μ (xᵢ²). -/
365theorem partialDeriv_v2_x_sq (μ i : Fin 4) (x : Fin 4 → ℝ) :
366 partialDeriv_v2 (fun y => y i ^ 2) μ x = 2 * x i * (if i = μ then 1 else 0) := by
367 unfold partialDeriv_v2
368 simp only [coordRay_apply]
369 let f_i := fun t => x i + t * basisVec μ i
370 have h_f : DifferentiableAt ℝ f_i 0 := differentiableAt_coordRay_i x μ i
371 rw [show (fun t => (x i + t * basisVec μ i) ^ 2) = f_i ^ 2 by rfl]
372 rw [deriv_pow h_f 2]
373 simp only [f_i]
374 split_ifs with h_eq
375 · subst h_eq
376 simp only [basisVec_self, mul_one]
377 rw [deriv_const_add, deriv_id'']
378 ring
379 · simp only [basisVec_ne h_eq, mul_zero, add_zero]
380 rw [deriv_const]
381 ring
382
383theorem deriv_coordRay_i (x : Fin 4 → ℝ) (i : Fin 4) :
384 deriv (fun t => (coordRay x i t) i) 0 = 1 := by
385 simp only [coordRay_apply, basisVec_self, mul_one]
386 rw [deriv_const_add, deriv_id'']
387
388theorem deriv_coordRay_j (x : Fin 4 → ℝ) (i j : Fin 4) (h : j ≠ i) :
389 deriv (fun t => (coordRay x i t) j) 0 = 0 := by
390 simp only [coordRay_apply, basisVec_ne h, mul_zero, add_zero]
391 exact deriv_const 0 (x j)
392
393/-- **THEOREM**: Functional derivative of spatialNormSq.
394 ∂_μ (∑ x_i²) = 2 x_μ for μ ∈ {1,2,3}, else 0.
395
396 **Derivation**: Using the chain rule and ∂_μ(x_i²) = 2x_i δ_{iμ}, we get:
397 ∂_μ(x₁² + x₂² + x₃²) = 2x₁δ_{1μ} + 2x₂δ_{2μ} + 2x₃δ_{3μ} = 2x_μ for μ ∈ {1,2,3}. -/
398theorem partialDeriv_v2_spatialNormSq (μ : Fin 4) (x : Fin 4 → ℝ) :
399 partialDeriv_v2 spatialNormSq μ x =
400 if μ = 0 then 0 else 2 * x μ := by
401 -- Each component x_i² gives 2x_i δ_{iμ}
402 have hd1 := partialDeriv_v2_x_sq μ 1 x
403 have hd2 := partialDeriv_v2_x_sq μ 2 x
404 have hd3 := partialDeriv_v2_x_sq μ 3 x
405 -- Enumerate all 4 cases for μ
406 fin_cases μ <;> simp_all [partialDeriv_v2, spatialNormSq, coordRay_apply, basisVec, deriv_const_add]
407
408/-- Differentiability of spatialNormSq along a coordinate ray. -/
409theorem differentiableAt_coordRay_spatialNormSq (x : Fin 4 → ℝ) (μ : Fin 4) :
410 DifferentiableAt ℝ (fun t => spatialNormSq (coordRay x μ t)) 0 := by
411 unfold spatialNormSq
412 apply DifferentiableAt.add
413 · apply DifferentiableAt.add
414 · exact differentiableAt_coordRay_i_sq x μ 1
415 · exact differentiableAt_coordRay_i_sq x μ 2
416 · exact differentiableAt_coordRay_i_sq x μ 3
417
418/-- Differentiability of spatialRadius along a coordinate ray. -/
419theorem differentiableAt_coordRay_spatialRadius (x : Fin 4 → ℝ) (μ : Fin 4) (hx : spatialRadius x ≠ 0) :
420 DifferentiableAt ℝ (fun t => spatialRadius (coordRay x μ t)) 0 := by
421 unfold spatialRadius
422 have h_sn_ne_zero : spatialNormSq (coordRay x μ 0) ≠ 0 := by
423 simp only [coordRay_zero]
424 exact (spatialRadius_ne_zero_iff x).mp hx
425 apply DifferentiableAt.sqrt (differentiableAt_coordRay_spatialNormSq x μ) h_sn_ne_zero
426
427/-- Differentiability of radialInv along a coordinate ray. -/
428theorem differentiableAt_coordRay_radialInv (n : ℕ) (x : Fin 4 → ℝ) (μ : Fin 4) (hx : spatialRadius x ≠ 0) :
429 DifferentiableAt ℝ (fun t => radialInv n (coordRay x μ t)) 0 := by
430 unfold radialInv
431 apply DifferentiableAt.div (differentiableAt_const (1 : ℝ))
432 · apply DifferentiableAt.pow (differentiableAt_coordRay_spatialRadius x μ hx)
433 · have h_pos : 0 < spatialRadius x := by
434 unfold spatialRadius
435 apply Real.sqrt_pos.mpr
436 have h_nonneg := spatialNormSq_nonneg x
437 have h_ne_zero := (spatialRadius_ne_zero_iff x).mp hx
438 exact lt_of_le_of_ne h_nonneg h_ne_zero.symm
439 simp only [coordRay_zero]
440 exact (pow_pos h_pos n).ne'
441
442theorem spatialRadius_coordRay_ne_zero_eventually {x : Fin 4 → ℝ} (hx : spatialRadius x ≠ 0) (μ : Fin 4) :
443 ∀ᶠ t in 𝓝 0, spatialRadius (coordRay x μ t) ≠ 0 := by
444 have h_cont : Continuous (fun t => spatialRadius (coordRay x μ t)) := by
445 unfold spatialRadius spatialNormSq coordRay basisVec
446 fun_prop
447 apply h_cont.continuousAt.eventually_ne
448 simp [coordRay_zero, hx]
449
450/-- Directional derivative of `spatialRadius` in coordinates.
451
452 For temporal direction (μ = 0), the spatial radius is invariant along the
453 coordinate ray, so the derivative is 0. For a spatial direction (μ ≠ 0),
454 we compose the chain rule for `Real.sqrt` (Mathlib `HasDerivAt.sqrt`)
455 with the derivative `∂_μ ‖x‖² = 2 x_μ` (lifted from
456 `partialDeriv_v2_spatialNormSq`), giving `(2 x_μ) / (2 √‖x‖²) = x_μ / r`.
457
458 Closes one of the seven §XXIII.B′ Mathlib calculus axioms. -/
459theorem partialDeriv_v2_spatialRadius (μ : Fin 4) (x : Fin 4 → ℝ) (hx : spatialRadius x ≠ 0) :
460 partialDeriv_v2 spatialRadius μ x =
461 if μ = 0 then 0 else x μ / spatialRadius x := by
462 by_cases hμ : μ = 0
463 · -- Temporal direction: `spatialRadius` is invariant along `coordRay x 0 _`.
464 simp only [hμ, ↓reduceIte]
465 unfold partialDeriv_v2
466 have h : ∀ t, spatialRadius (coordRay x 0 t) = spatialRadius x :=
467 spatialRadius_coordRay_temporal x
468 simp_rw [h]
469 exact deriv_const 0 _
470 · -- Spatial direction: chain rule for `Real.sqrt`.
471 simp only [hμ, ↓reduceIte]
472 unfold partialDeriv_v2 spatialRadius
473 -- `spatialNormSq x ≠ 0` since `spatialRadius x ≠ 0`.
474 have h_sn_ne : spatialNormSq x ≠ 0 := (spatialRadius_ne_zero_iff x).mp hx
475 -- Differentiability of `t ↦ ‖coordRay x μ t‖²` at 0.
476 have h_sn_da : DifferentiableAt ℝ (fun t => spatialNormSq (coordRay x μ t)) 0 :=
477 differentiableAt_coordRay_spatialNormSq x μ
478 -- Its derivative at 0 is `2 x_μ` (from `partialDeriv_v2_spatialNormSq`).
479 have h_sn_deriv : deriv (fun t => spatialNormSq (coordRay x μ t)) 0 = 2 * x μ := by
480 have := partialDeriv_v2_spatialNormSq μ x
481 simp only [partialDeriv_v2, hμ, ↓reduceIte] at this
482 exact this
483 -- Lift to `HasDerivAt`.
484 have h_sn_hda : HasDerivAt (fun t => spatialNormSq (coordRay x μ t)) (2 * x μ) 0 := by
485 have : HasDerivAt (fun t => spatialNormSq (coordRay x μ t))
486 (deriv (fun t => spatialNormSq (coordRay x μ t)) 0) 0 := h_sn_da.hasDerivAt
487 simpa [h_sn_deriv] using this
488 -- Value at 0 needed for `HasDerivAt.sqrt`.
489 have h_sn_eq_at_0 : spatialNormSq (coordRay x μ 0) = spatialNormSq x := by
490 simp [coordRay_zero]
491 -- Apply Mathlib's chain rule for sqrt.
492 have h_sqrt_hda : HasDerivAt (fun t => Real.sqrt (spatialNormSq (coordRay x μ t)))
493 (2 * x μ / (2 * Real.sqrt (spatialNormSq (coordRay x μ 0)))) 0 :=
494 h_sn_hda.sqrt (by rw [h_sn_eq_at_0]; exact h_sn_ne)
495 rw [h_sqrt_hda.deriv, h_sn_eq_at_0]
496 -- Goal: `2 x_μ / (2 √‖x‖²) = x_μ / √‖x‖²`.
497 have hr_ne : Real.sqrt (spatialNormSq x) ≠ 0 :=
498 fun h => h_sn_ne (Real.sqrt_eq_zero (spatialNormSq_nonneg x) |>.mp h)
499 field_simp
500
501/-- Directional derivative of `radialInv` in coordinates.
502
503 For temporal direction (μ = 0), `radialInv` is invariant along the ray,
504 so the derivative is 0. For a spatial direction (μ ≠ 0), we use the quotient
505 rule `HasDerivAt.div` on `1 / r^n`, composing with `HasDerivAt.pow` and
506 the derivative `∂_μ r = x_μ / r` (lifted from `partialDeriv_v2_spatialRadius`).
507
508 Closes one of the six remaining §XXIII.B′ Mathlib calculus axioms. -/
509theorem partialDeriv_v2_radialInv (n : ℕ) (μ : Fin 4) (x : Fin 4 → ℝ) (hx : spatialRadius x ≠ 0) :
510 partialDeriv_v2 (radialInv n) μ x =
511 if μ = 0 then 0 else - (n : ℝ) * x μ / (spatialRadius x) ^ (n + 2) := by
512 unfold partialDeriv_v2 radialInv
513 by_cases hμ : μ = 0
514 · simp only [hμ, ↓reduceIte]
515 have h : ∀ t, spatialRadius (coordRay x 0 t) = spatialRadius x :=
516 spatialRadius_coordRay_temporal x
517 have h2 : (fun t => 1 / spatialRadius (coordRay x 0 t) ^ n) =
518 (fun _ => 1 / spatialRadius x ^ n) := by
519 funext t; rw [h]
520 simp_rw [h2]; exact deriv_const 0 _
521 · simp only [hμ, ↓reduceIte]
522 cases n with
523 | zero => simp
524 | succ k =>
525 have hr_pos : 0 < spatialRadius x := spatialRadius_pos_of_ne_zero x hx
526 have h_r_da : DifferentiableAt ℝ (fun t => spatialRadius (coordRay x μ t)) 0 :=
527 differentiableAt_coordRay_spatialRadius x μ hx
528 have h_r_pow_da : DifferentiableAt ℝ (fun t => spatialRadius (coordRay x μ t) ^ (k + 1)) 0 :=
529 h_r_da.pow (k + 1)
530 have h_r_pow_ne : spatialRadius (coordRay x μ 0) ^ (k + 1) ≠ 0 := by
531 simp only [coordRay_zero]
532 exact pow_ne_zero (k + 1) hx
533 have h_r_deriv : deriv (fun t => spatialRadius (coordRay x μ t)) 0 = x μ / spatialRadius x := by
534 have := partialDeriv_v2_spatialRadius μ x hx
535 simp only [partialDeriv_v2, hμ, ↓reduceIte] at this
536 exact this
537 have h_r_hda : HasDerivAt (fun t => spatialRadius (coordRay x μ t)) (x μ / spatialRadius x) 0 := by
538 have : HasDerivAt (fun t => spatialRadius (coordRay x μ t))
539 (deriv (fun t => spatialRadius (coordRay x μ t)) 0) 0 := h_r_da.hasDerivAt
540 simpa [h_r_deriv] using this
541 have h_rpow_hda : HasDerivAt (fun t => spatialRadius (coordRay x μ t) ^ (k + 1))
542 (↑(k + 1) * spatialRadius x ^ k * (x μ / spatialRadius x)) 0 := by
543 have h1 := h_r_hda.pow (k + 1)
544 simp only [coordRay_zero] at h1
545 convert h1 using 2
546 have h_rinv_hda : HasDerivAt (fun t => (spatialRadius (coordRay x μ t) ^ (k + 1))⁻¹)
547 (-(↑(k + 1) * spatialRadius x ^ k * (x μ / spatialRadius x)) /
548 (spatialRadius x ^ (k + 1)) ^ 2) 0 := by
549 have h2 := h_rpow_hda.inv h_r_pow_ne
550 simp only [coordRay_zero] at h2
551 exact h2
552 have h_inv : (fun t => 1 / spatialRadius (coordRay x μ t) ^ (k + 1)) = fun t => (spatialRadius (coordRay x μ t) ^ (k + 1))⁻¹ := by funext t; exact one_div _
553 rw [h_inv, h_rinv_hda.deriv]
554 have hr_ne : spatialRadius x ≠ 0 := hx
555 have h_pow1 : (spatialRadius x ^ (k + 1)) ^ 2 = spatialRadius x ^ (2 * k + 2) := by
556 rw [← pow_mul]; congr 1; omega
557 have h_pow2 : k + 1 + 2 = k + 3 := by omega
558 rw [div_eq_div_iff]
559 · change -(↑(k + 1) * spatialRadius x ^ k * (x μ / spatialRadius x)) * spatialRadius x ^ (k + 1 + 2) = -↑(k + 1) * x μ * (spatialRadius x ^ (k + 1)) ^ 2
560 rw [h_pow1]
561 rw [h_pow2]
562 calc -(↑(k + 1) * spatialRadius x ^ k * (x μ / spatialRadius x)) * spatialRadius x ^ (k + 3)
563 _ = -(↑(k + 1) * spatialRadius x ^ k * (x μ * (spatialRadius x)⁻¹)) * spatialRadius x ^ (k + 3) := by rw [div_eq_mul_inv]
564 _ = -↑(k + 1) * x μ * (spatialRadius x ^ k * (spatialRadius x)⁻¹ * spatialRadius x ^ (k + 3)) := by ring
565 _ = -↑(k + 1) * x μ * (spatialRadius x ^ k * spatialRadius x ^ (k + 3) * (spatialRadius x)⁻¹) := by ring
566 _ = -↑(k + 1) * x μ * (spatialRadius x ^ (2 * k + 3) * (spatialRadius x)⁻¹) := by
567 have : spatialRadius x ^ k * spatialRadius x ^ (k + 3) = spatialRadius x ^ (2 * k + 3) := by rw [← pow_add]; congr 1; omega
568 rw [this]
569 _ = -↑(k + 1) * x μ * (spatialRadius x ^ (2 * k + 2) * spatialRadius x * (spatialRadius x)⁻¹) := by
570 have : spatialRadius x ^ (2 * k + 2) * spatialRadius x = spatialRadius x ^ (2 * k + 3) := by
571 calc spatialRadius x ^ (2 * k + 2) * spatialRadius x
572 _ = spatialRadius x ^ (2 * k + 2) * spatialRadius x ^ 1 := by rw [pow_one]
573 _ = spatialRadius x ^ (2 * k + 2 + 1) := by rw [← pow_add]
574 _ = spatialRadius x ^ (2 * k + 3) := by rfl
575 rw [← this]
576 _ = -↑(k + 1) * x μ * (spatialRadius x ^ (2 * k + 2) * (spatialRadius x * (spatialRadius x)⁻¹)) := by ring
577 _ = -↑(k + 1) * x μ * (spatialRadius x ^ (2 * k + 2) * 1) := by rw [mul_inv_cancel₀ hr_ne]
578 _ = -↑(k + 1) * x μ * spatialRadius x ^ (2 * k + 2) := by ring
579 · exact pow_ne_zero 2 (pow_ne_zero (k + 1) hr_ne)
580 · exact pow_ne_zero (k + 1 + 2) hr_ne
581
582/-- Differentiability of `s ↦ ∂(1/r^n)/∂x_μ` along a coordinate ray.
583
584 For temporal direction (μ = 0), `partialDeriv_v2 (radialInv n) 0 y = 0` whenever
585 `spatialRadius y ≠ 0`, so the function is locally constant 0 near `s = 0`.
586
587 For spatial direction (μ ≠ 0), `partialDeriv_v2 (radialInv n) μ y =
588 -n · y_μ / r(y)^(n+2)` whenever `spatialRadius y ≠ 0`, so the function is
589 locally a quotient of differentiable functions:
590 - numerator `s ↦ -n · (coordRay x ν s)_μ` is linear in `s`
591 - denominator `s ↦ r(coordRay x ν s)^(n+2)` is the `(n+2)`-power of the
592 already-proved differentiable `s ↦ r(coordRay x ν s)`
593 - denominator is nonzero at `s = 0` since `r(x) ≠ 0`.
594
595 Closes one of the §XXIII.B′ Mathlib calculus axioms. -/
596theorem differentiableAt_coordRay_partialDeriv_v2_radialInv
597 (n : ℕ) (x : Fin 4 → ℝ) (μ ν : Fin 4) (hx : spatialRadius x ≠ 0) :
598 DifferentiableAt ℝ (fun s => partialDeriv_v2 (radialInv n) μ (coordRay x ν s)) 0 := by
599 by_cases hμ : μ = 0
600 · -- Temporal: `partialDeriv_v2 (radialInv n) 0 y = 0` for `r y ≠ 0`,
601 -- so the function is eventually 0 near `s = 0`.
602 apply (differentiableAt_const (0 : ℝ)).congr_of_eventuallyEq
603 apply (spatialRadius_coordRay_ne_zero_eventually hx ν).mono
604 intro s hs
605 show partialDeriv_v2 (radialInv n) μ (coordRay x ν s) = 0
606 rw [partialDeriv_v2_radialInv n μ (coordRay x ν s) hs]
607 simp [hμ]
608 · -- Spatial: smooth quotient formula `-n · x_μ / r^(n+2)` is differentiable,
609 -- and `partialDeriv_v2` agrees with it where `r ≠ 0`.
610 have h_num_diff : DifferentiableAt ℝ
611 (fun s : ℝ => -(n : ℝ) * (coordRay x ν s) μ) 0 :=
612 (differentiableAt_coordRay_i x ν μ).const_mul (-(n : ℝ))
613 have h_denom_diff : DifferentiableAt ℝ
614 (fun s : ℝ => spatialRadius (coordRay x ν s) ^ (n + 2)) 0 :=
615 (differentiableAt_coordRay_spatialRadius x ν hx).pow _
616 have h_denom_ne : spatialRadius (coordRay x ν 0) ^ (n + 2) ≠ 0 := by
617 simp only [coordRay_zero]; exact pow_ne_zero (n + 2) hx
618 have h_smooth_diff : DifferentiableAt ℝ
619 (fun s : ℝ => -(n : ℝ) * (coordRay x ν s) μ /
620 (spatialRadius (coordRay x ν s)) ^ (n + 2)) 0 :=
621 h_num_diff.div h_denom_diff h_denom_ne
622 apply h_smooth_diff.congr_of_eventuallyEq
623 apply (spatialRadius_coordRay_ne_zero_eventually hx ν).mono
624 intro s hs
625 show partialDeriv_v2 (radialInv n) μ (coordRay x ν s) =
626 -(n : ℝ) * (coordRay x ν s) μ / (spatialRadius (coordRay x ν s)) ^ (n + 2)
627 rw [partialDeriv_v2_radialInv n μ (coordRay x ν s) hs]
628 simp [hμ]
629
630/-- Second directional derivative of `radialInv n`:
631 `∂_ν ∂_μ (1/r^n) = n · ((n+2) x_μ x_ν / r^{n+4} - δ_{μν} / r^{n+2})` for `μ, ν ∈ {1,2,3}`,
632 and `0` whenever either index is `0` (temporal).
633
634 Proof structure:
635 - `μ = 0`: `partialDeriv_v2 (radialInv n) 0 = 0` (already proved), so the outer derivative
636 is `deriv 0 = 0`.
637 - `μ ≠ 0, ν = 0`: spatial partial derivative is invariant along the temporal ray
638 (since both `x_μ` and `r` are unchanged), so the outer derivative is `deriv const = 0`.
639 - `μ, ν ≠ 0`: apply the quotient rule `HasDerivAt.div` to the smooth formula
640 `-n · x_μ / r^{n+2}`, then simplify the resulting algebraic expression
641 `(f' g - f g') / g²` to match the target. The simplification uses
642 `r^{n+2-1} = r^{n+1}` (natural-number subtraction), `(r^{n+2})² = r^{2(n+2)} = r^{n+4} · r^n`,
643 and `field_simp` plus `ring` to clear denominators.
644
645 Closes one of the §XXIII.B′ Mathlib calculus axioms. -/
646theorem secondDeriv_radialInv (n : ℕ) (μ ν : Fin 4) (x : Fin 4 → ℝ) (hx : spatialRadius x ≠ 0) :
647 secondDeriv (radialInv n) μ ν x =
648 if μ = 0 ∨ ν = 0 then 0 else
649 (n : ℝ) *
650 ((n + 2 : ℝ) * x μ * x ν / (spatialRadius x) ^ (n + 4) -
651 (if μ = ν then 1 else 0) / (spatialRadius x) ^ (n + 2)) := by
652 unfold secondDeriv
653 by_cases hμ : μ = 0
654 · -- μ = 0: outer function is eventually 0.
655 simp only [hμ, true_or, ↓reduceIte]
656 have h_ev : (fun s => partialDeriv_v2 (radialInv n) 0 (coordRay x ν s)) =ᶠ[𝓝 0]
657 (fun _ => (0 : ℝ)) := by
658 apply (spatialRadius_coordRay_ne_zero_eventually hx ν).mono
659 intro s hs
660 show partialDeriv_v2 (radialInv n) 0 (coordRay x ν s) = 0
661 rw [partialDeriv_v2_radialInv n 0 (coordRay x ν s) hs]
662 simp
663 rw [h_ev.deriv_eq]; exact deriv_const 0 _
664 · by_cases hν : ν = 0
665 · -- μ ≠ 0, ν = 0: outer function is constant in s.
666 simp only [hν, hμ, false_or, ↓reduceIte, or_true]
667 have h_const : ∀ s, partialDeriv_v2 (radialInv n) μ (coordRay x 0 s) =
668 partialDeriv_v2 (radialInv n) μ x := by
669 intro s
670 have hr_s : spatialRadius (coordRay x 0 s) = spatialRadius x :=
671 spatialRadius_coordRay_temporal x s
672 have h_coord : (coordRay x 0 s) μ = x μ := coordRay_temporal_spatial x s μ hμ
673 have hr_ne_s : spatialRadius (coordRay x 0 s) ≠ 0 := by rw [hr_s]; exact hx
674 rw [partialDeriv_v2_radialInv n μ (coordRay x 0 s) hr_ne_s,
675 partialDeriv_v2_radialInv n μ x hx]
676 simp only [hμ, ↓reduceIte]
677 rw [hr_s, h_coord]
678 simp_rw [h_const]
679 exact deriv_const 0 _
680 · -- Main case: μ ≠ 0 and ν ≠ 0.
681 simp only [hμ, hν, false_or, ↓reduceIte]
682 have hr_pos : 0 < spatialRadius x := spatialRadius_pos_of_ne_zero x hx
683 have hr_ne : spatialRadius x ≠ 0 := hx
684 -- Near 0, the outer function equals the smooth formula `-n · y_μ / r(y)^(n+2)`.
685 have h_ev : (fun s => partialDeriv_v2 (radialInv n) μ (coordRay x ν s)) =ᶠ[𝓝 0]
686 (fun s => -(n : ℝ) * (coordRay x ν s) μ /
687 (spatialRadius (coordRay x ν s)) ^ (n + 2)) := by
688 apply (spatialRadius_coordRay_ne_zero_eventually hx ν).mono
689 intro s hs
690 show partialDeriv_v2 (radialInv n) μ (coordRay x ν s) =
691 -(n : ℝ) * (coordRay x ν s) μ /
692 (spatialRadius (coordRay x ν s)) ^ (n + 2)
693 rw [partialDeriv_v2_radialInv n μ (coordRay x ν s) hs]
694 simp [hμ]
695 rw [h_ev.deriv_eq]
696 -- Compute deriv via HasDerivAt.div on the smooth formula.
697 -- Numerator: c(s) = -n · (x_μ + s · δ_{μ,ν}); c(0) = -n · x_μ; c'(0) = -n · δ_{μ,ν}.
698 have h_num_hda : HasDerivAt (fun s => -(n : ℝ) * (coordRay x ν s) μ)
699 (-(n : ℝ) * basisVec ν μ) 0 := by
700 simp only [coordRay_apply]
701 have h_inner : HasDerivAt (fun s => x μ + s * basisVec ν μ) (basisVec ν μ) 0 := by
702 have := ((hasDerivAt_id (0 : ℝ)).mul_const (basisVec ν μ)).const_add (x μ)
703 simpa using this
704 exact h_inner.const_mul (-(n : ℝ))
705 -- Denominator inner: r(coordRay x ν s); deriv at 0 is `x_ν / r`.
706 have h_r_da : DifferentiableAt ℝ (fun s => spatialRadius (coordRay x ν s)) 0 :=
707 differentiableAt_coordRay_spatialRadius x ν hx
708 have h_r_deriv : deriv (fun s => spatialRadius (coordRay x ν s)) 0 =
709 x ν / spatialRadius x := by
710 have := partialDeriv_v2_spatialRadius ν x hx
711 simp only [partialDeriv_v2, hν, ↓reduceIte] at this
712 exact this
713 have h_r_hda : HasDerivAt (fun s => spatialRadius (coordRay x ν s))
714 (x ν / spatialRadius x) 0 := by
715 have := h_r_da.hasDerivAt
716 simpa [h_r_deriv] using this
717 -- Denominator: g(s) = r(coordRay x ν s)^(n+2); g(0) = r^(n+2);
718 -- g'(0) = (n+2) · r^(n+1) · (x_ν / r) (note `n + 2 - 1 = n + 1` via Nat subtraction).
719 have h_den_hda : HasDerivAt (fun s => (spatialRadius (coordRay x ν s)) ^ (n + 2))
720 ((n + 2 : ℕ) * spatialRadius x ^ (n + 2 - 1) * (x ν / spatialRadius x)) 0 := by
721 have h := h_r_hda.pow (n + 2)
722 simp only [coordRay_zero] at h
723 exact h
724 have h_den_ne : (spatialRadius (coordRay x ν 0)) ^ (n + 2) ≠ 0 := by
725 simp only [coordRay_zero]; exact pow_ne_zero (n + 2) hr_ne
726 -- Apply the quotient rule.
727 have h_quot : HasDerivAt
728 (fun s => -(n : ℝ) * (coordRay x ν s) μ /
729 (spatialRadius (coordRay x ν s)) ^ (n + 2))
730 ((-(n : ℝ) * basisVec ν μ * spatialRadius x ^ (n + 2) -
731 -(n : ℝ) * x μ *
732 (↑(n + 2) * spatialRadius x ^ (n + 2 - 1) * (x ν / spatialRadius x))) /
733 (spatialRadius x ^ (n + 2)) ^ 2) 0 := by
734 have h := h_num_hda.div h_den_hda h_den_ne
735 simp only [coordRay_zero] at h
736 exact h
737 rw [h_quot.deriv]
738 -- Algebraic finish.
739 -- Step 1: Reduce `r^(n+2-1)` (Nat subtraction) to `r^(n+1)`.
740 have h_n_sub : n + 2 - 1 = n + 1 := by omega
741 rw [h_n_sub]
742 -- Step 2: Identify `basisVec ν μ` with `if μ = ν then 1 else 0` (definitional).
743 have h_basisVec : basisVec ν μ = if μ = ν then 1 else 0 := rfl
744 rw [h_basisVec]
745 -- Step 3: Prepare power identities for the final ring step.
746 have h_r_pow_ne : ∀ k, spatialRadius x ^ k ≠ 0 := fun k => pow_ne_zero k hr_ne
747 have h_pow_succ_1 : spatialRadius x ^ (n + 1) = spatialRadius x ^ n * spatialRadius x := by
748 rw [pow_succ]
749 have h_pow_succ_2 : spatialRadius x ^ (n + 2) =
750 spatialRadius x ^ n * spatialRadius x * spatialRadius x := by
751 rw [show n + 2 = (n + 1) + 1 from rfl, pow_succ, h_pow_succ_1]
752 have h_pow_n_plus_4 : spatialRadius x ^ (n + 4) =
753 spatialRadius x ^ n * spatialRadius x * spatialRadius x *
754 spatialRadius x * spatialRadius x := by
755 rw [show n + 4 = ((n + 2) + 1) + 1 from rfl, pow_succ, pow_succ, h_pow_succ_2]
756 -- Step 4: Clear denominators (field_simp), then use power identities and ring.
757 field_simp
758 rw [h_pow_succ_1, h_pow_succ_2, h_pow_n_plus_4]
759 push_cast
760 ring
761
762/-- Laplacian of `1/r` vanishes (the radial inverse is harmonic away from the origin).
763
764 `∇²(1/r) = ∂₁²(1/r) + ∂₂²(1/r) + ∂₃²(1/r)`.
765 Each diagonal second-derivative is `3 x_i² / r^5 - 1/r^3` from `secondDeriv_radialInv` at `n = 1`.
766 The sum is `3 ‖x‖² / r^5 - 3/r^3 = 3 r² / r^5 - 3/r^3 = 3/r^3 - 3/r^3 = 0`.
767
768 Closes one of the §XXIII.B′ Mathlib calculus axioms. -/
769theorem laplacian_radialInv_zero_no_const
770 (x : Fin 4 → ℝ) (hx : spatialRadius x ≠ 0) :
771 laplacian (radialInv 1) x = 0 := by
772 unfold laplacian
773 have hr_pos : 0 < spatialRadius x := spatialRadius_pos_of_ne_zero x hx
774 -- Use `(1 : Fin 4)`, `(2 : Fin 4)`, `(3 : Fin 4)` to match `spatialNormSq` syntactically.
775 have h1 := secondDeriv_radialInv 1 (1 : Fin 4) (1 : Fin 4) x hx
776 have h2 := secondDeriv_radialInv 1 (2 : Fin 4) (2 : Fin 4) x hx
777 have h3 := secondDeriv_radialInv 1 (3 : Fin 4) (3 : Fin 4) x hx
778 simp only [show ((1 : Fin 4) ≠ 0) from by decide,
779 show ((2 : Fin 4) ≠ 0) from by decide,
780 show ((3 : Fin 4) ≠ 0) from by decide,
781 or_self, ↓reduceIte] at h1 h2 h3
782 -- The `laplacian` definition uses `⟨i, _⟩ : Fin 4` form, but these are defeq to `(i : Fin 4)`.
783 -- Convert via `show` to align indices.
784 show secondDeriv (radialInv 1) (1 : Fin 4) (1 : Fin 4) x +
785 secondDeriv (radialInv 1) (2 : Fin 4) (2 : Fin 4) x +
786 secondDeriv (radialInv 1) (3 : Fin 4) (3 : Fin 4) x = 0
787 rw [h1, h2, h3]
788 -- Algebraic finish: 3 (x₁² + x₂² + x₃²)/r^5 - 3/r^3 = 0 since x₁² + x₂² + x₃² = r².
789 have hr_sq : spatialRadius x ^ 2 = x 1 ^ 2 + x 2 ^ 2 + x 3 ^ 2 := by
790 unfold spatialRadius
791 rw [Real.sq_sqrt (spatialNormSq_nonneg x)]; rfl
792 have h_pow5 : spatialRadius x ^ 5 = spatialRadius x ^ 3 * spatialRadius x ^ 2 := by
793 rw [← pow_add]
794 have h_r3_ne : spatialRadius x ^ 3 ≠ 0 := pow_ne_zero 3 hx
795 have h_r5_ne : spatialRadius x ^ 5 ≠ 0 := pow_ne_zero 5 hx
796 field_simp
797 rw [h_pow5, hr_sq]
798 ring
799
800/-- Laplacian of `1/r^n` for general `n`:
801 `∇²(1/r^n) = n(n-1)/r^(n+2)` (away from the origin).
802
803 Sum of diagonal `secondDeriv_radialInv n i i x` over `i ∈ {1,2,3}`:
804 `n · ((n+2) ‖x‖² / r^(n+4) - 3/r^(n+2)) = n · ((n+2)/r^(n+2) - 3/r^(n+2))`
805 `= n(n-1)/r^(n+2)` (using `‖x‖² = r²`).
806
807 Closes one of the §XXIII.B′ Mathlib calculus axioms. -/
808theorem laplacian_radialInv_n
809 (n : ℕ) (x : Fin 4 → ℝ) (hx : spatialRadius x ≠ 0) :
810 laplacian (radialInv n) x =
811 (n : ℝ) * ((n : ℝ) - 1) / (spatialRadius x) ^ (n + 2) := by
812 unfold laplacian
813 have hr_pos : 0 < spatialRadius x := spatialRadius_pos_of_ne_zero x hx
814 have h1 := secondDeriv_radialInv n (1 : Fin 4) (1 : Fin 4) x hx
815 have h2 := secondDeriv_radialInv n (2 : Fin 4) (2 : Fin 4) x hx
816 have h3 := secondDeriv_radialInv n (3 : Fin 4) (3 : Fin 4) x hx
817 simp only [show ((1 : Fin 4) ≠ 0) from by decide,
818 show ((2 : Fin 4) ≠ 0) from by decide,
819 show ((3 : Fin 4) ≠ 0) from by decide,
820 or_self, ↓reduceIte] at h1 h2 h3
821 show secondDeriv (radialInv n) (1 : Fin 4) (1 : Fin 4) x +
822 secondDeriv (radialInv n) (2 : Fin 4) (2 : Fin 4) x +
823 secondDeriv (radialInv n) (3 : Fin 4) (3 : Fin 4) x =
824 (n : ℝ) * ((n : ℝ) - 1) / spatialRadius x ^ (n + 2)
825 rw [h1, h2, h3]
826 -- Algebraic finish via linear_combination using `r² = x₁² + x₂² + x₃²`.
827 have hr_sq : spatialRadius x ^ 2 = x 1 ^ 2 + x 2 ^ 2 + x 3 ^ 2 := by
828 unfold spatialRadius
829 rw [Real.sq_sqrt (spatialNormSq_nonneg x)]; rfl
830 have h_pow_n_plus_4 : spatialRadius x ^ (n + 4) =
831 spatialRadius x ^ (n + 2) * spatialRadius x ^ 2 := by
832 rw [← pow_add]
833 have h_r_n2_ne : spatialRadius x ^ (n + 2) ≠ 0 := pow_ne_zero _ hx
834 have h_r_n4_ne : spatialRadius x ^ (n + 4) ≠ 0 := pow_ne_zero _ hx
835 field_simp
836 rw [h_pow_n_plus_4]
837 linear_combination
838 (-(n : ℝ) * ((n : ℝ) + 2) * spatialRadius x ^ (n + 2)) * hr_sq
839