Library one_d_wave_equation

Require Export Differential.
Require Export R_n.
Require Export ZArith.
Open Local Scope nat_scope.
Open Local Scope R_scope.

Section Wave.
Variable xmin xmax T:R.

Hypothesis xminLtxmax: xmin < xmax.
Hypothesis T_pos: 0 < T.

Variable partial_derive_firstvar: (R * R -> R) -> (R * R -> R).
Variable partial_derive_secondvar: (R * R -> R) -> (R * R -> R).

Definition Xjk (dX : R * R) (j: Z) (k:nat) : R * R :=
  (xmin+(IZR j)*(fst dX), (INR k)*(snd dX)).

Definition Xxt (X dX:R*R) : R*R :=
  Xjk dX (jx (fst X-xmin) (fst dX)) (kt (snd X) (snd dX)).

Notation ni x := (ni xmin xmax x).
Notation finite_dot_product_dx := (fun dx u v => (finite_dot_product_dx xmin xmax dx u v)).
Notation finite_squared_norm_dx := (fun dx u => (finite_squared_norm_dx xmin xmax dx u)).
Notation finite_norm_dx := (fun dx u => (finite_norm_dx xmin xmax dx u)).

Variable c:R.
Variable u0 u1: R->R.
Variable f : R*R -> R.
Variable usol : R * R -> R.

Variable u0h u1h: R*R -> Z-> R.
Variable fh: R*R -> Z -> nat -> R.
Variable unh : R * R -> Z -> nat -> R.
Variable eps:R.

Definition L (u : R * R -> R) (X : R * R) : R :=
  (partial_derive_secondvar_n partial_derive_secondvar 2%nat u) X -
  c ^ 2 * (partial_derive_firstvar_n partial_derive_firstvar 2%nat u) X.

Definition L1 (u: R*R -> R) (X:R*R) : R :=
   (partial_derive_secondvar u) X.

Definition L0 (u: R*R -> R) (X:R*R) : R :=
  u X.

Definition is_solution (u: R*R -> R) :=
  (forall X:R*R, xmin <= fst X <= xmax -> L u X=f X) /\
  (forall x:R, xmin <= x <= xmax -> L1 u (x,0) = u1 x) /\
  (forall x:R, xmin <= x <= xmax -> L0 u (x,0) = u0 x) /\
  (forall t:R, u(xmin,t) = 0) /\
  (forall t:R, u(xmax,t) = 0).

Hypothesis usol_is_sol: is_solution usol.

Definition discretization (u : R * R -> R) (dX : R * R) : Z -> nat-> R :=
  fun j k => u (Xjk dX j k).

Definition ubar (dX : R * R) (j:Z) (k : nat) : R :=
   discretization usol dX j k.

Definition Ah (unh : Z -> R) (dx : R) (j : Z) : R :=
  - (c *c) * (unh (j + 1)%Z - 2 * unh j + unh (j - 1)%Z) / (dx *dx).

Definition Lh (uh: Z -> nat -> R) (dX : R * R) (j:Z) (n : nat) : R :=
   match n with
     O => uh j 0%nat
   | 1 => (uh j 1%nat - uh j 0%nat)/ (snd dX) + (snd dX) / 2 * (Ah (fun j' => uh j' 0%nat) (fst dX) j)
   | S (S m) => (uh j (m + 2)%nat - 2 * uh j (m+1)%nat + uh j m) / (snd dX *snd dX) +
       (Ah (fun j' => uh j' (m+1)%nat) (fst dX) j)
   end.

Definition scheme (uh: R*R -> Z -> nat -> R) (dX : R * R) (j:Z) (n : nat) : R :=
   match (Z_lt_le_dec 0 j) with
   | left _ => match (Z_lt_le_dec j (ni (fst dX))) with
         | left _ => match n with
                 O => Lh (uh dX) dX j 0 - u0h dX j
               | 1 => Lh (uh dX) dX j 1 - u1h dX j
               | S (S m) => Lh (uh dX) dX j (S (S m)) - fh dX j (S m)
               end
         | right _ => uh dX j n
         end
   | right _ => uh dX j n
   end.

Definition is_discrete_solution (dX:R*R) := snd dX <> 0 ->
 forall j n, (0 <= j <= ni (fst dX))%Z -> scheme unh dX j n = 0.

Definition LL1h (v: R*R->R) (X dX:R*R) : R:=
   (v (fst X , snd X+snd dX) - v X)/ (snd dX) - partial_derive_secondvar v X -
   (snd dX) / 2 * (c*c) *(v (fst X +fst dX, snd X) - 2*v X +
   v (fst X -fst dX, snd X))/(fst dX*fst dX).

Definition energy (uh : Z -> nat -> R) (dX : R * R) (n : nat) :=
  1/2 *
    finite_squared_norm_dx (fst dX)
      (scalar_multiply_n (1 / snd dX)
        (subtract_n
          (fun j => uh j (n + 1)%nat)
          (fun j => uh j n))) +
  1/2 *
    finite_dot_product_dx (fst dX)
      (fun j => (Ah (fun j' => uh j' n) (fst dX) j))
      (fun j => uh j (n + 1)%nat).

Definition truncation_error := (scheme ubar).

Definition convergence_error (dX : R * R) (j:Z) (n : nat) : R :=
  (ubar dX j n) - (unh dX j n).

Lemma jx_Betw: forall (x:{x| Betw xmin xmax x}) dx, 0 < dx ->
   (0 <= jx (projT1 x-xmin) dx <= ni dx)%Z.

Lemma Xjk_Betw: forall (x:{x| Betw xmin xmax x}) dX i, 0 < fst dX ->
   xmin <= fst (Xjk dX (jx (projT1 x - xmin) (fst dX)) i) <= xmax.

Lemma Xjk_Betw_j: forall dX j i, 0 < fst dX ->
   (0 <= j <= ni (fst dX))%Z ->
   xmin <= fst (Xjk dX j i) <= xmax.

Lemma scheme_eq_1: forall uh dX i n, (i <= 0)%Z ->
  scheme uh dX i n = uh dX i n.

Lemma scheme_eq_2: forall uh dX i n, (ni (fst dX) <= i)%Z ->
  scheme uh dX i n = uh dX i n.

Lemma scheme_eq_3: forall uh dX j n, (0 < j < ni (fst dX))%Z ->
  scheme uh dX j n = match n with
                 O => Lh (uh dX) dX j 0 - u0h dX j
               | 1 => Lh (uh dX) dX j 1 - u1h dX j
               | S (S m) => Lh (uh dX) dX j (S (S m)) - fh dX j (S m)
               end.

Lemma Ah_is_linear: forall uh1 uh2 dx j,
  Ah (fun j0 => uh1 j0 - uh2 j0) dx j = Ah uh1 dx j - Ah uh2 dx j.

Lemma Lh_is_linear: forall uh1 uh2 dX j n,
  Lh (fun j0 n0 => uh1 j0 n0 - uh2 j0 n0) dX j n = Lh uh1 dX j n - Lh uh2 dX j n.

Lemma Green_formula: forall dx u v,
  0 < dx ->
 v 0%Z = 0 -> v (ni dx) = 0 ->
 finite_dot_product_dx dx (fun n => u (n+1)%Z -2* u n + u (n-1)%Z) v
  = - finite_dot_product_dx dx (fun n => u (n+1)%Z - u n)
                                                     (fun n => v (n+1)%Z - v n).

Lemma Ah_is_symmetrical :
  forall dx u v,
  0 < dx ->
  u 0%Z = 0 -> u (ni dx) = 0 -> v 0%Z = 0 -> v (ni dx) = 0 ->
    finite_dot_product_dx dx (Ah u dx) v
   = finite_dot_product_dx dx u (Ah v dx).

Lemma Ah_eq: forall uh1 uh2 dx, (forall j, (0 <= j <= ni dx)%Z -> uh1 j=uh2 j) ->
  forall j, (0 < j < ni dx)%Z -> Ah uh1 dx j= Ah uh2 dx j.

Lemma Ah_cont: forall dx vh,
  0 < dx ->
  vh 0%Z = 0 -> vh (ni dx) = 0 ->
  finite_dot_product_dx dx (Ah vh dx) vh
      <= 4*c*c/(dx*dx)*finite_squared_norm_dx dx vh.

Lemma Ah_pos: forall dx vh,
  0 < dx -> vh 0%Z = 0 -> vh (ni dx) = 0 ->
  0 <= finite_dot_product_dx dx (Ah vh dx) vh.

Hypothesis eps_lt: 0 < eps < 1.

Definition CFL (dX:R*R) := c * (snd dX) / (fst dX).

Definition CFL_condition_r (dX:R*R) := CFL dX <= 1 - eps.

Definition ni_exact (dX:R*R) := (IZR (ni (fst dX)) = (xmax-xmin) / (fst dX))%R.

Notation Oupse := (OuP (fun dX => 0 < fst dX /\ 0 < snd dX /\ ni_exact dX)).

Notation OupsCFLe := (OuP (fun dX => 0 < fst dX /\ 0 < snd dX /\ CFL_condition_r dX /\ ni_exact dX)).

Definition scheme_is_consistent_of_order (p q : nat) :=
    Oupse (fun (t: {t| Under T t}) dX => finite_squared_norm_dx (fst dX)
                          (fun j => truncation_error dX j (kt (projT1 t) (snd dX))))
      (fun dX => Rsqr ((fst dX) ^ p + (snd dX) ^ q)).

Definition scheme_is_convergent_of_order (p q : nat) :=
    OupsCFLe (fun (t: {t| Under T t}) dX => finite_norm_dx (fst dX)
                          (fun j => convergence_error dX j (kt (projT1 t) (snd dX))))
      (fun dX => (fst dX) ^ p + (snd dX) ^ q).

End Wave.

Section zeroness.

Variable xmin xmax:R.
Hypothesis xminLtxmax: xmin < xmax.
Variable c:R.
Variable dX:R*R.
Variable u0h u1h: R*R -> Z-> R.
Variable fh: R*R -> Z -> nat -> R.
Variable unh : R * R -> Z -> nat -> R.

Hypothesis unh_is_discrete_sol: is_discrete_solution xmin xmax c u0h u1h fh unh dX.

Notation ni := (fun x => (ni xmin xmax x)).

Lemma zeroness_unh_1: forall (n:nat), 0 < fst dX -> snd dX <> 0 ->
    unh dX 0 n = 0.

Lemma zeroness_unh_2: forall (n:nat), 0 < fst dX -> snd dX <> 0 ->
    unh dX (ni (fst dX)) n = 0.

End zeroness.

Section zeroness2.

Variable xmin xmax:R.
Hypothesis xminLtxmax: xmin < xmax.

Variable partial_derive_firstvar: (R * R -> R) -> (R * R -> R).
Variable partial_derive_secondvar: (R * R -> R) -> (R * R -> R).

Variable c:R.
Variable u0 u1: R->R.
Variable f : R*R -> R.
Variable dX:R*R.
Variable usol : R * R -> R.

Hypothesis usol_is_sol: is_solution xmin xmax partial_derive_firstvar partial_derive_secondvar c u0 u1 f usol.

Notation ni x := (ni xmin xmax x).

Lemma zeroness_usol_xmin_ni : forall t, fst dX <> 0 -> ni_exact xmin xmax dX ->
 usol (xmin + IZR (ni (fst dX)) * fst dX, t) =0.

Lemma zeroness_ubar_1: forall (n:nat),
    ubar xmin usol dX 0 n = 0.

Lemma zeroness_ubar_2: forall (n:nat),
    fst dX <> 0 -> ni_exact xmin xmax dX ->
    ubar xmin usol dX (ni (fst dX)) n = 0.

End zeroness2.

Section my_unh.

Variable xmin xmax c:R.
Variable u0h u1h: R*R -> Z-> R.
Variable fh: R*R -> Z -> nat -> R.

Fixpoint my_unh (dX:R*R) (j:Z) (n:nat) {struct n} : R :=
   if Z_lt_le_dec 0 j then
   if Z_lt_le_dec j (ni xmin xmax (fst dX)) then
    match n with
     | 0%nat => u0h dX j
     | S (0%nat as n1) => my_unh dX j n1 - snd dX*(snd dX/2*Ah c (fun j'=> my_unh dX j' n1) (fst dX)j - u1h dX j)
     | S (S n2 as n1) => (2*my_unh dX j n1 - my_unh dX j n2)
         - snd dX*snd dX*(Ah c (fun j'=> my_unh dX j' n1) (fst dX) j -fh dX j n1)
     end
  else 0%R
  else 0%R.

Lemma my_unh_correct: forall dX, is_discrete_solution xmin xmax c u0h u1h fh my_unh dX.

Lemma scheme_eq: forall unh1 unh2 nmax dX, snd dX <> 0 ->
    (forall j n, (0 <= j <= ni xmin xmax (fst dX))%Z -> (n <= nmax)%nat ->
         scheme xmin xmax c u0h u1h fh unh1 dX j n = scheme xmin xmax c u0h u1h fh unh2 dX j n)
 -> (forall j n, (0 <= j <= ni xmin xmax (fst dX))%Z -> (n <= nmax)%nat -> unh1 dX j n = unh2 dX j n).

End my_unh.