CFEM.shape: Lagrange shape functions, and some supporting theory
From mathcomp Require Import all_ssreflect ssralg ssrnum archimedean finfun.
From mathcomp Require Import all_algebra all_field all_analysis all_reals.
Import Order.TTheory GRing.Theory Num.Theory.
Unset Implicit Arguments.
Local Open Scope R_scope.
Local Open Scope order_scope.
Local Open Scope ring_scope.
From mathcomp.algebra_tactics Require Import ring lra.
Section S.
Context {R : realType}.
Definition continuously_differentiable [n] (f: 'rV[R]_n → R) : Prop :=
(∀ x, differentiable f x) ∧ ∀ i, continuous (derive f ^~ ('e_i)).
Definition convex_combination' [d n] (vtx: 'M[R]_(d,n)) (y: 'cV_d) : Prop :=
∃ x: 'cV_n, (∀ i, x i 0 ≥ 0)
∧ \sum_i x i 0 = 1
∧ vtx ×m x = y.
Definition convex_combination [d n] (vtx: 'M[R]_(d,n)) (y: 'cV_d) : Prop :=
∃ x: 'cV_n, (∀ i, x i 0 ≥ 0)
∧ col_mx vtx (const_mx (m:=1) 1) ×m x = col_mx y (const_mx 1).
Lemma convex_combination_e: ∀ d n vtx y,
@convex_combination d n vtx y ↔ @convex_combination' d n vtx y.
Definition convex_hull [n d] (vtx: 'M[R]_(d,n)) := sig (convex_combination vtx).
End S.
Import GRing.
Module Shape.
Section S.
Context {R : realType}.
Record shape : Type := {
d : nat;
nsh: nat;
θ: 'rV_d → 'rV_nsh;
vtx: 'M[R]_(nsh,d);
K := convex_hull vtx;
lagrangian: ∀ i j, θ (row i vtx) 0 j = if i==j then 1 else 0;
diff: ∀ j: 'I_nsh, continuously_differentiable (fun i ⇒ θ i 0 j)
}.
End S.
End Shape.
Definition single {R: realType} (x: R) := @const_mx R 1 1 x.
From Stdlib Require Import Lia.
From Stdlib Require Import FunctionalExtensionality.
Ltac case_splitP j :=
tryif clearbody j then fail "case_splitP requires a variable, but got a local definition" j
else tryif is_var j then idtac else fail "case_splitP requires a variable, but got" j;
match type of j with 'I_(addn ?a ?b) ⇒
let i := fresh "j" in let H := fresh in
destruct (splitP j) as [i H | i H];
[replace j with (@lshift a b i); [ | apply ord_inj; simpl; lia]
|replace j with (@rshift a b i); [ | apply ord_inj; simpl; lia]];
clear j H; rename i into j
end.
Lemma eq_differentiable: ∀ {K: numDomainType} {V W: normedModType K}
(f g: NormedModule.sort V → NormedModule.sort W) x,
f =1 g → differentiable f x = differentiable g x.
Lemma eq_continuously_differentiable: ∀ {R: realType} [n] (f g: ('rV_n → Real.sort R) ),
f =1 g → continuously_differentiable f = continuously_differentiable g.
Lemma eq_locked_continuously_differentiable: ∀ {R: realType} [n] (f g: ('rV_n → Real.sort R) ),
f =1 g → locked continuously_differentiable n f = locked continuously_differentiable n g.
Section S.
Context {R : realType}.
Lemma derive_comp_mx: ∀ [n] (f: R → Real.sort R),
(∀ x, differentiable f x) →
∀ i : 'I_n.+1,
'D_('e_i) (f \o (fun m : 'rV[R]_n.+1 ⇒ m 0 i)) = ('D_1 f) \o (fun m: 'rV[R]_n.+1 ⇒ m 0 i).
Lemma derive_comp_mx_neq: ∀ [n] (f: R → Real.sort R),
(∀ x, differentiable f x) →
∀ i j : 'I_n.+1,
j != i →
'D_('e_i) (f \o (fun m : 'rV[R]_n.+1 ⇒ m 0 j)) = (fun⇒0).
Definition at00 (m: 'cV[R]_1) := m 0 0.
Lemma derive_comp_mx1: ∀ (f: Real.sort R → Real.sort R),
(∀ x, differentiable f x) →
'D_(1%:M) (f \o at00) = ('D_1 f) \o at00.
Inductive is_multivariate_polynomial [d] : ('rV[R]_d → R) → Prop :=
| IMP_const: ∀ (c: R), is_multivariate_polynomial (fun _ ⇒ c)
| IMP_sum: ∀ a b, is_multivariate_polynomial a → is_multivariate_polynomial b →
is_multivariate_polynomial (fun x ⇒ a x + b x)
| IMP_opp: ∀ a, is_multivariate_polynomial a →
is_multivariate_polynomial (fun x ⇒ - a x)
| IMP_prod: ∀ a b, is_multivariate_polynomial a → is_multivariate_polynomial b →
is_multivariate_polynomial (fun x ⇒ a x × b x)
| IMP_mono1: ∀ i (k: nat), is_multivariate_polynomial (fun x ⇒ x 0 i)
| IMP_mono: ∀ i (k: nat), is_multivariate_polynomial (fun x ⇒ x 0 i ^ k).
Lemma continuously_differentiable_cst: ∀ [d] c, @continuously_differentiable R d.+1 (fun⇒c).
Lemma continuously_differentiable_vacuous: ∀ f, @continuously_differentiable R O f.
Lemma continuously_differentiable_add:
∀ [d] (a b: 'rV_d.+1 → Real.sort R),
continuously_differentiable a →
continuously_differentiable b →
continuously_differentiable (a \+ b)%E.
Lemma continuously_differentiable_mul:
∀ [d] (a b: 'rV_d.+1 → Real.sort R),
continuously_differentiable a →
continuously_differentiable b →
continuously_differentiable (fun x ⇒ a x × b x).
Lemma continously_differentiable_coord:
∀ [d] (i: 'I_d.+1),
continuously_differentiable (fun x : 'rV[R]_d.+1 ⇒ fun_of_matrix x 0 i).
Lemma multivariate_polynomial_continuously_differentiable:
∀ [d] (f: 'rV[R]_d → R),
is_multivariate_polynomial f →
continuously_differentiable f.
Ltac prove_continuously_differentiable :=
let j := fresh "j" in
intro j;
lazymatch goal with
| j: 'I_?n |- continuously_differentiable (fun i ⇒ fun_of_matrix (?f i) 0 j) ⇒ try unfold f
| |- _ ⇒ fail "prove_continuously_differentiable must be applied to a goal of the form, continuously_differentiable (fun i => ...)"
end;
repeat match type of j with context [S (S ?A)] ⇒ change (S (S ?A)) with (1 + S A)%nat in j end;
repeat case_splitP j; rewrite !ord1 /=;
rewrite [continuously_differentiable]lock;
repeat (under eq_locked_continuously_differentiable do [ rewrite /single /at00 ?mxE /split /= ];
repeat (destruct (ltnP _ _); try discriminate));
rewrite -lock;
apply multivariate_polynomial_continuously_differentiable; repeat constructor.
Ltac prove_lagrangian :=
let i := fresh "i" in let j := fresh "j" in
intros i j;
lazymatch goal with
| i: 'I_?n, j: 'I_?n |- fun_of_matrix (?f (row i _)) 0 j = _ ⇒ try unfold f
end;
rewrite /single /at00 !mxE;
repeat match type of i with context [S (S ?A)] ⇒ change (S (S ?A)) with (1 + S A)%nat in i end;
repeat match type of j with context [S (S ?A)] ⇒ change (S (S ?A)) with (1 + S A)%nat in j end;
repeat case_splitP j; repeat case_splitP i;
repeat (rewrite ?ord1 ?mxE /= /split; repeat (destruct (ltnP _ _); try discriminate));
rewrite ?mxE /=; try nra.
Definition shapes1dP1_function : 'rV_1 → 'rV_(1 + 1) :=
((fun x ⇒ row_mx (single (1/2 × (1-x))) (single (1/2 × (1 + x)))) \o at00).
Definition shapes1dP1_vertices : 'cV[R]_2 := col_mx (single (-1)) (single 1).
Definition shapes1dP1 : @Shape.shape R.
Defined.
Definition shapes1dP2_vertices : 'cV[R]_3 := col_mx (single (-1)) (col_mx (single 0) (single 1)).
Definition shapes1dP2_function : 'rV_1 → 'rV_3 :=
((fun x ⇒ row_mx (single (-(1/2)*(1-x)*x)) (row_mx (single ((1-x)*(1+x))) (single ((1/2)*x*(1+x))))) \o at00).
Definition shapes1dP2 : @Shape.shape R.
Defined.
Definition shapes2dT1_vertices : 'M[R]_(3,2) :=
col_mx (row_mx (single 0) (single 0))
(col_mx (row_mx (single 1) (single 0))
(row_mx (single 0) (single 1))).
Definition shapes2dT1_function : 'rV[R]_2 → 'rV[R]_3 :=
(fun x ⇒ row_mx (single (1 - x 0 0 - x 0 1))
(row_mx (single (x 0 0))
(single (x 0 1)))).
Definition shapes2dT1 : @Shape.shape R.
Defined.
End S.
From mathcomp Require Import all_algebra all_field all_analysis all_reals.
Import Order.TTheory GRing.Theory Num.Theory.
Unset Implicit Arguments.
Local Open Scope R_scope.
Local Open Scope order_scope.
Local Open Scope ring_scope.
From mathcomp.algebra_tactics Require Import ring lra.
Section S.
Context {R : realType}.
Definition continuously_differentiable [n] (f: 'rV[R]_n → R) : Prop :=
(∀ x, differentiable f x) ∧ ∀ i, continuous (derive f ^~ ('e_i)).
Definition convex_combination' [d n] (vtx: 'M[R]_(d,n)) (y: 'cV_d) : Prop :=
∃ x: 'cV_n, (∀ i, x i 0 ≥ 0)
∧ \sum_i x i 0 = 1
∧ vtx ×m x = y.
Definition convex_combination [d n] (vtx: 'M[R]_(d,n)) (y: 'cV_d) : Prop :=
∃ x: 'cV_n, (∀ i, x i 0 ≥ 0)
∧ col_mx vtx (const_mx (m:=1) 1) ×m x = col_mx y (const_mx 1).
Lemma convex_combination_e: ∀ d n vtx y,
@convex_combination d n vtx y ↔ @convex_combination' d n vtx y.
Definition convex_hull [n d] (vtx: 'M[R]_(d,n)) := sig (convex_combination vtx).
End S.
Import GRing.
Module Shape.
Section S.
Context {R : realType}.
Record shape : Type := {
d : nat;
nsh: nat;
θ: 'rV_d → 'rV_nsh;
vtx: 'M[R]_(nsh,d);
K := convex_hull vtx;
lagrangian: ∀ i j, θ (row i vtx) 0 j = if i==j then 1 else 0;
diff: ∀ j: 'I_nsh, continuously_differentiable (fun i ⇒ θ i 0 j)
}.
End S.
End Shape.
Definition single {R: realType} (x: R) := @const_mx R 1 1 x.
From Stdlib Require Import Lia.
From Stdlib Require Import FunctionalExtensionality.
Ltac case_splitP j :=
tryif clearbody j then fail "case_splitP requires a variable, but got a local definition" j
else tryif is_var j then idtac else fail "case_splitP requires a variable, but got" j;
match type of j with 'I_(addn ?a ?b) ⇒
let i := fresh "j" in let H := fresh in
destruct (splitP j) as [i H | i H];
[replace j with (@lshift a b i); [ | apply ord_inj; simpl; lia]
|replace j with (@rshift a b i); [ | apply ord_inj; simpl; lia]];
clear j H; rename i into j
end.
Lemma eq_differentiable: ∀ {K: numDomainType} {V W: normedModType K}
(f g: NormedModule.sort V → NormedModule.sort W) x,
f =1 g → differentiable f x = differentiable g x.
Lemma eq_continuously_differentiable: ∀ {R: realType} [n] (f g: ('rV_n → Real.sort R) ),
f =1 g → continuously_differentiable f = continuously_differentiable g.
Lemma eq_locked_continuously_differentiable: ∀ {R: realType} [n] (f g: ('rV_n → Real.sort R) ),
f =1 g → locked continuously_differentiable n f = locked continuously_differentiable n g.
Section S.
Context {R : realType}.
Lemma derive_comp_mx: ∀ [n] (f: R → Real.sort R),
(∀ x, differentiable f x) →
∀ i : 'I_n.+1,
'D_('e_i) (f \o (fun m : 'rV[R]_n.+1 ⇒ m 0 i)) = ('D_1 f) \o (fun m: 'rV[R]_n.+1 ⇒ m 0 i).
Lemma derive_comp_mx_neq: ∀ [n] (f: R → Real.sort R),
(∀ x, differentiable f x) →
∀ i j : 'I_n.+1,
j != i →
'D_('e_i) (f \o (fun m : 'rV[R]_n.+1 ⇒ m 0 j)) = (fun⇒0).
Definition at00 (m: 'cV[R]_1) := m 0 0.
Lemma derive_comp_mx1: ∀ (f: Real.sort R → Real.sort R),
(∀ x, differentiable f x) →
'D_(1%:M) (f \o at00) = ('D_1 f) \o at00.
Inductive is_multivariate_polynomial [d] : ('rV[R]_d → R) → Prop :=
| IMP_const: ∀ (c: R), is_multivariate_polynomial (fun _ ⇒ c)
| IMP_sum: ∀ a b, is_multivariate_polynomial a → is_multivariate_polynomial b →
is_multivariate_polynomial (fun x ⇒ a x + b x)
| IMP_opp: ∀ a, is_multivariate_polynomial a →
is_multivariate_polynomial (fun x ⇒ - a x)
| IMP_prod: ∀ a b, is_multivariate_polynomial a → is_multivariate_polynomial b →
is_multivariate_polynomial (fun x ⇒ a x × b x)
| IMP_mono1: ∀ i (k: nat), is_multivariate_polynomial (fun x ⇒ x 0 i)
| IMP_mono: ∀ i (k: nat), is_multivariate_polynomial (fun x ⇒ x 0 i ^ k).
Lemma continuously_differentiable_cst: ∀ [d] c, @continuously_differentiable R d.+1 (fun⇒c).
Lemma continuously_differentiable_vacuous: ∀ f, @continuously_differentiable R O f.
Lemma continuously_differentiable_add:
∀ [d] (a b: 'rV_d.+1 → Real.sort R),
continuously_differentiable a →
continuously_differentiable b →
continuously_differentiable (a \+ b)%E.
Lemma continuously_differentiable_mul:
∀ [d] (a b: 'rV_d.+1 → Real.sort R),
continuously_differentiable a →
continuously_differentiable b →
continuously_differentiable (fun x ⇒ a x × b x).
Lemma continously_differentiable_coord:
∀ [d] (i: 'I_d.+1),
continuously_differentiable (fun x : 'rV[R]_d.+1 ⇒ fun_of_matrix x 0 i).
Lemma multivariate_polynomial_continuously_differentiable:
∀ [d] (f: 'rV[R]_d → R),
is_multivariate_polynomial f →
continuously_differentiable f.
Ltac prove_continuously_differentiable :=
let j := fresh "j" in
intro j;
lazymatch goal with
| j: 'I_?n |- continuously_differentiable (fun i ⇒ fun_of_matrix (?f i) 0 j) ⇒ try unfold f
| |- _ ⇒ fail "prove_continuously_differentiable must be applied to a goal of the form, continuously_differentiable (fun i => ...)"
end;
repeat match type of j with context [S (S ?A)] ⇒ change (S (S ?A)) with (1 + S A)%nat in j end;
repeat case_splitP j; rewrite !ord1 /=;
rewrite [continuously_differentiable]lock;
repeat (under eq_locked_continuously_differentiable do [ rewrite /single /at00 ?mxE /split /= ];
repeat (destruct (ltnP _ _); try discriminate));
rewrite -lock;
apply multivariate_polynomial_continuously_differentiable; repeat constructor.
Ltac prove_lagrangian :=
let i := fresh "i" in let j := fresh "j" in
intros i j;
lazymatch goal with
| i: 'I_?n, j: 'I_?n |- fun_of_matrix (?f (row i _)) 0 j = _ ⇒ try unfold f
end;
rewrite /single /at00 !mxE;
repeat match type of i with context [S (S ?A)] ⇒ change (S (S ?A)) with (1 + S A)%nat in i end;
repeat match type of j with context [S (S ?A)] ⇒ change (S (S ?A)) with (1 + S A)%nat in j end;
repeat case_splitP j; repeat case_splitP i;
repeat (rewrite ?ord1 ?mxE /= /split; repeat (destruct (ltnP _ _); try discriminate));
rewrite ?mxE /=; try nra.
Definition shapes1dP1_function : 'rV_1 → 'rV_(1 + 1) :=
((fun x ⇒ row_mx (single (1/2 × (1-x))) (single (1/2 × (1 + x)))) \o at00).
Definition shapes1dP1_vertices : 'cV[R]_2 := col_mx (single (-1)) (single 1).
Definition shapes1dP1 : @Shape.shape R.
Defined.
Definition shapes1dP2_vertices : 'cV[R]_3 := col_mx (single (-1)) (col_mx (single 0) (single 1)).
Definition shapes1dP2_function : 'rV_1 → 'rV_3 :=
((fun x ⇒ row_mx (single (-(1/2)*(1-x)*x)) (row_mx (single ((1-x)*(1+x))) (single ((1/2)*x*(1+x))))) \o at00).
Definition shapes1dP2 : @Shape.shape R.
Defined.
Definition shapes2dT1_vertices : 'M[R]_(3,2) :=
col_mx (row_mx (single 0) (single 0))
(col_mx (row_mx (single 1) (single 0))
(row_mx (single 0) (single 1))).
Definition shapes2dT1_function : 'rV[R]_2 → 'rV[R]_3 :=
(fun x ⇒ row_mx (single (1 - x 0 0 - x 0 1))
(row_mx (single (x 0 0))
(single (x 0 1)))).
Definition shapes2dT1 : @Shape.shape R.
Defined.
End S.
some experiments with mathcomp.algebra's Lagrange polynomials
From mathcomp Require Import algebra.qpoly.
Ltac simpl_bigop_index_enum :=
match goal with |- context [@bigop.body ?t1 ?t2 ?idx (index_enum ?s) ?f] ⇒ pattern (@bigop.body t1 t2 idx (index_enum s) f) end;
match goal with |- ?hideme _ ⇒ let hide := fresh "hide" in
set hide := hideme;
rewrite bigop.unlock locked_withE Finite.enum.unlock /= /ord_enum /= /reducebig /Option.apply !insubT /Sub /=;
subst hide; cbv beta
end.
Lemma bigop_2_1: ∀ [R] [op: R → R → R] [idx: R] (f: 'I_2 → R), \big[op/idx]_(i < 2 | i != 0) f i = op (f 1) idx.
Lemma bigop_2_0: ∀ [R] [op: R → R → R] [idx: R] (f: 'I_2 → R), \big[op/idx]_(i < 2 | i != 1) f i = op (f 0) idx.
Section S.
Context {R : realType}.
Definition R_of_nat : nat → Nmodule.sort (reals_Real__to__GRing_Nmodule R) := @natmul R 1.
Lemma injective_R_of_nat : injective R_of_nat.
Hint Resolve injective_R_of_nat : core.
Definition the_points : nat → Nmodule.sort (reals_Real__to__GRing_Nmodule R) := fun i ⇒ -1 + natmul 2 i.
Lemma injective_the_points: injective the_points.
Hint Resolve injective_the_points : core.
Definition lagrange_2 := lagrange 2 the_points.
Definition lagrange_2_0 : {poly_2 R} := tnth lagrange_2 0.
Lemma sub_natmulE1: ∀ (x: R) (i j: nat), is_true (j≤i)%N → add (natmul x i) (opp (natmul x j)) = (natmul x (subn i j)).
Lemma sub_natmulE2: ∀ (x: R) (i j: nat), is_true (i≤j)%N → add (natmul x i) (opp (natmul x j)) = opp (natmul x (subn j i)).
Ltac rsimpl :=
repeat (rewrite ?opprK ?invrN ?subr0 ?invr1 ?subr0 ?mul1r ?mulr1 ?mulN1r ?opprB ?addSn ?add0n ?mulN1r;
try (rewrite modn_small; [ | lia]);
try (rewrite sub_natmulE1; [ | lia]);
try (rewrite sub_natmulE2; [ | lia])).
Goal ∀ x, horner (tnth lagrange_2 0) x = (2^-1) × (1 - x).
Goal ∀ x, horner (tnth lagrange_2 1) x = (2^-1) × (1 + x).
Goal ∀ i j : 'I_2, horner (tnth lagrange_2 i) (the_points j) = (i==j)%:R.
Goal ∀ x, horner (tnth (lagrange 1 R_of_nat) 0) x = 1.
Goal ∀ x, horner (tnth (lagrange 3 R_of_nat) 0) x = 2^-1 × (x - 1) × (x - 2) .
Goal ∀ x, horner (tnth (lagrange 3 R_of_nat) 1) x = - x × (x - 2).
Goal ∀ x, horner (tnth (lagrange 3 R_of_nat) 2) x = 2^-1 × x × (x - 1).
End S.