DiscreteGreen.lean

1 import Mathlib.LinearAlgebra.Matrix.NonsingularInverse
2 import Mathlib.Data.Matrix.Basic
3 import Mathlib.Tactic
4 
5 /-! Author: Komorebi P. The literal Dirichlet tridiagonal inverse at every grid size. -/
6 namespace KomorebiPaperVI.DiscreteGreen
7 open scoped BigOperators
8 set_option backward.isDefEq.respectTransparency false
9 
10 noncomputable def green (N : ℝ) (i j : ℕ) : ℝ :=
11   min (i:ℝ) (j:ℝ)-(i:ℝ)*(j:ℝ)/N
12 
13 theorem green_zero (N : ℝ) (j : ℕ) : green N 0 j=0 := by simp [green]
14 
15 theorem green_end (n j : ℕ) (hj : j  ≤  n+1) : green (n+1) (n+1) j=0 := by
16   have h : (j:ℝ)  ≤  n+1 := by exact_mod_cast hj
17   have hn : (n:ℝ)+1≠0 := by positivity
18   simp only [green,Nat.cast_add,Nat.cast_one,min_eq_right h]
19   field_simp
20   ring
21 
22 theorem green_second_difference (N : ℝ) (i j : ℕ) :
23     2*green N (i+1) j-green N i j-green N (i+2) j=if i+1=j then 1 else 0 := by
24   unfold green
25   by_cases ha : j  ≤  i
26   · have h0 : (j:ℝ) ≤ i := by exact_mod_cast ha
27     have h1 : (j:ℝ) ≤ (i+1:ℕ) := by exact_mod_cast (show j ≤ i+1 by omega)
28     have h2 : (j:ℝ) ≤ (i+2:ℕ) := by exact_mod_cast (show j ≤ i+2 by omega)
29     rw [min_eq_right h1,min_eq_right h0,min_eq_right h2,ite_eq_right (by omega)]
30     push_cast
31     ring
32   · by_cases hb : j=i+1
33     · subst j
34       have h0 : (i:ℝ) ≤ (i+1:ℕ) := by exact_mod_cast (show i ≤ i+1 by omega)
35       have h2 : ((i+1:ℕ):ℝ) ≤ (i+2:ℕ) := by exact_mod_cast (show i+1 ≤ i+2 by omega)
36       rw [min_self,min_eq_left h0,min_eq_right h2,ite_eq_left rfl]
37       push_cast
38       ring
39     · have h0 : (i:ℝ) ≤ j := by exact_mod_cast (show i ≤ j by omega)
40       have h1 : ((i+1:ℕ):ℝ) ≤ j := by exact_mod_cast (show i+1 ≤ j by omega)
41       have h2 : ((i+2:ℕ):ℝ) ≤ j := by exact_mod_cast (show i+2 ≤ j by omega)
42       rw [min_eq_left h1,min_eq_left h0,min_eq_left h2,ite_eq_right (by omega)]
43       push_cast
44       ring
45 
46 noncomputable def tridiagonal (n : ℕ) : Matrix (Fin n) (Fin n) ℝ := fun i j =>
47   (if i=j then 2 else 0)-(if j.val=i.val+1 then 1 else 0)-(if j.val+1=i.val then 1 else 0)
48 
49 noncomputable def inverseKernel (n : ℕ) : Matrix (Fin n) (Fin n) ℝ :=
50   fun i j => green (n+1) (i.val+1) (j.val+1)
51 
52 theorem sparse_sum (n q : ℕ) (v : Fin n → ℝ) :
53     (∑ k : Fin n, if k.val=q then v k else 0)=if h:q < n then v ⟨q,h⟩ else 0 := by
54   classical
55   split_ifs with h
56   · rw [Finset.sum_eq_single (⟨q,h⟩ : Fin n)]
57     · simp
58     · intro k _ hk
59       exact ite_eq_right (by intro he; apply hk; exact Fin.ext he)
60     · simp
61   · apply Finset.sum_eq_zero
62     intro k _
63     exact ite_eq_right (by omega)
64 
65 theorem tridiagonal_mul (n : ℕ) (M : Matrix (Fin n) (Fin n) ℝ) (i j : Fin n) :
66     (tridiagonal n*M) i j=2*M i j-
67       (if h:i.val+1 < n then M ⟨i.val+1,h⟩ j else 0)-
68       (if h:0 < i.val then M ⟨i.val-1,by omega⟩ j else 0) := by
69   classical
70   simp only [Matrix.mul_apply,tridiagonal,sub_mul,Finset.sum_sub_distrib,ite_mul,one_mul,zero_mul]
71   rw [Finset.sum_ite_eq,sparse_sum]
72   simp only [Finset.mem_univ,ite_true]
73   congr 1
74   by_cases h:0 < i.val
75   · rw [dite_eq_left h,Finset.sum_eq_single (⟨i.val-1,by omega⟩ : Fin n)]
76     · simp [show i.val-1+1=i.val by omega]
77     · intro k _ hk
78       exact ite_eq_right (by intro he; apply hk; apply Fin.ext; change k.val=i.val-1; omega)
79     · simp
80   · rw [dite_eq_right h]
81     apply Finset.sum_eq_zero
82     intro k _
83     exact ite_eq_right (by omega)
84 
85 theorem tridiagonal_inverse_kernel (n : ℕ) : tridiagonal n*inverseKernel n=1 := by
86   ext i j
87   rw [tridiagonal_mul]
88   have h := green_second_difference (n+1) i.val (j.val+1)
89   have hu : (if h:i.val+1 < n then inverseKernel n ⟨i.val+1,h⟩ j else 0)=
90       green (n+1) (i.val+2) (j.val+1) := by
91     split_ifs with hi
92     · rfl
93     · have he : i.val+2=n+1 := by omega
94       rw [he,green_end n (j.val+1) (by omega)]
95   have hd : (if h:0 < i.val then inverseKernel n ⟨i.val-1,by omega⟩ j else 0)=
96       green (n+1) i.val (j.val+1) := by
97     split_ifs with hi
98     · simp only [inverseKernel]
99       congr 1
100       omega
101     · have he : i.val=0 := by omega
102       rw [he,green_zero]
103   rw [hu,hd]
104   simp only [inverseKernel,Matrix.one_apply]
105   have he : (i.val+1=j.val+1) ↔ i=j := by constructor  <;> intro hh; exact Fin.ext (by omega); subst j; rfl
106   simp only [he] at h
107   linarith
108 
109 theorem inverse_eq_kernel (n : ℕ) : (tridiagonal n)⁻¹=inverseKernel n := by
110   exact Matrix.inv_eq_right_inv (tridiagonal_inverse_kernel n)
111 
112 noncomputable def laplacian (n : ℕ) : Matrix (Fin n) (Fin n) ℝ :=
113   ((n:ℝ)+1)^2 • tridiagonal n
114 
115 theorem laplacian_inverse (n : ℕ) :
116     (laplacian n)⁻¹=(((n:ℝ)+1)^2)⁻¹ • inverseKernel n := by
117   apply Matrix.inv_eq_right_inv
118   rw [laplacian,smul_mul_smul_comm,tridiagonal_inverse_kernel]
119   have hn : ((n:ℝ)+1)^2≠0 := by positivity
120   simp [hn]
121 
122 theorem grid_reconstruction (n : ℕ) (i j : Fin n) :
123     ((n:ℝ)+1)*(laplacian n)⁻¹ i j=
124       min (((i.val:ℝ)+1)/((n:ℝ)+1)) (((j.val:ℝ)+1)/((n:ℝ)+1))-
125         (((i.val:ℝ)+1)/((n:ℝ)+1))*(((j.val:ℝ)+1)/((n:ℝ)+1)) := by
126   rw [laplacian_inverse]
127   simp only [Matrix.smul_apply,smul_eq_mul,inverseKernel,green,Nat.cast_add,Nat.cast_one]
128   rw [min_div_div_right (by positivity)]
129   field_simp
130 
131 end KomorebiPaperVI.DiscreteGreen

Back to source controls