Library FEM.poly_LagP1k

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 geom_simplex geom_transf_affine.
From FEM Require Import multi_index poly_Pdk poly_LagPd1_ref.

Local Open Scope R_scope.

Section Poly_LagP1k.

Context {k : nat}.

[RR9557v1] Def 1448, p. 58.
LagP1k i = Prod{j<>i} (x - node j) / Prod{j<>i} (node i - node j).
Definition LagP1k_cur_aux vtx_cur (i : 'I_(pbinom 1 k).+1) (x : 'R^1) : R :=
  prod_R (constF _ (x ord0) - (fun jskipF (node_cur vtx_cur) i j ord0))%G.

Definition LagP1k_cur vtx_cur : '('R^1 R)^((pbinom 1 k).+1) :=
  fun i xLagP1k_cur_aux vtx_cur i x / LagP1k_cur_aux vtx_cur i (node_cur vtx_cur i).

Definition LagP1k_ref : '('R^1 R)^((pbinom 1 k).+1) :=
  LagP1k_cur vtx_ref.

[RR9557v1] Lem 1470, pp. 61-62.
Lemma LagP1k_cur_eq : {vtx_cur : 'R^{2,1}}
  (Hvtx: affine_independent vtx_cur),
 LagP1k_cur vtx_cur = fun i xLagP1k_ref i (T_geom_inv Hvtx x).
Proof.
intros vtx_cur Hvtx.
apply fun_ext; intros i.
unfold LagP1k_ref.
apply fun_ext; intros x.
unfold LagP1k_cur, LagP1k_cur_aux.
rewrite !prod_R_div; f_equal.
apply fun_ext; intros j.
rewrite !fct_minus_eq.
rewrite !constF_correct.
rewrite -{1}(T_geom_inv_can_r Hvtx x).
pose (x_ref := T_geom_inv Hvtx x).
fold x_ref.
rewrite -{1}(T_geom_inv_can_r Hvtx (node_cur vtx_cur i)).
rewrite -(T_geom_inv_can_r Hvtx (skipF (node_cur vtx_cur) i j)).
rewrite !T_geom_inv_transports_node; try easy.
rewrite -node_ref_cur.
rewrite !T_geom_1k; try easy.
rewrite !fct_plus_eq !fct_scal_r_eq !fct_minus_eq.
rewrite -!minus_plus_r_eq.
rewrite -!scal_minus_distr_r.
unfold scal, minus, plus, mult, opp; simpl.
unfold mult; simpl; unfold opp; simpl; unfold skipF.
unfold Rdiv; rewrite Rmult_assoc; f_equal.
rewrite 2!Rinv_mult Rmult_comm 2!Rmult_assoc; f_equal.
rewrite Rinv_inv; field.
apply Rminus_eq_contra.
apply vtx_cur_1k_neq; easy.
apply inF_not; intros i0.
apply nesym.
apply Rminus_eq_contra.
apply node_cur_neq.
apply vtx_ref_affine_independent.
apply inF_not; intros i0.
apply nesym.
apply Rminus_eq_contra.
apply node_cur_neq; easy.
Qed.

[RR9557v1] Lem 1469, p. 61.
See also Lem 1449, Eq (8.5), p. 58.
Lemma LagP1k_cur_kron_node : (i : 'I_(pbinom 1 k).+1) {vtx_cur : 'R^{2,1}},
   affine_independent vtx_cur
  (LagP1k_cur vtx_cur)^~ (node_cur vtx_cur i) = kronR i.
Proof.
intros i vtx_cur Hvtx.
apply extF; intros j.
unfold LagP1k_cur, LagP1k_cur_aux.
destruct (ord_eq_dec i j) as [Hj | Hj].
rewrite prod_R_div.
rewrite Hj kronR_is_1; try easy.
rewrite prod_R_one_compat; try easy.
apply extF; intros l.
rewrite zeroF; unfold zero; simpl.
apply Rinv_r.
rewrite fct_minus_eq.
apply Rminus_eq_contra.
apply node_cur_neq; try easy.
apply inF_not; intros i0.
apply nesym.
apply Rminus_eq_contra.
apply node_cur_neq; easy.
rewrite prod_R_div.
rewrite kronR_is_0.
rewrite prod_R_zero; try easy.
unfold inF.
(insert_ord Hj).
apply eq_sym.
unfold Rdiv.
assert (H : (constF (pbinom 1 k) (node_cur vtx_cur i ord0) -
 (skipF (node_cur vtx_cur) j)^~ ord0)%G (insert_ord Hj) = 0).
apply Rminus_diag_eq.
unfold constF, skipF.
f_equal.
rewrite skip_insert_ord; easy.
rewrite H; lra.
contradict Hj.
apply ord_inj; easy.
apply inF_not; intros i0.
apply nesym.
apply Rminus_eq_contra.
apply node_cur_neq; easy.
Qed.

Lemma LagP1k_cur_kron_node_alt : (i j : 'I_(pbinom 1 k).+1) {vtx_cur : 'R^{2,1}},
  affine_independent vtx_cur
  LagP1k_cur vtx_cur j (node_cur vtx_cur i) =
    kronR i j.
Proof.
intros i j vtx_cur Hvtx.
generalize LagP1k_cur_kron_node.
intros H.
specialize (H i vtx_cur Hvtx).
rewrite fun_ext_equiv in H; easy.
Qed.

Lemma LagP1k_cur_lin_indep : {vtx_cur : 'R^{2,1}},
  affine_independent vtx_cur
    lin_indep (LagP1k_cur vtx_cur).
Proof.
movevtx_cur Hvtx L /fun_ext_equiv HL.
rewrite -(lc_kron_r L).
apply extF; intros i.
specialize (HL (node_cur vtx_cur i)).
rewrite zeroF.
rewrite fct_zero_eq in HL.
rewrite fct_lc_r_eq in HL.
rewrite LagP1k_cur_kron_node in HL; try easy.
rewrite -HL.
rewrite fct_lc_r_eq.
f_equal.
apply fun_ext; intros l.
apply kron_sym.
Qed.

Lemma LagP1k_ref_kron_node : (i j : 'I_(pbinom 1 k).+1),
  LagP1k_ref j (node_ref i) = kronR i j.
Proof.
intros i j.
unfold LagP1k_ref; rewrite node_ref_cur.
apply LagP1k_cur_kron_node_alt.
apply vtx_ref_affine_independent.
Qed.

Lemma LagP1k_ref_lin_indep : lin_indep LagP1k_ref.
Proof.
unfold LagP1k_ref.
apply LagP1k_cur_lin_indep.
apply vtx_ref_affine_independent; easy.
Qed.

End Poly_LagP1k.

Section Poly_LagP11.

Lemma LagP1k_cur_0 : (i : 'I_(pbinom 1 1).+1) (x : 'R^1),
   i = ord0
    LagP1k_cur vtx_ref i x = 1 - sum x.
Proof.
intros i x Hi.
rewrite Hi sum_1.
unfold LagP1k_cur, LagP1k_cur_aux.
rewrite prod_R_div.
2 : {apply inF_not; intros j.
apply nesym.
apply Rminus_eq_contra.
apply node_cur_neq; auto.
apply vtx_ref_affine_independent. }
rewrite -node_ref_cur; unfold node_ref.
rewrite prod_R_singleF_alt.
rewrite castF_id.
rewrite !fct_minus_eq !constF_correct.
replace (skipF (fun (ipk : 'I_(pbinom 1 1).+1) (i0 : 'I_1) ⇒
  INR (Adk 1 1 ipk i0) / INR 1)
   ord0 ord0 ord0) with 1.
rewrite A1k_eq; simpl.
unfold minus, opp, plus, one; simpl; field.
unfold skipF; rewrite Rdiv_1_r; easy.
Qed.

Lemma LagP1k_cur_1 : (i : 'I_(pbinom 1 1).+1) (x : 'R^1)
  (Hi : cast_ord (pbinomS_1_r 1) i ord0),
    LagP1k_cur vtx_ref i x = x (lower_S Hi).
Proof.
intros i x Hi.
assert (i = ord_max).
destruct (ord2_dec i) as [Hi1 | Hi1]; try easy.
contradict Hi; rewrite Hi1; apply ord_inj; easy.
unfold LagP1k_cur, LagP1k_cur_aux.
rewrite prod_R_div.
2 : {apply inF_not; intros j.
apply nesym.
apply Rminus_eq_contra.
apply node_cur_neq; auto.
apply vtx_ref_affine_independent. }
rewrite -node_ref_cur; unfold node_ref.
rewrite prod_R_singleF_alt.
rewrite castF_id.
rewrite !fct_minus_eq !constF_correct.
replace (skipF
   (fun (ipk : 'I_(pbinom 1 1).+1) (i0 : 'I_1) ⇒
    INR (Adk 1 1 ipk i0) / INR 1) i ord0 ord0) with 0.
rewrite !minus_zero_r Rdiv_1_r A1k_eq constF_correct.
replace (INR i) with 1.
rewrite Rdiv_1_r; f_equal.
apply ord_inj; unfold ord0, lower_S; simpl.
rewrite H; simpl; auto with arith.
rewrite H; simpl.
replace 1 with (INR 1); try easy; f_equal.
unfold skipF; rewrite Rdiv_1_r A1k_eq constF_correct H; easy.
Qed.

Lemma LagP11_ref_eq :
  LagP1k_ref =
    castF (eq_sym (pbinomS_1_r 1)) LagPd1_ref.
Proof.
unfold LagP1k_ref.
unfold castF; rewrite eq_sym_involutive.
apply extF; intros i.
apply fun_ext; intros x; simpl.
destruct (ord_eq_dec i ord0) as [H | H].
rewrite LagPd1_ref_0.
apply LagP1k_cur_0; easy.
rewrite H.
apply ord_inj; easy.
rewrite LagPd1_ref_S.
contradict H.
apply (f_equal (@nat_of_ord 2)) in H.
simpl in H.
apply ord_inj; easy.
intros H0; apply LagP1k_cur_1.
Qed.

End Poly_LagP11.

Section P1k_Facts.

Context {k : nat}.

Lemma LagP1k_ref_is_poly : inclF (@LagP1k_ref k) (Pdk 1 k).
Proof.
intros i.
unfold LagP1k_ref, LagP1k_cur.
eapply Pdk_ext.
2: apply lin_span_scal_closed.
intros x; unfold Rdiv; rewrite Rmult_comm; easy.
unfold LagP1k_cur_aux.
apply Pdk_mult_iter with (k:=fun _S O).
intros j.
eapply Pdk_ext.
2: apply lin_span_minus_closed.
2: apply (Pdk_mul_var _ (fun _ ⇒ 1) ord0).
2: apply Pd0_eq; easy.
2: apply Pdk_monot with O; auto with arith.
2: apply Pd0_eq.
intros y; unfold minus; simpl; unfold opp, plus, fct_opp; simpl; unfold fct_opp, fct_plus; simpl.
rewrite constF_correct Rmult_1_r; easy.
easy.
rewrite sum_constF_nat Nat.mul_1_r.
rewrite pbinom_1_l; apply Nat.le_refl.
Qed.

[RR9557v1] Lem 1455, p. 59.
See also Lem 1449, p. 58.
Lemma LagP1k_ref_basis : basis (Pdk 1 k) (@LagP1k_ref k).
Proof.
apply lin_indep_basis.
apply Pdk_has_dim.
apply LagP1k_ref_is_poly.
apply LagP1k_ref_lin_indep; easy.
Qed.

Lemma P1k_lin_span_LagP1k_cur : {vtx_cur : 'R^{2,1}},
  affine_independent vtx_cur
  Pdk 1 k = lin_span (@LagP1k_cur k vtx_cur).
Proof.
intros vtx_cur Hvtx.
rewrite LagP1k_cur_eq; try easy.
rewrite (Pdk_am_comp_basis k LagP1k_ref (T_geom_inv Hvtx)); try easy.
apply LagP1k_ref_basis.
now apply T_geom_inv_am.
apply f_inv_bij.
Qed.

Lemma LagP1k_cur_basis : {vtx_cur},
  affine_independent vtx_cur
  basis (Pdk 1 k) (@LagP1k_cur k vtx_cur).
Proof.
intros vtx_cur Hvtx.
rewrite (P1k_lin_span_LagP1k_cur Hvtx); try easy.
apply basis_lin_span_equiv.
apply LagP1k_cur_lin_indep; easy.
Qed.

[RR9557v1] Lem 1469, p. 61.
See also Lem 1449, Eq (8.6), p. 58.
Lemma LagP1k_cur_sum_1 : {vtx_cur : 'R^{2,1}} (x : 'R^1),
    affine_independent vtx_cur
    sum (fun i : 'I_(pbinom 1 k).+1LagP1k_cur vtx_cur i x) = 1.
Proof.
intros vtx_cur x Hvtx.
assert (H:Pdk 1 k 1%K).
apply Pdk_monot with O; try auto with arith.
apply Pd0_eq; easy.
rewrite (P1k_lin_span_LagP1k_cur Hvtx) in H.
inversion H.
rewrite -lc_ones_l.
generalize (fun_ext_rev H1); intros H2; specialize (H2 x).
unfold one in H2; simpl in H2; unfold fct_one in H2.
rewrite fct_lc_r_eq in H2.
unfold one in H2; simpl in H2.
rewrite -H2.
apply lc_ext_l; intros i; unfold ones.
rewrite constF_correct; apply eq_sym.
generalize (fun_ext_rev H1 (node_cur vtx_cur i)); intros H3.
unfold one in H3; simpl in H3; unfold fct_one in H3; simpl in H3.
rewrite -H3.
rewrite fct_lc_r_eq.
erewrite lc_eq_r.
2: rewrite LagP1k_cur_kron_node; try easy.
rewrite -{1}(lc_kron_r L).
rewrite fct_lc_r_eq.
apply lc_ext_r; intros j.
apply kron_sym.
Qed.

Lemma LagP1k_ref_sum_1 : x,
  sum ((@LagP1k_ref k)^~ x) = 1.
Proof.
intros x.
unfold LagP1k_ref.
apply LagP1k_cur_sum_1.
apply (vtx_ref_affine_independent); easy.
Qed.

End P1k_Facts.