Library FEM.geom_simplex

This file is part of the Coq Numerical Analysis library
Copyright (C) Boldo, Clément, Martin, Mayero, Mouhcine
This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the COPYING file for more details.

Bibliography

[RR9557v1] François Clément and Vincent Martin, Finite element method. Detailed proofs to be formalized in Coq, RR-9557, first version, Inria Paris, 2024, https://inria.hal.science/hal-04713897v1.

From Requisite Require Import stdlib_wR ssr_wMC.

From Coquelicot Require Import Hierarchy.
Close Scope R_scope.

From Algebra Require Import Algebra_wDep.
From FEM Require Import multi_index.

Local Open Scope R_scope.

Section R_compl.

Lemma INR_n0 : {n}, (0 < n)%coq_nat INR n 0.
Proof. move=>> /Nat.lt_neq H; apply not_0_INR; easy. Qed.

Lemma INR_pred : {n}, (0 < n)%coq_nat INR n.-1 = INR n - 1.
Proof. intros; rewrite -Nat.sub_1_r minus_INR; auto with arith. Qed.

Lemma INR_inv_eq :
   {n}, (0 < n)%coq_nat / INR n = (1 - INR n.-1 / INR n)%G.
Proof.
intros; rewrite INR_pred//; unfold minus, opp, plus; simpl.
rewrite Rdiv_plus_distr Rdiv_diag;
    [lra | apply not_0_INR, not_eq_sym, Nat.lt_neq; easy].
Qed.

Lemma INR_inv_pred :
   {n}, (0 < n)%coq_nat (/ INR n - 1%K)%G = - (INR n.-1 / INR n)%G.
Proof. intros; rewrite INR_inv_eq// minus2_l; easy. Qed.

Lemma INR_inv_decomp :
   {n}, (0 < n)%coq_nat (/ INR n + INR n.-1 / INR n)%M = 1.
Proof. intros; rewrite INR_inv_eq// plus_minus_l; easy. Qed.

Lemma INR_inv_n0 : {n}, n 0%N / INR n 0.
Proof. intros; apply Rinv_neq_0_compat, not_0_INR; easy. Qed.

End R_compl.

Section Vertices_Nodes_dk.

Context {d k : nat}.

[RR9557v1] Def 1434, p. 55.
See also Def 1365, p. 40.
Definition vtx_ref : 'R^{d.+1,d} :=
  fun imatch (ord_eq_dec i ord0) with
    | left _zero
    | right _kronR (i - 1)%coq_nat
    end.

Lemma vtx_ref_0 : i, i = ord0 vtx_ref i = 0%M.
Proof. intros; unfold vtx_ref; case (ord_eq_dec _ _); easy. Qed.

Lemma vtx_ref_S0 :
   i (j : 'I_d), i ord0 (i - 1)%coq_nat j vtx_ref i j = 0.
Proof.
intros i j Hi Hj; unfold vtx_ref;
    destruct (ord_eq_dec _ _); [| apply kron_is_0]; easy.
Qed.

Lemma vtx_ref_S1 :
   i (j : 'I_d), i ord0 (i - 1)%coq_nat = j vtx_ref i j = 1.
Proof.
intros i j Hi Hj; unfold vtx_ref;
    destruct (ord_eq_dec _ _); [| apply kron_is_1]; easy.
Qed.

Lemma vtx_ref_ge_0 : i j, 0 vtx_ref i j.
Proof.
intros i j; destruct (ord_eq_dec i ord0) as [Hi | Hi].
rewrite vtx_ref_0; easy.
destruct (Nat.eq_dec (i - 1)%coq_nat j) as [Hj | Hj].
rewrite vtx_ref_S1//; apply Rle_0_1.
rewrite vtx_ref_S0; easy.
Qed.

Lemma vtx_ref_kron :
   {i} (Hi : i ord0), vtx_ref i = kronR (lower_S Hi).
Proof.
intros i Hi; apply extF; intros j;
    destruct (Nat.eq_dec (i - 1)%coq_nat j) as [Hj | Hj].
rewrite vtx_ref_S1// kronR_is_1//.
rewrite kronR_is_0// vtx_ref_S0//; easy.
Qed.

Lemma vtx_ref_liftF_S : i, liftF_S vtx_ref i = kronR i.
Proof.
intros; apply extF; intro; unfold liftF_S.
rewrite (vtx_ref_kron (lift_S_not_first _)) lower_lift_S; easy.
Qed.

[RR9557v1] Lem 1436, p. 55.
Lemma vtx_ref_affine_independent : affine_independent vtx_ref.
Proof.
intros L HL; rewrite -(lc_kron_r L) -HL; apply lc_ext_r; intros i.
apply extF; intros j; rewrite 2!fct_minus_eq -(minus_zero_r (kron _ _ _)).
rewrite vtx_ref_liftF_S; f_equal.
rewrite constF_correct; rewrite vtx_ref_0; easy.
Qed.

[RR9557v1] Lem 1591, Eq (9.84), p. 98.
a_alpha = (1 - |alpha|/k) v_0 + \sum (alpha_i / k) v_i.
This definition uses total function Rinv. It will be used only for k>0, and there is a specific definition for k=0 down below.
Definition node_cur (vtx_cur : 'R^{d.+1,d}) : 'R^{(pbinom d k).+1,d} :=
  fun ipk
    (scal (1 - (sum (scal (/ INR k) (mapF INR (Adk d k ipk)))))%G
          (vtx_cur ord0) +
    lin_comb (scal (/ INR k) (mapF INR (Adk d k ipk))) (liftF_S vtx_cur))%M.

Definition node_ref : 'R^{(pbinom d k).+1,d} :=
  fun ipk iINR (Adk d k ipk i) / INR k.

[RR9557v1] Lem 1599, Eq (9.92), p. 100.
Lemma node_ref_cur : node_ref = node_cur vtx_ref.
Proof.
unfold node_ref, node_cur; apply extF; intros ipk; apply extF; intros i.
rewrite (lc_ext_r _ vtx_ref_liftF_S) lc_kron_r.
rewrite fct_plus_eq fct_scal_r_eq vtx_ref_0; [| easy].
rewrite scal_zero_r plus_zero_l fct_scal_eq mapF_correct.
rewrite -inv_eq_R -div_eq_R; unfold div; rewrite mult_comm_R; easy.
Qed.

Lemma node_ref_eq :
   ipk, node_ref ipk = scal (/ INR k) (mapF INR (Adk d k ipk)).
Proof.
intros; unfold node_ref; apply extF; intro.
unfold Rdiv; rewrite Rmult_comm; easy.
Qed.

[RR9557v1] Def 1588, Eq (9.81), p. 97.
a_alpha = v_0 + \sum (alpha_i / k) (v_i - v_0).
The choice for 0 is arbitrary, it could be any index k<=d.
Lemma node_cur_eq :
   (vtx_cur : 'R^{d.+1,d}),
    node_cur vtx_cur = fun ipk
      (vtx_cur ord0 + lin_comb (scal (/ INR k) (mapF INR (Adk d k ipk)))
        (liftF_S vtx_cur - constF d (vtx_cur ord0))%G)%M.
Proof.
intros vtx_cur; unfold node_cur; apply extF; intros ipk.
rewrite scal_minus_l scal_one lc_minus_r lc_constF_r.
rewrite minus_eq -!plus_assoc; apply plus_compat_l, plus_comm.
Qed.

Lemma node_cur_inj :
   (vtx_cur : 'R^{d.+1,d}),
    affine_independent vtx_cur injective (node_cur vtx_cur).
Proof.
case (Nat.eq_dec k 0); intros Hk vtx_cur Hvtx.
intros [i Hi] [j Hj] _; subst; apply ord_inj; simpl.
rewrite pbinom_0_r in Hi, Hj; move: Hi ⇒ /ltP Hi; move: Hj ⇒ /ltP Hj; lia.
intros i j; rewrite node_cur_eq.
move⇒ /plus_reg_l /lc_minus_zero_l_equiv.
rewrite -scal_minus_r lc_scal_l.
move⇒ /(scal_zero_reg_r_R _ _ (INR_inv_n0 Hk)) /Hvtx
    /minus_zero_equiv /(mapF_inj _ _ _ INR_eq).
apply Adk_inj.
Qed.

[RR9557v1] Lem 1592, p. 98.
Lemma vtx_cur_node_cur :
   vtx_cur, (0 < k)%coq_nat invalF vtx_cur (node_cur vtx_cur).
Proof.
intros vtx_cur Hk ipk; destruct (ord_eq_dec ipk ord0) as [-> | Hipk].
ord0;
    rewrite node_cur_eq (Adk_0 d k) scal_zero_r lc_zero_l plus_zero_r; easy.
assert (H :
    scal (/ INR k) (fun j : 'I_dINR k × kronR (lower_S Hipk) j) =
      fun jkronR (lower_S Hipk) j).
  apply extF; intros j; rewrite fct_scal_r_eq;
  unfold scal; simpl; unfold mult; simpl.
  rewrite -Rmult_assoc (Rinv_l _ (INR_n0 _))// Rmult_1_l; easy.
(Adk_inv d k (itemF d k (lower_S Hipk))).
rewrite node_cur_eq Adk_kron (lc_eq_l _ (kronR (lower_S Hipk)))//.
rewrite lc_kron_l_in_r fct_minus_eq constF_correct liftF_lower_S.
rewrite plus_minus_r; easy.
Qed.

Lemma vtx_ref_node_ref : (0 < k)%coq_nat invalF vtx_ref node_ref.
Proof. rewrite node_ref_cur; apply vtx_cur_node_cur. Qed.

End Vertices_Nodes_dk.

Section Vertices_Nodes_1k.

Context {k : nat}.
Context {vtx_cur : 'R^{2,1}}.
Hypothesis Hvtx : affine_independent vtx_cur.

Lemma vtx_cur_1k_neq :
  affine_independent vtx_cur vtx_cur ord1 ord0 vtx_cur ord0 ord0.
Proof.
move⇒ /lin_indep_1_equiv.
rewrite fct_minus_eq liftF_S_0 constF_correct -contra_equiv fun_ext_equiv.
intros H i; rewrite I_1_is_unit; apply minus_zero_equiv; easy.
Qed.

Lemma node_cur_neq :
   (i : 'I_(pbinom 1 k).+1) j,
    node_cur vtx_cur i ord0 skipF (node_cur vtx_cur) i j ord0.
Proof.
move: Hvtx; move⇒ /node_cur_inj H; specialize (H k); move: H ⇒ /inj_contra H.
intros i j; specialize (H _ _ (not_eq_sym (skip_ord_correct_m i j))).
contradict H; apply extF; intros l; rewrite I_1_is_unit; easy.
Qed.

End Vertices_Nodes_1k.

Section Vertices_Nodes_d1.

Context {d : nat}.

Lemma vtx_cur_node_cur_d1 :
   (vtx_cur : 'R^{d.+1,d}),
    vtx_cur = castF (pbinomS_1_r d) (node_cur vtx_cur).
Proof.
intros vtx_cur; apply extF; intros i; rewrite node_cur_eq; unfold castF.
rewrite Rinv_1 scal_one; destruct (ord_eq_dec i ord0) as [-> | Hi].
rewrite cast_ordS_0// Adk_0 lc_zero_l plus_zero_r; easy.
assert (Hi' : ¬ (i < 1)%coq_nat) by now apply ord_n0_nlt_equiv.
rewrite Ad1_eq; unfold castF.
rewrite concatF_correct_r mapF_itemF_0// itemF_kron_eq; simpl.
replace (fun j ⇒ 1 × _) with (fun j : 'I_dkronR (i - 1) j)
    by now apply extF; intros j; rewrite Rmult_1_l.
rewrite (lc_kron_l_in_r (lower_S Hi)) fct_minus_eq liftF_lower_S.
rewrite constF_correct plus_minus_r; easy.
Qed.

Lemma vtx_ref_node_ref_d1 : vtx_ref = castF (pbinomS_1_r d) node_ref.
Proof. rewrite node_ref_cur; apply vtx_cur_node_cur_d1. Qed.

End Vertices_Nodes_d1.

Section Sub_vertices_Def.

Context {d : nat}.
Variable k : nat.

Variable vtx_cur : 'R^{d.+1,d}.

[RR9557v1] Def 1594, Eq (9.88), p. 98. sub_vtx 0 = \underline{v}0 = v_0, sub_vtx i = \underline{v}i = a{(k-1)e^i}, for i <> 0.
Definition sub_vtx : 'R^{d.+1,d} :=
  fun imatch (ord_eq_dec i ord0) with
    | left _vtx_cur ord0
    | right Hinode_cur vtx_cur (Adk_inv d k (itemF d k.-1 (lower_S Hi)))
    end.

Lemma sub_vtx_0 : {i}, i = ord0 sub_vtx i = vtx_cur ord0.
Proof. intros; unfold sub_vtx; case (ord_eq_dec _ _); easy. Qed.

Lemma sub_vtx_S :
   {i} (Hi : i ord0),
     sub_vtx i = node_cur vtx_cur (Adk_inv d k (itemF d k.-1 (lower_S Hi))).
Proof.
intros; unfold sub_vtx; case (ord_eq_dec _ _);
    [| intros; do 3 f_equal; apply ord_inj]; easy.
Qed.

Lemma sub_vtx_d_1 : i, k = 1%nat sub_vtx i = vtx_cur ord0.
Proof.
intros i Hk; unfold sub_vtx; destruct (ord_eq_dec _ _) as [Hi | Hi]; [easy |].
subst; rewrite node_cur_eq Adk_inv_correct_r; [| rewrite sum_itemF; lia].
rewrite lc_scal_l mapF_itemF_0// lc_itemF_l
    scal_zero_l scal_zero_r plus_zero_r; easy.
Qed.

[RR9557v1] Lem 1595, p. 99. \underline{v}i = 1/k v_0 + (k-1)/k v_i.
Lemma sub_vtx_eq :
  (0 < k)%coq_nat
  liftF_S sub_vtx =
    (constF d (scal (/ INR k) (vtx_cur ord0)) +
      scal (INR k.-1 / INR k) (liftF_S vtx_cur))%M.
Proof.
intros Hk; apply extF; intros i; unfold liftF_S.
rewrite (sub_vtx_S (lift_S_not_first i)); unfold node_cur.
rewrite -lc_ones_r Adk_inv_correct_r; [| rewrite sum_itemF; lia].
rewrite !lc_scal_l mapF_itemF_0// !lc_itemF_l.
rewrite ones_correct liftF_lower_S !scal_assoc scal_one_r//.
unfold mult; simpl; rewrite Rmult_comm -(INR_inv_decomp Hk).
unfold Rdiv; rewrite minus_plus_r; easy.
Qed.

End Sub_vertices_Def.

Section Sub_vertices_Facts.

Context {d k : nat}.
Hypothesis Hk : (1 < k)%coq_nat.

Context {vtx_cur : 'R^{d.+1,d}}.
Hypothesis Hvtx : affine_independent vtx_cur.

[RR9557v1] Lem 1597, pp. 99-100.
Lemma sub_vtx_affine_independent : affine_independent (sub_vtx k vtx_cur).
Proof.
assert (Hk0 : (0 < k)%coq_nat) by lia.
intros L; rewrite sub_vtx_0// sub_vtx_eq// plus_comm -plus_minus_assoc.
rewrite -constF_minus -{2}(scal_one (vtx_cur _)) -scal_minus_l INR_inv_pred//.
rewrite constF_scal scal_opp_l -minus_eq -scal_minus_r lc_scal_r.
assert (Hk1 : invertible (INR k.-1 / INR k)).
  unfold Rdiv; apply invertible_equiv_R, is_integral_R;
      [apply INR_n0 | apply INR_inv_n0]; lia.
move⇒ /(scal_zero_reg_r _ _ Hk1); apply Hvtx.
Qed.

End Sub_vertices_Facts.

Section Sub_nodes_Def.

Context {d : nat}.
Variable k : nat.

Variable vtx_cur : 'R^{d.+1,d}.

[RR9557v1] Lem 1598, p. 100.
Definition sub_node : 'R^{(pbinom d k.-1).+1,d} :=
  node_cur (sub_vtx k vtx_cur).

End Sub_nodes_Def.

Section Sub_nodes_Facts.

Context {d k : nat}.
Hypothesis Hk : (1 < k)%coq_nat.

Variable vtx_cur : 'R^{d.+1,d}.

[RR9557v1] Lem 1598, p. 100.
Lemma sub_node_cur_eq :
  sub_node k vtx_cur = widenF (pbinomS_monot_pred d k) (node_cur vtx_cur).
Proof.
assert (Hk0 : (0 < k)%coq_nat) by lia.
apply extF; intros i; unfold sub_node, node_cur, widenF.
rewrite sub_vtx_0// sub_vtx_eq//. rewrite lc_plus_r lc_constF_r.
rewrite plus_assoc scal_assoc -scal_distr_r -!scal_sum_distr_l.
rewrite lc_scal_r -lc_scal_l scal_assoc -minus2.
rewrite -mult_comm_R mult_assoc -mult_minus_distr_r; unfold scal; simpl.
repeat f_equal.
2,4: apply Adk_previous_layer; easy.
1,2: rewrite INR_pred//; unfold mult; simpl; unfold minus, opp, plus; simpl;
    field; rewrite -INR_pred//; split; apply INR_n0; lia.
Qed.

End Sub_nodes_Facts.