author Manuel Eberl Fri Jul 13 18:27:18 2018 +0100 (2 days ago) changeset 68625 2ec84498f562 parent 68624 205d352ed727 child 68628 958511f53de8
HOL-Analysis: Volume of a simplex, Heron's theorem
```     1.1 --- a/src/HOL/Analysis/Analysis.thy	Fri Jul 13 16:54:36 2018 +0100
1.2 +++ b/src/HOL/Analysis/Analysis.thy	Fri Jul 13 18:27:18 2018 +0100
1.3 @@ -25,6 +25,7 @@
1.4    Gamma_Function
1.5    Change_Of_Vars
1.6    Lipschitz
1.7 +  Simplex_Content
1.8  begin
1.9
1.10  end
```
```     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/HOL/Analysis/Simplex_Content.thy	Fri Jul 13 18:27:18 2018 +0100
2.3 @@ -0,0 +1,263 @@
2.4 +(*
2.5 +   File:      Analysis/Simplex_Content.thy
2.6 +   Author:    Manuel Eberl <eberlm@in.tum.de>
2.7 +
2.8 +   The content of an n-dimensional simplex, including the formula for the content of a
2.9 +   triangle and Heron's formula.
2.10 +*)
2.11 +section \<open>Volume of a simplex\<close>
2.12 +theory Simplex_Content
2.13 +  imports Change_Of_Vars
2.14 +begin
2.15 +
2.16 +lemma fact_neq_top_ennreal [simp]: "fact n \<noteq> (top :: ennreal)"
2.17 +  by (induction n) (auto simp: ennreal_mult_eq_top_iff)
2.18 +
2.19 +lemma ennreal_fact: "ennreal (fact n) = fact n"
2.20 +  by (induction n) (auto simp: ennreal_mult algebra_simps ennreal_of_nat_eq_real_of_nat)
2.21 +
2.22 +context
2.23 +  fixes S :: "'a set \<Rightarrow> real \<Rightarrow> ('a \<Rightarrow> real) set"
2.24 +  defines "S \<equiv> (\<lambda>A t. {x. (\<forall>i\<in>A. 0 \<le> x i) \<and> sum x A \<le> t})"
2.25 +begin
2.26 +
2.27 +lemma emeasure_std_simplex_aux_step:
2.28 +  assumes "b \<notin> A" "finite A"
2.29 +  shows   "x(b := y) \<in> S (insert b A) t \<longleftrightarrow> y \<in> {0..t} \<and> x \<in> S A (t - y)"
2.30 +  using assms sum_nonneg[of A x] unfolding S_def
2.31 +  by (force simp: sum_delta_notmem algebra_simps)
2.32 +
2.33 +lemma emeasure_std_simplex_aux:
2.34 +  fixes t :: real
2.35 +  assumes "finite (A :: 'a set)" "t \<ge> 0"
2.36 +  shows   "emeasure (Pi\<^sub>M A (\<lambda>_. lborel))
2.37 +             (S A t \<inter> space (Pi\<^sub>M A (\<lambda>_. lborel))) = t ^ card A / fact (card A)"
2.38 +  using assms(1,2)
2.39 +proof (induction arbitrary: t rule: finite_induct)
2.40 +  case (empty t)
2.41 +  thus ?case by (simp add: PiM_empty S_def)
2.42 +next
2.43 +  case (insert b A t)
2.44 +  define n where "n = Suc (card A)"
2.45 +  have n_pos: "n > 0" by (simp add: n_def)
2.46 +  let ?M = "\<lambda>A. (Pi\<^sub>M A (\<lambda>_. lborel))"
2.47 +  {
2.48 +    fix A :: "'a set" and t :: real assume "finite A"
2.49 +    have "S A t \<inter> space (Pi\<^sub>M A (\<lambda>_. lborel)) =
2.50 +            Pi\<^sub>E A (\<lambda>_. {0..}) \<inter> (\<lambda>x. sum x A) -` {..t} \<inter> space (Pi\<^sub>M A (\<lambda>_. lborel))"
2.51 +      by (auto simp: S_def space_PiM)
2.52 +    also have "\<dots> \<in> sets (Pi\<^sub>M A (\<lambda>_. lborel))"
2.53 +      using \<open>finite A\<close> by measurable
2.54 +    finally have "S A t \<inter> space (Pi\<^sub>M A (\<lambda>_. lborel)) \<in> sets (Pi\<^sub>M A (\<lambda>_. lborel))" .
2.55 +  } note meas [measurable] = this
2.56 +
2.57 +  interpret product_sigma_finite "\<lambda>_. lborel"
2.58 +    by standard
2.59 +  have "emeasure (?M (insert b A)) (S (insert b A) t \<inter> space (?M (insert b A))) =
2.60 +        nn_integral (?M (insert b A))
2.61 +          (\<lambda>x. indicator (S (insert b A) t \<inter> space (?M (insert b A))) x)"
2.62 +    using insert.hyps by (subst nn_integral_indicator) auto
2.63 +  also have "\<dots> = (\<integral>\<^sup>+ y. \<integral>\<^sup>+ x. indicator (S (insert b A) t \<inter> space (?M (insert b A)))
2.64 +                    (x(b := y)) \<partial>?M A \<partial>lborel)"
2.65 +    using insert.prems insert.hyps by (intro product_nn_integral_insert_rev) auto
2.66 +  also have "\<dots> = (\<integral>\<^sup>+ y. \<integral>\<^sup>+ x. indicator {0..t} y * indicator (S A (t - y) \<inter> space (?M A)) x
2.67 +                    \<partial>?M A \<partial>lborel)"
2.68 +    using insert.hyps insert.prems emeasure_std_simplex_aux_step[of b A]
2.69 +    by (intro nn_integral_cong)
2.70 +       (auto simp: fun_eq_iff indicator_def space_PiM PiE_def extensional_def)
2.71 +  also have "\<dots> = (\<integral>\<^sup>+ y. indicator {0..t} y * (\<integral>\<^sup>+ x. indicator (S A (t - y) \<inter> space (?M A)) x
2.72 +                    \<partial>?M A) \<partial>lborel)" using \<open>finite A\<close>
2.73 +    by (subst nn_integral_cmult) auto
2.74 +  also have "\<dots> = (\<integral>\<^sup>+ y. indicator {0..t} y * emeasure (?M A) (S A (t - y) \<inter> space (?M A)) \<partial>lborel)"
2.75 +    using \<open>finite A\<close> by (subst nn_integral_indicator) auto
2.76 +  also have "\<dots> = (\<integral>\<^sup>+ y. indicator {0..t} y * (t - y) ^ card A / ennreal (fact (card A)) \<partial>lborel)"
2.77 +    using insert.IH by (intro nn_integral_cong) (auto simp: indicator_def divide_ennreal)
2.78 +  also have "\<dots> = (\<integral>\<^sup>+ y. indicator {0..t} y * (t - y) ^ card A \<partial>lborel) / ennreal (fact (card A))"
2.79 +    using \<open>finite A\<close> by (subst nn_integral_divide) auto
2.80 +  also have "(\<integral>\<^sup>+ y. indicator {0..t} y * (t - y) ^ card A \<partial>lborel) =
2.81 +               (\<integral>\<^sup>+y\<in>{0..t}. ennreal ((t - y) ^ (n - 1)) \<partial>lborel)"
2.82 +    by (intro nn_integral_cong) (auto simp: indicator_def n_def)
2.83 +  also have "((\<lambda>x. - ((t - x) ^ n / n)) has_real_derivative (t - x) ^ (n - 1)) (at x)"
2.84 +    if "x \<in> {0..t}" for x by (rule derivative_eq_intros refl | simp add: n_pos)+
2.85 +  hence "(\<integral>\<^sup>+y\<in>{0..t}. ennreal ((t - y) ^ (n - 1)) \<partial>lborel) =
2.86 +           ennreal (-((t - t) ^ n / n) - (-((t - 0) ^ n / n)))"
2.87 +    using insert.prems insert.hyps by (intro nn_integral_FTC_Icc) auto
2.88 +  also have "\<dots> = ennreal (t ^ n / n)" using n_pos by (simp add: zero_power)
2.89 +  also have "\<dots> / ennreal (fact (card A)) = ennreal (t ^ n / n / fact (card A))"
2.90 +    using n_pos \<open>t \<ge> 0\<close> by (subst divide_ennreal) auto
2.91 +  also have "t ^ n / n / fact (card A) = t ^ n / fact n"
2.92 +    by (simp add: n_def)
2.93 +  also have "n = card (insert b A)"
2.94 +    using insert.hyps by (subst card_insert) (auto simp: n_def)
2.95 +  finally show ?case .
2.96 +qed
2.97 +
2.98 +end
2.99 +
2.100 +lemma emeasure_std_simplex:
2.101 +  "emeasure lborel (convex hull (insert 0 Basis :: 'a :: euclidean_space set)) =
2.102 +     ennreal (1 / fact DIM('a))"
2.103 +proof -
2.104 +  have "emeasure lborel {x::'a. (\<forall>i\<in>Basis. 0 \<le> x \<bullet> i) \<and> sum ((\<bullet>) x) Basis \<le> 1} =
2.105 +               emeasure (distr (Pi\<^sub>M Basis (\<lambda>b. lborel)) borel (\<lambda>f. \<Sum>b\<in>Basis. f b *\<^sub>R b))
2.106 +                 {x::'a. (\<forall>i\<in>Basis. 0 \<le> x \<bullet> i) \<and> sum ((\<bullet>) x) Basis \<le> 1}"
2.107 +    by (subst lborel_eq) simp
2.108 +  also have "\<dots> = emeasure (Pi\<^sub>M Basis (\<lambda>b. lborel))
2.109 +                    ({y::'a \<Rightarrow> real. (\<forall>i\<in>Basis. 0 \<le> y i) \<and> sum y Basis \<le> 1} \<inter>
2.110 +                      space (Pi\<^sub>M Basis (\<lambda>b. lborel)))"
2.111 +    by (subst emeasure_distr) auto
2.112 +  also have "\<dots> = ennreal (1 / fact DIM('a))"
2.113 +    by (subst emeasure_std_simplex_aux) auto
2.114 +  finally show ?thesis by (simp only: std_simplex)
2.115 +qed
2.116 +
2.117 +theorem content_std_simplex:
2.118 +  "measure lborel (convex hull (insert 0 Basis :: 'a :: euclidean_space set)) =
2.119 +     1 / fact DIM('a)"
2.120 +  by (simp add: measure_def emeasure_std_simplex)
2.121 +
2.122 +(* TODO: move to Change_of_Vars *)
2.123 +lemma measure_lebesgue_linear_transformation:
2.124 +  fixes A :: "(real ^ 'n :: {finite, wellorder}) set"
2.125 +  fixes f :: "_ \<Rightarrow> real ^ 'n :: {finite, wellorder}"
2.126 +  assumes "bounded A" "A \<in> sets lebesgue" "linear f"
2.127 +  shows   "measure lebesgue (f ` A) = \<bar>det (matrix f)\<bar> * measure lebesgue A"
2.128 +proof -
2.129 +  from assms have [intro]: "A \<in> lmeasurable"
2.130 +    by (intro bounded_set_imp_lmeasurable) auto
2.131 +  hence [intro]: "f ` A \<in> lmeasurable"
2.132 +    by (intro lmeasure_integral measurable_linear_image assms)
2.133 +  have "measure lebesgue (f ` A) = integral (f ` A) (\<lambda>_. 1)"
2.134 +    by (intro lmeasure_integral measurable_linear_image assms) auto
2.135 +  also have "\<dots> = integral (f ` A) (\<lambda>_. 1 :: real ^ 1) \$ 0"
2.136 +    by (subst integral_component_eq_cart [symmetric]) (auto intro: integrable_on_const)
2.137 +  also have "\<dots> = \<bar>det (matrix f)\<bar> * integral A (\<lambda>x. 1 :: real ^ 1) \$ 0"
2.138 +    using assms
2.139 +    by (subst integral_change_of_variables_linear)
2.140 +       (auto simp: o_def absolutely_integrable_on_def intro: integrable_on_const)
2.141 +  also have "integral A (\<lambda>x. 1 :: real ^ 1) \$ 0 = integral A (\<lambda>x. 1)"
2.142 +    by (subst integral_component_eq_cart [symmetric]) (auto intro: integrable_on_const)
2.143 +  also have "\<dots> = measure lebesgue A"
2.144 +    by (intro lmeasure_integral [symmetric]) auto
2.145 +  finally show ?thesis .
2.146 +qed
2.147 +
2.148 +theorem content_simplex:
2.149 +  fixes X :: "(real ^ 'n :: {finite, wellorder}) set" and f :: "'n :: _ \<Rightarrow> real ^ ('n :: _)"
2.150 +  assumes "finite X" "card X = Suc CARD('n)" and x0: "x0 \<in> X" and bij: "bij_betw f UNIV (X - {x0})"
2.151 +  defines "M \<equiv> (\<chi> i. \<chi> j. f j \$ i - x0 \$ i)"
2.152 +  shows "content (convex hull X) = \<bar>det M\<bar> / fact (CARD('n))"
2.153 +proof -
2.154 +  define g where "g = (\<lambda>x. M *v x)"
2.155 +  have [simp]: "M *v axis i 1 = f i - x0" for i :: 'n
2.156 +    by (simp add: M_def matrix_vector_mult_basis column_def vec_eq_iff)
2.157 +  define std where "std = (convex hull insert 0 Basis :: (real ^ 'n :: _) set)"
2.158 +  have compact: "compact std" unfolding std_def
2.159 +    by (intro finite_imp_compact_convex_hull) auto
2.160 +
2.161 +  have "measure lebesgue (convex hull X) = measure lebesgue (((+) (-x0)) ` (convex hull X))"
2.162 +    by (rule measure_translation [symmetric])
2.163 +  also have "((+) (-x0)) ` (convex hull X) = convex hull (((+) (-x0)) ` X)"
2.164 +    by (rule convex_hull_translation [symmetric])
2.165 +  also have "((+) (-x0)) ` X = insert 0 ((\<lambda>x. x - x0) ` (X - {x0}))"
2.166 +    using x0 by (auto simp: image_iff)
2.167 +  finally have eq: "measure lebesgue (convex hull X) = measure lebesgue (convex hull \<dots>)" .
2.168 +
2.169 +  from compact have "measure lebesgue (g ` std) = \<bar>det M\<bar> * measure lebesgue std"
2.170 +    by (subst measure_lebesgue_linear_transformation)
2.171 +       (auto intro: finite_imp_bounded_convex_hull dest: compact_imp_closed simp: g_def std_def)
2.172 +  also have "measure lebesgue std = content std" using compact
2.173 +    by (intro measure_completion) (auto dest: compact_imp_closed)
2.174 +  also have "content std = 1 / fact CARD('n)" unfolding std_def
2.175 +    by (simp add: content_std_simplex)
2.176 +  also have "g ` std = convex hull (g ` insert 0 Basis)" unfolding std_def
2.177 +    by (rule convex_hull_linear_image) (auto simp: g_def)
2.178 +  also have "g ` insert 0 Basis = insert 0 (g ` Basis)"
2.179 +    by (auto simp: g_def)
2.180 +  also have "g ` Basis = (\<lambda>x. x - x0) ` range f"
2.181 +    by (auto simp: g_def Basis_vec_def image_iff)
2.182 +  also have "range f = X - {x0}" using bij
2.183 +    using bij_betw_imp_surj_on by blast
2.184 +  also note eq [symmetric]
2.185 +  finally show ?thesis
2.186 +    using finite_imp_compact_convex_hull[OF \<open>finite X\<close>] by (auto dest: compact_imp_closed)
2.187 +qed
2.188 +
2.189 +theorem content_triangle:
2.190 +  fixes A B C :: "real ^ 2"
2.191 +  shows "content (convex hull {A, B, C}) =
2.192 +           \<bar>(C \$ 1 - A \$ 1) * (B \$ 2 - A \$ 2) - (B \$ 1 - A \$ 1) * (C \$ 2 - A \$ 2)\<bar> / 2"
2.193 +proof -
2.194 +  define M :: "real ^ 2 ^ 2" where "M \<equiv> (\<chi> i. \<chi> j. (if j = 1 then B else C) \$ i - A \$ i)"
2.195 +  define g where "g = (\<lambda>x. M *v x)"
2.196 +  define std where "std = (convex hull insert 0 Basis :: (real ^ 2) set)"
2.197 +  have [simp]: "M *v axis i 1 = (if i = 1 then B - A else C - A)" for i
2.198 +    by (auto simp: M_def matrix_vector_mult_basis column_def vec_eq_iff)
2.199 +  have compact: "compact std" unfolding std_def
2.200 +    by (intro finite_imp_compact_convex_hull) auto
2.201 +
2.202 +  have "measure lebesgue (convex hull {A, B, C}) =
2.203 +          measure lebesgue (((+) (-A)) ` (convex hull {A, B, C}))"
2.204 +    by (rule measure_translation [symmetric])
2.205 +  also have "((+) (-A)) ` (convex hull {A, B, C}) = convex hull (((+) (-A)) ` {A, B, C})"
2.206 +    by (rule convex_hull_translation [symmetric])
2.207 +  also have "((+) (-A)) ` {A, B, C} = {0, B - A, C - A}"
2.208 +    by (auto simp: image_iff)
2.209 +  finally have eq: "measure lebesgue (convex hull {A, B, C}) =
2.210 +                      measure lebesgue (convex hull {0, B - A, C - A})" .
2.211 +
2.212 +  from compact have "measure lebesgue (g ` std) = \<bar>det M\<bar> * measure lebesgue std"
2.213 +    by (subst measure_lebesgue_linear_transformation)
2.214 +       (auto intro: finite_imp_bounded_convex_hull dest: compact_imp_closed simp: g_def std_def)
2.215 +  also have "measure lebesgue std = content std" using compact
2.216 +    by (intro measure_completion) (auto dest: compact_imp_closed)
2.217 +  also have "content std = 1 / 2" unfolding std_def
2.218 +    by (simp add: content_std_simplex)
2.219 +  also have "g ` std = convex hull (g ` insert 0 Basis)" unfolding std_def
2.220 +    by (rule convex_hull_linear_image) (auto simp: g_def)
2.221 +  also have "g ` insert 0 Basis = insert 0 (g ` Basis)"
2.222 +    by (auto simp: g_def)
2.223 +  also have "(2 :: 2) \<noteq> 1" by auto
2.224 +  hence "\<not>(\<forall>y::2. y = 1)" by blast
2.225 +  hence "g ` Basis = {B - A, C - A}"
2.226 +    by (auto simp: g_def Basis_vec_def image_iff)
2.227 +  also note eq [symmetric]
2.228 +  finally show ?thesis
2.229 +    using finite_imp_compact_convex_hull[of "{A, B, C}"]
2.230 +    by (auto dest!: compact_imp_closed simp: det_2 M_def)
2.231 +qed
2.232 +
2.233 +theorem heron:
2.234 +  fixes A B C :: "real ^ 2"
2.235 +  defines "a \<equiv> dist B C" and "b \<equiv> dist A C" and "c \<equiv> dist A B"
2.236 +  defines "s \<equiv> (a + b + c) / 2"
2.237 +  shows   "content (convex hull {A, B, C}) = sqrt (s * (s - a) * (s - b) * (s - c))"
2.238 +proof -
2.239 +  have [simp]: "(UNIV :: 2 set) = {1, 2}"
2.240 +    using exhaust_2 by auto
2.241 +  have dist_eq: "dist (A :: real ^ 2) B ^ 2 = (A \$ 1 - B \$ 1) ^ 2 + (A \$ 2 - B \$ 2) ^ 2"
2.242 +    for A B by (simp add: dist_vec_def dist_real_def)
2.243 +  have nonneg: "s * (s - a) * (s - b) * (s - c) \<ge> 0"
2.244 +    using dist_triangle[of A B C] dist_triangle[of A C B] dist_triangle[of B C A]
2.245 +    by (intro mult_nonneg_nonneg) (auto simp: s_def a_def b_def c_def dist_commute)
2.246 +
2.247 +  have "16 * content (convex hull {A, B, C}) ^ 2 =
2.248 +          4 * ((C \$ 1 - A \$ 1) * (B \$ 2 - A \$ 2) - (B \$ 1 - A \$ 1) * (C \$ 2 - A \$ 2)) ^ 2"
2.249 +    by (subst content_triangle) (simp add: power_divide)
2.250 +  also have "\<dots> = (2 * (dist A B ^ 2 * dist A C ^ 2 + dist A B ^ 2 * dist B C ^ 2 +
2.251 +      dist A C ^ 2 * dist B C ^ 2) - (dist A B ^ 2) ^ 2 - (dist A C ^ 2) ^ 2 - (dist B C ^ 2) ^ 2)"
2.252 +    unfolding dist_eq unfolding power2_eq_square by algebra
2.253 +  also have "\<dots> = (a + b + c) * ((a + b + c) - 2 * a) * ((a + b + c) - 2 * b) *
2.254 +                    ((a + b + c) - 2 * c)"
2.255 +    unfolding power2_eq_square by (simp add: s_def a_def b_def c_def algebra_simps)
2.256 +  also have "\<dots> = 16 * s * (s - a) * (s - b) * (s - c)"
2.257 +    by (simp add: s_def divide_simps mult_ac)
2.258 +  finally have "content (convex hull {A, B, C}) ^ 2 = s * (s - a) * (s - b) * (s - c)"
2.259 +    by simp
2.260 +  also have "\<dots> = sqrt (s * (s - a) * (s - b) * (s - c)) ^ 2"
2.261 +    by (intro real_sqrt_pow2 [symmetric] nonneg)
2.262 +  finally show ?thesis using nonneg
2.263 +    by (subst (asm) power2_eq_iff_nonneg) auto
2.264 +qed
2.265 +
2.266 +end
```