Quelle  Ball_Volume.thy

(*
   File:     HOL/Analysis/Ball_Volume.thy
   Author:   Manuel Eberl, TU München
*)

section ‹The Volume of an ‹n›-Dimensional Ball›

theory Ball_Volume
  imports Gamma_Function Lebesgue_Integral_Substitution
begin

text ‹
 We define the volume of the unit ball in terms of the Gamma function. Note that the
 dimension need not be an integer; we also allow fractional dimensions, although we do
 not use this case or prove anything about it for now.
 ›
definition✐‹tag important› unit_ball_vol :: "real → real" where
  "unit_ball_vol n = pi powr (n / 2) / Gamma (n / 2 + 1)"

lemma unit_ball_vol_pos [simp]: "n ≥ 0 ==> unit_ball_vol n > 0"
  by (force simp: unit_ball_vol_def intro: divide_nonneg_pos)

lemma unit_ball_vol_nonneg [simp]: "n ≥ 0 ==> unit_ball_vol n ≥ 0"
  by (simp add: dual_order.strict_implies_order)

text ‹
 We first need the value of the following integral, which is at the core of
 computing the measure of an ‹n + 1›-dimensional ball in terms of the measure of an
 ‹n›-dimensional one.
 ›
lemma emeasure_cball_aux_integral:
  "(∫⁺x. indicator {-1..1} x * sqrt (1 - x²) ^ n ∂lborel) =
      ennreal (Beta (1 / 2) (real n / 2 + 1))"
proof -
  have "((λt. t powr (-1 / 2) * (1 - t) powr (real n / 2)) has_integral
          Beta (1 / 2) (real n / 2 + 1)) {0..1}"
    using has_integral_Beta_real[of "1/2" "n / 2 + 1"] by simp
  from nn_integral_has_integral_lebesgue[OF _ this] have
     "ennreal (Beta (1 / 2) (real n / 2 + 1)) =
        nn_integral lborel (λt. ennreal (t powr (-1 / 2) * (1 - t) powr (real n / 2) *
                                indicator {0^2..1^2} t))"
    by (simp add: mult_ac ennreal_mult' ennreal_indicator)
  also have "… = (∫⁺ x. ennreal (x² powr - (1 / 2) * (1 - x²) powr (real n / 2) * (2 * x) *
                          indicator {0..1} x) ∂lborel)"
    by (subst nn_integral_substitution[where g = "λx. x ^ 2" and g' = "λx. 2 * x"])
       (auto intro!: derivative_eq_intros continuous_intros simp: set_borel_measurable_def)
  also have "… = (∫⁺ x. 2 * ennreal ((1 - x²) powr (real n / 2) * indicator {0..1} x) ∂lborel)"
    by (intro nn_integral_cong_AE AE_I[of _ _ "{0}"])
       (auto simp: indicator_def powr_minus powr_half_sqrt field_split_simps ennreal_mult')
  also have "… = (∫⁺ x. ennreal ((1 - x²) powr (real n / 2) * indicator {0..1} x) ∂lborel) +
                    (∫⁺ x. ennreal ((1 - x²) powr (real n / 2) * indicator {0..1} x) ∂lborel)"
    (is "_ = ?I + _") by (simp add: mult_2 nn_integral_add)
  also have "?I = (∫⁺ x. ennreal ((1 - x²) powr (real n / 2) * indicator {-1..0} x) ∂lborel)"
    by (subst nn_integral_real_affine[of _ "-1" 0])
       (auto simp: indicator_def intro!: nn_integral_cong)
  hence "?I + ?I = … + ?I" by simp
  also have "… = (∫⁺ x. ennreal ((1 - x²) powr (real n / 2) *
                    (indicator {-1..0} x + indicator{0..1} x)) ∂lborel)"
    by (subst nn_integral_add [symmetric]) (auto simp: algebra_simps)
  also have "… = (∫⁺ x. ennreal ((1 - x²) powr (real n / 2) * indicator {-1..1} x) ∂lborel)"
    by (intro nn_integral_cong_AE AE_I[of _ _ "{0}"]) (auto simp: indicator_def)
  also have "… = (∫⁺ x. ennreal (indicator {-1..1} x * sqrt (1 - x²) ^ n) ∂lborel)"
    by (intro nn_integral_cong_AE AE_I[of _ _ "{1, -1}"])
       (auto simp: powr_half_sqrt [symmetric] indicator_def abs_square_le_1
          abs_square_eq_1 powr_def exp_of_nat_mult [symmetric] emeasure_lborel_countable)
  finally show ?thesis ..
qed

lemma real_sqrt_le_iff': "x ≥ 0 ==> y ≥ 0 ==> sqrt x ≤ y ⟷ x ≤ y ^ 2"
  using real_le_lsqrt sqrt_le_D by blast

text ‹
 Isabelle's type system makes it very difficult to do an induction over the dimension
 of a Euclidean space type, because the type would change in the inductive step. To avoid
 this problem, we instead formulate the problem in a more concrete way by unfolding the
 definition of the Euclidean norm.
 ›
lemma emeasure_cball_aux:
  assumes "finite A" "r > 0"
  shows   "emeasure (Pi_M A (λ_. lborel))
             ({f. sqrt (∑i∈A. (f i)²) ≤ r} ∩ space (Pi_M A (λ_. lborel))) =
             ennreal (unit_ball_vol (real (card A)) * r ^ card A)"
  using assms
proof (induction arbitrary: r)
  case (empty r)
  thus ?case
    by (simp add: unit_ball_vol_def space_PiM)
next
  case (insert i A r)
  interpret product_sigma_finite "λ_. lborel"
    by standard
  have "emeasure (Pi_M (insert i A) (λ_. lborel))
            ({f. sqrt (∑i∈insert i A. (f i)²) ≤ r} ∩ space (Pi_M (insert i A) (λ_. lborel))) =
        nn_integral (Pi_M (insert i A) (λ_. lborel))
          (indicator ({f. sqrt (∑i∈insert i A. (f i)²) ≤ r} ∩
          space (Pi_M (insert i A) (λ_. lborel))))"
    by (subst nn_integral_indicator) auto
  also have "… = (∫⁺ y. ∫⁺ x. indicator ({f. sqrt ((f i)² + (∑i∈A. (f i)²)) ≤ r} ∩
                                space (Pi_M (insert i A) (λ_. lborel))) (x(i := y))
                   ∂Pi_M A (λ_. lborel) ∂lborel)"
    using insert.prems insert.hyps by (subst product_nn_integral_insert_rev) auto
  also have "… = (∫⁺ (y::real). ∫⁺ x. indicator {-r..r} y * indicator ({f. sqrt ((∑i∈A. (f i)²)) ≤
               sqrt (r ^ 2 - y ^ 2)} ∩ space (Pi_M A (λ_. lborel))) x ∂Pi_M A (λ_. lborel) ∂lborel)"
  proof (intro nn_integral_cong, goal_cases)
    case (1 y f)
    have *: "y ∈ {-r..r}" if "y ^ 2 + c ≤ r ^ 2" "c ≥ 0" for c
    proof -
      have "y ^ 2 ≤ y ^ 2 + c" using that by simp
      also have "… ≤ r ^ 2" by fact
      finally show ?thesis
        using ‹r > 0› by (simp add: power2_le_iff_abs_le abs_if split: if_splits)
    qed
    have "(∑x∈A. (if x = i then y else f x)²) = (∑x∈A. (f x)²)"
      using insert.hyps by (intro sum.cong) auto
    thus ?case using 1 ‹r > 0›
      by (auto simp: sum_nonneg real_sqrt_le_iff' indicator_def PiE_def space_PiM dest!: *)
  qed
  also have "… = (∫⁺ (y::real). indicator {-r..r} y * (∫⁺ x. indicator ({f. sqrt ((∑i∈A. (f i)²))
                                   ≤ sqrt (r ^ 2 - y ^ 2)} ∩ space (Pi_M A (λ_. lborel))) x
                  ∂Pi_M A (λ_. lborel)) ∂lborel)" by (subst nn_integral_cmult) auto
  also have "… = (∫⁺ (y::real). indicator {-r..r} y * emeasure (PiM A (λ_. lborel))
      ({f. sqrt ((∑i∈A. (f i)²)) ≤ sqrt (r ^ 2 - y ^ 2)} ∩ space (Pi_M A (λ_. lborel))) ∂lborel)"
    using ‹finite A› by (intro nn_integral_cong, subst nn_integral_indicator) auto
  also have "… = (∫⁺ (y::real). indicator {-r..r} y * ennreal (unit_ball_vol (real (card A)) *
                                  (sqrt (r ^ 2 - y ^ 2)) ^ card A) ∂lborel)"
  proof (intro nn_integral_cong_AE, goal_cases)
    case 1
    have "AE y in lborel. y ∉ {-r,r}"
      by (intro AE_not_in countable_imp_null_set_lborel) auto
    thus ?case
    proof eventually_elim
      case (elim y)
      show ?case
      proof (cases "y ∈ {-r<..<r}")
        case True
        hence "y² < r²" by (subst real_sqrt_less_iff [symmetric]) auto
        thus ?thesis by (subst insert.IH) (auto)
      qed (insert elim, auto)
    qed
  qed
  also have "… = ennreal (unit_ball_vol (real (card A))) *
                    (∫⁺ (y::real). indicator {-r..r} y * (sqrt (r ^ 2 - y ^ 2)) ^ card A ∂lborel)"
    by (subst nn_integral_cmult [symmetric])
       (auto simp: mult_ac ennreal_mult' [symmetric] indicator_def intro!: nn_integral_cong)
  also have "(∫⁺ (y::real). indicator {-r..r} y * (sqrt (r ^ 2 - y ^ 2)) ^ card A ∂lborel) =
               (∫⁺ (y::real). r ^ card A * indicator {-1..1} y * (sqrt (1 - y ^ 2)) ^ card A
               ∂(distr lborel borel ((*) (1/r))))" using ‹r > 0›
    by (subst nn_integral_distr)
       (auto simp: indicator_def field_simps real_sqrt_divide intro!: nn_integral_cong)
  also have "… = (∫⁺ x. ennreal (r ^ Suc (card A)) *
               (indicator {- 1..1} x * sqrt (1 - x²) ^ card A) ∂lborel)" using ‹r > 0›
    by (subst lborel_distr_mult) (auto simp: nn_integral_density ennreal_mult' [symmetric] mult_ac)
  also have "… = ennreal (r ^ Suc (card A)) * (∫⁺ x. indicator {- 1..1} x *
                    sqrt (1 - x²) ^ card A ∂lborel)"
    by (subst nn_integral_cmult) auto
  also note emeasure_cball_aux_integral
  also have "ennreal (unit_ball_vol (real (card A))) * (ennreal (r ^ Suc (card A)) *
                 ennreal (Beta (1/2) (card A / 2 + 1))) =
               ennreal (unit_ball_vol (card A) * Beta (1/2) (card A / 2 + 1) * r ^ Suc (card A))"
    using ‹r > 0› by (simp add: ennreal_mult' [symmetric] mult_ac)
  also have "unit_ball_vol (card A) * Beta (1/2) (card A / 2 + 1) = unit_ball_vol (Suc (card A))"
    by (auto simp: unit_ball_vol_def Beta_def Gamma_eq_zero_iff field_simps
          Gamma_one_half_real powr_half_sqrt [symmetric] powr_add [symmetric])
  also have "Suc (card A) = card (insert i A)" using insert.hyps by simp
  finally show ?case .
qed


text ‹
 We now get the main theorem very easily by just applying the above lemma.
 ›
context
  fixes c :: "'a :: euclidean_space" and r :: real
  assumes r: "r ≥ 0"
begin

theorem✐‹tag unimportant› emeasure_cball:
  "emeasure lborel (cball c r) = ennreal (unit_ball_vol (DIM('a)) * r ^ DIM('a))"
proof (cases "r = 0")
  case False
  with r have r: "r > 0" by simp
  have "(lborel :: 'a measure) =
          distr (Pi_M Basis (λ_. lborel)) borel (λf. ∑b∈Basis. f b *_R b)"
    by (rule lborel_eq)
  also have "emeasure … (cball 0 r) =
               emeasure (Pi_M Basis (λ_. lborel))
               ({y. dist 0 (∑b∈Basis. y b *_R b :: 'a) ≤ r} ∩ space (Pi_M Basis (λ_. lborel)))"
    by (subst emeasure_distr) (auto simp: cball_def)
  also have "{f. dist 0 (∑b∈Basis. f b *_R b :: 'a) ≤ r} = {f. sqrt (∑i∈Basis. (f i)²) ≤ r}"
    by (subst euclidean_dist_l2) (auto simp: L2_set_def)
  also have "emeasure (Pi_M Basis (λ_. lborel)) (… ∩ space (Pi_M Basis (λ_. lborel))) =
               ennreal (unit_ball_vol (real DIM('a)) * r ^ DIM('a))"
    using r by (subst emeasure_cball_aux) simp_all
  also have "emeasure lborel (cball 0 r :: 'a set) =
               emeasure (distr lborel borel (λx. c + x)) (cball c r)"
    by (subst emeasure_distr) (auto simp: cball_def dist_norm norm_minus_commute)
  also have "distr lborel borel (λx. c + x) = lborel"
    using lborel_affine[of 1 c] by (simp add: density_1)
  finally show ?thesis .
qed auto

corollary✐‹tag unimportant› content_cball:
  "content (cball c r) = unit_ball_vol (DIM('a)) * r ^ DIM('a)"
  by (simp add: measure_def emeasure_cball r)

corollary✐‹tag unimportant› emeasure_ball:
  "emeasure lborel (ball c r) = ennreal (unit_ball_vol (DIM('a)) * r ^ DIM('a))"
proof -
  from negligible_sphere[of c r] have "sphere c r ∈ null_sets lborel"
    by (auto simp: null_sets_completion_iff negligible_iff_null_sets negligible_convex_frontier)
  hence "emeasure lborel (ball c r ∪ sphere c r :: 'a set) = emeasure lborel (ball c r :: 'a set)"
    by (intro emeasure_Un_null_set) auto
  also have "ball c r ∪ sphere c r = (cball c r :: 'a set)" by auto
  also have "emeasure lborel … = ennreal (unit_ball_vol (real DIM('a)) * r ^ DIM('a))"
    by (rule emeasure_cball)
  finally show ?thesis ..
qed

corollary✐‹tag important› content_ball:
  "content (ball c r) = unit_ball_vol (DIM('a)) * r ^ DIM('a)"
  by (simp add: measure_def r emeasure_ball)

end


text ‹
 Lastly, we now prove some nicer explicit formulas for the volume of the unit balls in
 the cases of even and odd integer dimensions.
 ›
lemma unit_ball_vol_even:
  "unit_ball_vol (real (2 * n)) = pi ^ n / fact n"
  by (simp add: unit_ball_vol_def add_ac powr_realpow Gamma_fact)

lemma unit_ball_vol_odd':
        "unit_ball_vol (real (2 * n + 1)) = pi ^ n / pochhammer (1 / 2) (Suc n)"
  and unit_ball_vol_odd:
        "unit_ball_vol (real (2 * n + 1)) =
           (2 ^ (2 * Suc n) * fact (Suc n)) / fact (2 * Suc n) * pi ^ n"
proof -
  have "unit_ball_vol (real (2 * n + 1)) =
          pi powr (real n + 1 / 2) / Gamma (1 / 2 + real (Suc n))"
    by (simp add: unit_ball_vol_def field_simps)
  also have "pochhammer (1 / 2) (Suc n) = Gamma (1 / 2 + real (Suc n)) / Gamma (1 / 2)"
    by (intro pochhammer_Gamma) auto
  hence "Gamma (1 / 2 + real (Suc n)) = sqrt pi * pochhammer (1 / 2) (Suc n)"
    by (simp add: Gamma_one_half_real)
  also have "pi powr (real n + 1 / 2) / … = pi ^ n / pochhammer (1 / 2) (Suc n)"
    by (simp add: powr_add powr_half_sqrt powr_realpow)
  finally show "unit_ball_vol (real (2 * n + 1)) = …" .
  also have "pochhammer (1 / 2 :: real) (Suc n) =
               fact (2 * Suc n) / (2 ^ (2 * Suc n) * fact (Suc n))"
    using fact_double[of "Suc n", where ?'a = real] by (simp add: divide_simps mult_ac)
  also have "pi ^n / … = (2 ^ (2 * Suc n) * fact (Suc n)) / fact (2 * Suc n) * pi ^ n"
    by simp
  finally show "unit_ball_vol (real (2 * n + 1)) = …" .
qed

lemma unit_ball_vol_numeral:
  "unit_ball_vol (numeral (Num.Bit0 n)) = pi ^ numeral n / fact (numeral n)" (is ?th1)
  "unit_ball_vol (numeral (Num.Bit1 n)) = 2 ^ (2 * Suc (numeral n)) * fact (Suc (numeral n)) /
    fact (2 * Suc (numeral n)) * pi ^ numeral n" (is ?th2)
proof -
  have "numeral (Num.Bit0 n) = (2 * numeral n :: nat)"
    by (simp only: numeral_Bit0 mult_2 ring_distribs)
  also have "unit_ball_vol … = pi ^ numeral n / fact (numeral n)"
    by (rule unit_ball_vol_even)
  finally show ?th1 by simp
next
  have "numeral (Num.Bit1 n) = (2 * numeral n + 1 :: nat)"
    by (simp only: numeral_Bit1 mult_2)
  also have "unit_ball_vol … = 2 ^ (2 * Suc (numeral n)) * fact (Suc (numeral n)) /
                                  fact (2 * Suc (numeral n)) * pi ^ numeral n"
    by (rule unit_ball_vol_odd)
  finally show ?th2 by simp
qed

lemmas eval_unit_ball_vol = unit_ball_vol_numeral fact_numeral


text ‹
 Just for fun, we compute the volume of unit balls for a few dimensions.
 ›
lemma unit_ball_vol_0 [simp]: "unit_ball_vol 0 = 1"
  using unit_ball_vol_even[of 0] by simp

lemma unit_ball_vol_1 [simp]: "unit_ball_vol 1 = 2"
  using unit_ball_vol_odd[of 0] by simp

corollary✐‹tag unimportant›
          unit_ball_vol_2: "unit_ball_vol 2 = pi"
      and unit_ball_vol_3: "unit_ball_vol 3 = 4 / 3 * pi"
      and unit_ball_vol_4: "unit_ball_vol 4 = pi² / 2"
      and unit_ball_vol_5: "unit_ball_vol 5 = 8 / 15 * pi²"
  by (simp_all add: eval_unit_ball_vol)

corollary✐‹tag unimportant› circle_area:
  "r ≥ 0 ==> content (ball c r :: (real ^ 2) set) = r ^ 2 * pi"
  by (simp add: content_ball unit_ball_vol_2)

corollary✐‹tag unimportant› sphere_volume:
  "r ≥ 0 ==> content (ball c r :: (real ^ 3) set) = 4 / 3 * r ^ 3 * pi"
  by (simp add: content_ball unit_ball_vol_3)

text ‹
 Useful equivalent forms
 ›
corollary✐‹tag unimportant› content_ball_eq_0_iff [simp]: "content (ball c r) = 0 ⟷ r ≤ 0"
proof -
  have "r > 0 ==> content (ball c r) > 0"
    by (simp add: content_ball unit_ball_vol_def)
  then show ?thesis
    by (fastforce simp: ball_empty)
qed

corollary✐‹tag unimportant› content_ball_gt_0_iff [simp]: "0 < content (ball z r) ⟷ 0 < r"
  by (auto simp: zero_less_measure_iff)

corollary✐‹tag unimportant› content_cball_eq_0_iff [simp]: "content (cball c r) = 0 ⟷ r ≤ 0"
proof (cases "r = 0")
  case False
  moreover have "r > 0 ==> content (cball c r) > 0"
    by (simp add: content_cball unit_ball_vol_def)
  ultimately show ?thesis
    by fastforce
qed auto

corollary✐‹tag unimportant› content_cball_gt_0_iff [simp]: "0 < content (cball z r) ⟷ 0 < r"
  by (auto simp: zero_less_measure_iff)

end