theory Example_van_der_Pol
imports
  "../Numerics/Euler_Affine"
  "../Numerics/Optimize_Float"
begin

declare [[ cd_patterns "(_ has_derivative ?f x) (at x within _)"]]
lemmas [simp] = split_beta' inverse_eq_divide

subsection {* Van der Pol oscillator *}
text {*\label{sec:vanderpol}*}

abbreviation "vanderpol_real \<equiv> \<lambda>(x::real, y::real). (y, y * (1 + - x*x) + - x)"
approximate_affine vanderpol vanderpol_real

schematic_lemma vanderpol_fderiv:
  "((\<lambda>(x::real, y::real). (y, y * (1 + - x*x) + - x))
    has_derivative (case x of (x, y) \<Rightarrow> \<lambda>(dx, dy). ?D x y dx dy))
   (at x within X)"
  by (auto intro!: derivative_eq_intros)

concrete_definition vanderpol_fderiv uses vanderpol_fderiv

lemma vanderpol_fderiv':
  "((\<lambda>(x::real, y::real). (y, y * (1 + - x*x) + - x))
    has_derivative
   (case x of (x, y) \<Rightarrow> \<lambda>(dx, dy). (dy, - (y * (2 * (dx * x))) + dy * (1 + - (x * x)) + - dx)))
   (at x within X)"
  by (auto intro!: derivative_eq_intros simp: split_beta inverse_eq_divide)

abbreviation "vanderpol_d_real \<equiv> \<lambda>(x::real, y) (dx, dy). (dy, - (y * (2 * (dx * x))) + dy * (1 + - (x * x)) + - dx)"
approximate_affine vanderpol_d vanderpol_d_real

abbreviation "vanderpol_dd_real \<equiv>
  \<lambda>(x::real, y::real). \<lambda>(dx, dy). \<lambda>(dx', dy'). (0 * dx', - (dy * (2 * (dx' * x))) + -
    (y * (2 * (dx * dx')) + dy' * (2 * (dx * x))))"
approximate_affine vanderpol_dd vanderpol_dd_real

lemma vanderpol_fderiv2: "((\<lambda>x. vanderpol_d_real x dx) has_derivative vanderpol_dd_real x dx) (at x)"
  by (auto intro!: derivative_eq_intros simp: split_beta')

abbreviation "vanderpol_ivp \<equiv> \<lambda>optns args. uncurry_options vanderpol optns (hd args) (tl args)"
abbreviation "vanderpol_d_ivp \<equiv> \<lambda>optns args. uncurry_options vanderpol_d optns (hd args) (hd (tl args)) (tl (tl args))"
abbreviation "vanderpol_dd_ivp \<equiv> \<lambda>optns args. uncurry_options vanderpol_dd optns (hd args) (hd (tl args)) (hd (tl (tl args))) (tl (tl (tl args)))"

interpretation vanderpol!: aform_approximate_ivp
  vanderpol_ivp
  vanderpol_d_ivp
  vanderpol_dd_ivp
  vanderpol_real
  vanderpol_d_real
  vanderpol_dd_real
  apply unfold_locales
  apply (rule Joints2_JointsI)
  unfolding list.sel
  apply (rule vanderpol, assumption, assumption)
  apply (drule length_set_of_apprs, simp)
  apply (rule vanderpol_fderiv')
  apply (rule vanderpol_fderiv2)
  apply (rule Joints2_JointsI)
  apply (rule vanderpol_d, assumption, assumption)
  apply (drule length_set_of_apprs, simp)
  apply (rule Joints2_JointsI)
  apply (rule vanderpol_dd, assumption, assumption)
  apply (drule length_set_of_apprs, simp)
  apply (auto intro!: continuous_intros simp: split_beta')
  done

definition "constrain x y = (if x = y then () else ())"
lemma [code]: "constrain x y = ()" by (simp add: constrain_def)

definition "decimal x ex = (if ex \<ge> 0 then x * 10 ^ nat ex else real_divl 53 x (10^nat (-ex)))"

definition "vanderpol_optns =
    \<lparr>
    precision = 53,
    tolerance = FloatR 1 (- 4),
    stepsize  = FloatR 1 (- 8),
    iterations = 40,
    halve_stepsizes = 10,
    widening_mod = 40,
    rk2_param = FloatR 1 (- 1),
    max_tdev_thres = FloatR 1 (- 3),
    presplit_summary_tolerance = FloatR 1 (- 2),
    collect_granularity = FloatR 1 (- 3),
    start_section = Sctn (0, -1) (- FloatR 100 6),
    checkpoint_sections = [],
    stop_sections = [Sctn (0, -1) (- decimal 225 (-2))],
    safe_appr = \<lambda>_. True,
    snap_sections = 1,
    min_section_distance = 1/2/2/2,
    next_section_distance_factor = 25,
    next_section_turn_distance_factor = 0,
    printing_fun = (let re = print_rectangle 0 0 0; af = print_aform; _ = constrain re af;
      both = (\<lambda>c. let _ = re c in af c) in both),
    irect_mod = 8,
    result_fun = ivls_result 20 40,
    max_intersection_step = FloatR 1 (- 10),
    adaptive_atol = FloatR 1 (-30),
    method_id = 2
  \<rparr>"

(* width(splitting) Euler, atol, reduction *)
definition vanderpol_optns'::"int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> (real \<times> real, (real \<times> real) aform, _) options"
  where "vanderpol_optns' w m a r = vanderpol_optns \<lparr>
      collect_granularity := FloatR 1 w,
      max_tdev_thres := FloatR 1 w,
      adaptive_atol := FloatR 1 a,
      start_section := if r\<noteq>0 then Sctn (0, -1) (- FloatR 100 (- 6)) else
        Sctn (0, -1) (- FloatR 1 100),
      next_section_distance_factor := if r ~= 0 then 100 else FloatR 1 100,
      next_section_turn_distance_factor := if r ~= 0 then 0 else FloatR 1 100,
      max_intersection_step := FloatR 1 (- 2),
      method_id := nat m\<rparr>"

definition
  vanderpol_limit::"[int, int, int, int, int] \<Rightarrow>_" where
  "vanderpol_limit w e a r err = poincare_distance_d_dres (vanderpol_ivp) (vanderpol_d_ivp) (vanderpol_dd_ivp)
  (vanderpol_optns' w e a r)
  [msum_aform' (aform_of_point (decimal 125 (-2), decimal 225 (-2))) (aform_of_ivl 0 (decimal err (-2), 0))]"

subsection {* Extraction of Code *}
text {* Code is extracted and can be used when imported using vanderpol.ML *}
export_code vanderpol_limit nat_of_integer int_of_integer in SML module_name Vanderpol file "vanderpol_core.ML"

subsection {* Evaluation *}

text {* The output can also be checked inside Isabelle/jEdit: open the "Raw Output" panel during
  the evaluation: *}

definition "vanderpol_limit' = do
  {
    res' \<leftarrow> vanderpol_limit (-2) 2 (-11) 1 50;
    let res = (List.maps fst res');
    dRETURN (map Inf_aform res, map Sup_aform res, map snd res')
  }"

value [code] vanderpol_limit'

lemma  vdpl_res: "vanderpol_limit' = dRETURN
  ([(FloatR 24759730456532176 (- 54), FloatR 20266198323167228 (- 53))],
   [(FloatR 26264698040279640 (- 54), FloatR 20266198323167236 (- 53))],
   [Sctn (FloatR 0 0, FloatR (- 1) 0) (FloatR (- 9) (- 2))])"
  by eval

context approximate_ivp begin

lemma poincare_distance_d_dRETURN:
  assumes [refine_vcg]: "0 < stepsize optns" "0 < rk2_param optns" "rk2_param optns \<le> 1"
  assumes start_below: "\<forall>x\<in>set X0S. local.appr_below_sctn optns x (start_section optns)"
  assumes start_safe: "\<forall>x\<in>set X0S. safe_appr optns x"
  assumes "poincare_distance_d_dres optns X0S = dRETURN stops"
  shows
      "(\<exists>CXS. Ball CXS (safe_appr optns) \<and> flowsto
      (UNION (set X0S) set_of_appr)
      (UNION (CXS) set_of_appr)
      (\<Union>(XS, sctn)\<in>set stops. (UNION (set XS) set_of_appr) \<inter> plane_of sctn))"
proof -
  note poincare_distance_d_dres.refine[of optns X0S]
  also note poincare_distance_d[OF assms(1-5)]
  finally show ?thesis
    unfolding assms by simp
qed

end

lemma
  vanderpol_invar:
  shows  "\<exists>CXS E. vanderpol.flowsto {(1.25, 2.25) .. (1.75, 2.25)} CXS E \<and> E \<subseteq> {(1.25, 2.25) .. (1.75, 2.25)}"
proof -
  from vdpl_res[unfolded vanderpol_limit'_def]
  obtain res fres where res: "vanderpol_limit (- 2) 2 (- 11) 1 50 = dRETURN [res]"
    "snd res = Sctn (0, -1) ((FloatR (- 9) (- 2)))"
    "fst res = [fres]"
    "Inf_aform fres = (FloatR 24759730456532176 (- 54), FloatR 20266198323167228 (- 53))"
      (is "_ = (?ix, ?iy)")
    "Sup_aform fres = (FloatR 26264698040279640 (- 54), FloatR 20266198323167236 (- 53))"
    (is "_ = ((?sx, ?sy))")
    by (cases "vanderpol_limit (- 2) 2 (- 11) 1 50") (auto simp: Let_def List.maps_def)
  have decs: "decimal 50 (- 2) = FloatR 1 (-1)"
    "decimal 125 (- 2) = FloatR 5 (- 2)"
    "decimal 175 (- 2) = FloatR 7 (- 2)"
    "decimal 225 (- 2) = FloatR 9 (- 2)"
    by eval+
  have dec_rat:
  "decimal 50 (- 2) = 0.5"
    "decimal 125 (- 2) = 1.25"
    "decimal 175 (- 2) = 1.75"
    "decimal 225 (- 2) = 2.25"
    by (auto simp: powr_minus powr_numeral decs)
  have "\<exists>CXS. Ball CXS (safe_appr (vanderpol_optns' (- 2) 2 (- 11) 1)) \<and>
      approximate_ivp.flowsto vanderpol_real
       (UNION
         (set [msum_aform' (aform_of_point (decimal 125 (- 2), decimal 225 (- 2)))
                (aform_of_ivl 0 (decimal (real (50::int)) (- 2), 0))])
         Affine)
       (UNION CXS Affine)
       (\<Union>(XS, sctn)\<in>set [res].
           UNION (set XS) Affine \<inter> plane_of sctn)"
    apply (rule vanderpol.poincare_distance_d_dRETURN[OF _ _ _ _ _ res(1)[unfolded vanderpol_limit_def]])
    apply eval+
    done
  then obtain CXS where
    "vanderpol.flowsto
      (Affine (msum_aform' (aform_of_point (1.25, 2.25))
        (aform_of_ivl 0 (decimal (real (50::int)) (- 2), 0))))
      (UNION CXS Affine)
      (Affine fres \<inter> plane_of (snd res))"
    by (auto simp: res dec_rat)
  also have "Affine (msum_aform' (aform_of_point (1.25, 2.25))
      (aform_of_ivl 0 (decimal (real (50::int)) (- 2), 0))) = {(1.25, 2.25) .. (1.75, 2.25)}"
    apply (auto simp: msum_appr_eq aform_of_point_def Affine_aform_of_ivl zero_prod_def
      powr_minus powr_numeral dec_rat)
    apply (rule_tac x = "5/4" in exI)
    apply (rule_tac x = "a - 5/4" in exI)
    apply auto
    done
  finally have 1: "vanderpol.flowsto
    {(1.25, 2.25) .. (1.75, 2.25)}
    (UNION CXS Affine)
    (Affine fres \<inter> plane_of (snd res))"
    .
  have "Affine fres \<inter> plane_of (snd res) \<subseteq> {Inf_aform fres .. Sup_aform fres} \<inter> plane_of (snd res)"
    by auto
  also have "\<dots> = {(?ix, 2.25) .. (?sx, 2.25)}" (is "_ = ?end")
    unfolding res
    by (auto simp add: powr_minus powr_numeral plane_of_def )
  finally have "Affine fres \<inter> plane_of (snd res) \<subseteq> ?end" .
  from vanderpol.flowsto_subset[OF 1 this order_refl]
  have "vanderpol.flowsto {(1.25, 2.25) .. (1.75, 2.25)} (UNION CXS Affine) ?end"
    "?end \<subseteq> {(1.25, 2.25) .. (1.75, 2.25)}"
    by (auto simp: plane_of_def powr_minus powr_numeral)
  thus ?thesis
    by blast
qed

lemma vanderpol_flow_invar:
  "\<exists>E. E \<subseteq> {(1.25, 2.25) .. (1.75, 2.25)} \<and>
    (\<forall>x\<in>{(1.25, 2.25) .. (1.75, 2.25)}.
      \<exists>t\<ge>0. t \<in> vanderpol.existence_ivl x \<and>
            vanderpol.flow x t \<in> E)"
  using vanderpol_invar
  unfolding vanderpol.flowsto_def
  apply safe
  apply (rule_tac x = E in exI)
  by auto

end
