Theory Abstract_Reachability_Analysis
theory Abstract_Reachability_Analysis
imports
Abstract_Rigorous_Numerics
Affine_Arithmetic.Affine_Arithmetic
"../Refinement/Refine_String"
"../Refinement/Refine_Folds"
Ordinary_Differential_Equations.Flow
Runge_Kutta
begin
subsection ‹Misc›
lemma nth_concat_exists:
"∃k j. concat xs ! i = xs ! k ! j ∧ k < length xs ∧ j < length (xs ! k)"
if "i < length (concat xs)"
using that
proof (induction xs arbitrary: i)
case Nil
then show ?case by auto
next
case (Cons xs xss)
from Cons.prems consider "i < length xs"
| "i ≥ length xs" "i < length xs + length (concat xss)"
by (cases "i < length xs") auto
then show ?case
proof cases
case 1
then show ?thesis
by (force simp: nth_append intro: exI[where x=i] exI[where x=0])
next
case 2
then have "i - length xs < length (concat xss)" by arith
with Cons.IH[of "i - length xs"]
obtain k j where
"concat xss ! (i - length xs) = xss ! k ! j" "k < length xss" "j < length (xss ! k)"
by auto
then show ?thesis
using 2
by (fastforce simp: nth_append nth_Cons split: nat.splits
intro: exI[where x=j] exI[where x="k + 1"])
qed
qed
lemma nth_concatE:
assumes "i < length (concat xs)"
obtains k j where "concat xs ! i = xs ! k ! j" "k < length xs" "j < length (xs ! k)"
apply atomize_elim
using assms nth_concat_exists by blast
lemma max_Var_floatariths_concat:
"max_Var_floatariths (concat xs) ≤ k"
if "⋀x. x ∈ set xs ⟹ max_Var_floatariths x ≤ k"
using that max_Var_floatarith_le_max_Var_floatariths_nthI
by (fastforce simp: in_set_conv_nth intro!: max_Var_floatariths_leI
elim!: nth_concatE)
lemma max_Var_floatariths_list_update:
"max_Var_floatariths (xs[xa := y]) ≤ k"
if "max_Var_floatariths (xs) ≤ k"
and "max_Var_floatarith y ≤ k"
by (metis neq_le_trans linorder_le_cases list_update_beyond
max_Var_floatariths_list_updateI that)
lemma max_Var_floatarith_0[simp]: "max_Var_floatarith 0 = 0"
and max_Var_floatarith_1[simp]: "max_Var_floatarith 1 = 0"
by (auto simp: zero_floatarith_def one_floatarith_def)
lemma list_set_rel_br: "⟨Id⟩list_set_rel = br set distinct"
by (auto simp: list_set_rel_def)
lemma
br_list_relD:
shows "(x, y) ∈ ⟨br a i⟩list_set_rel ⟹ y = a ` set x ∧ list_all i x"
apply (auto simp: list_set_rel_def br_def list_rel_def)
subgoal premises prems for s t
using prems
by (induction arbitrary: y rule: list.rel_induct) auto
subgoal premises prems for s t
using prems
by (induction arbitrary: y rule: list.rel_induct) auto
subgoal premises prems for s
using prems
by (induction arbitrary: y rule: list.rel_induct) auto
done
lemma sctn_rel_br: "⟨br a I⟩sctn_rel = br (λx. case x of Sctn n p ⇒ Sctn (a n) p) (λx. I (normal x))"
apply (auto simp: sctn_rel_def br_def in_rel_def[abs_def] split: sctn.splits)
subgoal for b x1 x2 by (cases b) auto
subgoal for a b by (cases a; cases b) auto
done
lemma br_list_rel: "⟨br a I⟩list_rel = br (map a) (list_all I)"
by (fastforce simp: list_rel_def br_def list_all2_iff Ball_def in_set_zip list_all_length
intro!: nth_equalityI)
lemma list_set_rel_brp: "⟨br a I⟩list_set_rel = br (λxs. a ` set xs) (λxs. list_all I xs ∧ distinct (map a xs))"
unfolding list_set_rel_def br_list_rel br_chain o_def o_def
by (auto)
declare INF_cong_simp [cong] SUP_cong_simp [cong] image_cong_simp [cong del]
context auto_ll_on_open begin
definition "stable_on CX trap ⟷
(∀t x0. flow0 x0 t ∈ trap ⟶ t ∈ existence_ivl0 x0 ⟶ t > 0 ⟶
(∀s∈{0<..t}. flow0 x0 s ∈ CX) ⟶ x0 ∈ trap)"
lemma stable_onD:
"⋀t x0. flow0 x0 t ∈ trap ⟹ t ∈ existence_ivl0 x0 ⟹ t > 0 ⟹
(⋀s. 0 < s ⟹ s ≤ t ⟹ flow0 x0 s ∈ CX) ⟹
x0 ∈ trap"
if "stable_on CX trap"
using that by (auto simp: stable_on_def)
lemma nonneg_interval_mem_existence_ivlI[intro]:
"0 ≤ t1 ⟹ t1 ≤ t2 ⟹ t2 ∈ existence_ivl0 x0 ⟹ {t1..t2} ⊆ existence_ivl0 x0"
"t1 ≤ t2 ⟹ t2 ≤ 0 ⟹ t1 ∈ existence_ivl0 x0 ⟹ {t1..t2} ⊆ existence_ivl0 x0"
"t1 ≤ 0 ⟹ 0 ≤ t2 ⟹ t1 ∈ existence_ivl0 x0 ⟹ t2 ∈ existence_ivl0 x0 ⟹ {t1..t2} ⊆ existence_ivl0 x0"
apply auto
apply (drule ivl_subset_existence_ivl) apply auto
apply (drule ivl_subset_existence_ivl') apply auto
apply (drule segment_subset_existence_ivl, assumption)
apply (auto simp: closed_segment_eq_real_ivl)
done
lemma interval_subset_existence_ivl:
"t ∈ existence_ivl0 x0 ⟹ s ∈ existence_ivl0 x0 ⟹ t ≤ s ⟹ {t .. s} ⊆ existence_ivl0 x0"
using segment_subset_existence_ivl[of s x0 t]
by (auto simp: closed_segment_eq_real_ivl)
end
lemma(in c1_on_open_euclidean) diff_existence_ivl_iff[simp]:
"t2 - t1 ∈ existence_ivl0 (flow0 x0 t1) ⟷ t2 ∈ existence_ivl0 x0"
if "t1 ≤ t2" "t1 ∈ existence_ivl0 x0"
apply auto
apply (drule existence_ivl_trans[OF that(2)])
apply (auto intro!: diff_existence_ivl_trans that)
done
lemma (in auto_ll_on_open) flow_trans':
"flow0 (flow0 x0 t1) t2 = flow0 x0 (t1 + t2)"
if "t1 ∈ existence_ivl0 x0" "t1 + t2 ∈ existence_ivl0 x0"
apply (subst flow_trans)
using that
by (auto intro!: existence_ivl_trans')
context auto_ll_on_open begin
definition "flowpipe0 X0 hl hu CX X1 ⟷ 0 ≤ hl ∧ hl ≤ hu ∧ X0 ⊆ X ∧ CX ⊆ X ∧ X1 ⊆ X ∧
(∀(x0) ∈ X0. ∀h ∈ {hl .. hu}. h ∈ existence_ivl0 x0 ∧ (flow0 x0 h) ∈ X1 ∧ (∀h' ∈ {0 .. h}. (flow0 x0 h') ∈ CX))"
lemma flowpipe0D:
assumes "flowpipe0 X0 hl hu CX X1"
shows flowpipe0_safeD: "X0 ∪ CX ∪ X1 ⊆ X"
and flowpipe0_nonneg: "0 ≤ hl" "hl ≤ hu"
and flowpipe0_exivl: "hl ≤ h ⟹ h ≤ hu ⟹ (x0) ∈ X0 ⟹ h ∈ existence_ivl0 x0"
and flowpipe0_discrete: "hl ≤ h ⟹ h ≤ hu ⟹ (x0) ∈ X0 ⟹ (flow0 x0 h) ∈ X1"
and flowpipe0_cont: "hl ≤ h ⟹ h ≤ hu ⟹ (x0) ∈ X0 ⟹ 0 ≤ h' ⟹ h' ≤ h ⟹ (flow0 x0 h') ∈ CX"
using assms
by (auto simp: flowpipe0_def)
lemma flowpipe0_source_subset: "flowpipe0 X0 hl hu CX X1 ⟹ X0 ⊆ CX"
apply (auto dest: bspec[where x=hl] bspec[where x=0] simp: flowpipe0_def)
apply (drule bspec)
apply (assumption)
apply (drule bspec[where x=hl])
apply auto
apply (drule bspec[where x=0])
by (auto simp: flow_initial_time_if)
end
subsection ‹Options›
definition [refine_vcg_def]: "precision_spec = SPEC (λprec::nat. True)"
definition [refine_vcg_def]: "adaptive_atol_spec = SPEC (λx::real. True)"
definition [refine_vcg_def]: "adaptive_rtol_spec = SPEC (λx::real. True)"
definition [refine_vcg_def]: "method_spec = SPEC (λm::nat. True)"
definition [refine_vcg_def]: "start_stepsize_spec = SPEC (λx::real. x > 0)"
definition [refine_vcg_def]: "iterations_spec = SPEC (λn::nat. True)"
definition [refine_vcg_def]: "halve_stepsizes_spec = SPEC (λn::nat. True)"
definition [refine_vcg_def]: "widening_mod_spec = SPEC (λn::nat. True)"
definition [refine_vcg_def]: "rk2_param_spec = SPEC (λr::real. 0 < r ∧ r ≤ 1)"
typedef ode_ops = "{(ode_e::floatarith list, safe_form::form).
open_form safe_form ∧
max_Var_floatariths ode_e ≤ length ode_e ∧
max_Var_form safe_form ≤ length ode_e}"
by (auto intro!: exI[where x="[floatarith.Num 0]"]
exI[where x="Less (floatarith.Num 0) (floatarith.Num 1)"])
setup_lifting type_definition_ode_ops
lift_definition ode_expression::"ode_ops ⇒ floatarith list" is fst .
lift_definition safe_form_expr::"ode_ops ⇒ form" is snd .
lemma open_form_ode_op[intro, simp]: "open_form (safe_form_expr odo)"
and max_Var_ode_expression: "max_Var_floatariths (ode_expression odo) ≤ length (ode_expression odo)"
and max_Var_form_safe_form_expr: "max_Var_form (safe_form_expr odo) ≤ length (ode_expression odo)"
by (transfer, auto)+
lift_definition (code_dt) mk_ode_ops::"floatarith list ⇒ form ⇒ ode_ops option" is
"λode_e safe_form.
if (open_form safe_form ∧ max_Var_floatariths ode_e ≤ length ode_e ∧ max_Var_form safe_form ≤ length ode_e)
then Some (ode_e, safe_form) else None"
by (auto simp:)
lemma
assumes "mk_ode_ops e s = Some odo"
shows ode_expression_mk_ode_ops: "ode_expression odo = e"
and safe_form_expr_mk_ode_ops: "safe_form_expr odo = s"
using assms
by (transfer, simp split: if_splits prod.splits)+
locale ode_operations = fixes ode_ops::ode_ops begin
definition "ode_e = ode_expression ode_ops"
definition "safe_form = safe_form_expr ode_ops"
definition ode::"'a ⇒ 'a::executable_euclidean_space"
where "ode x = eucl_of_list (interpret_floatariths ode_e (list_of_eucl x))"
definition "ode_d_expr_nth N n i =
FDERIV_floatarith
(FDERIV_n_floatarith (ode_e ! i) [0..<N] (map floatarith.Var [N..<2 * N]) n) [0..<N]
(map floatarith.Var [2 * N..<3 * N])"
definition "ode_d_expr N n =
(FDERIV_floatariths
(FDERIV_n_floatariths ode_e [0..<N] (map floatarith.Var [N..<2 * N]) n)
[0..<N]
(map floatarith.Var [2 * N..< 3 * N]))"
definition ode_d_raw::"nat ⇒ 'a ⇒ 'a ⇒ 'a ⇒ 'a::executable_euclidean_space"
where "ode_d_raw n x dn d =
eucl_of_list (interpret_floatariths (ode_d_expr DIM('a) n) (list_of_eucl x @ list_of_eucl dn @ list_of_eucl d))"
definition "ode_fa_nth xs i = subst_floatarith (λi. xs ! i) (ode_e ! i)"
definition "ode_fa xs = map (subst_floatarith (λi. xs ! i)) ode_e"
definition "ode_d_fa_nth n xs ds i = subst_floatarith (λi. (xs@ds@ds) ! i) (ode_d_expr_nth (length xs) n i)"
definition "ode_d_fa n xs ds = map (subst_floatarith (λi. (xs@ds@ds) ! i)) (ode_d_expr (length xs) n)"
definition safe::"'a::executable_euclidean_space ⇒ bool"
where "safe x ⟷
length ode_e = DIM('a) ∧
max_Var_floatariths ode_e ≤ DIM('a) ∧
open_form safe_form ∧
max_Var_form safe_form ≤ DIM('a) ∧
interpret_form safe_form (list_of_eucl x) ∧
isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x)"
definition "Csafe = Collect safe"
definition "euler_incr_fas_nth X0 h CX i = X0 ! i + h * (ode_fa_nth CX i)"
definition "euler_incr_fas X0 h CX = map (euler_incr_fas_nth X0 h CX) [0..<length X0]"
definition "euler_err_fas_nth X0 h CX i = ((h ^⇩e 2) / floatarith.Num 2) * ode_d_fa_nth 0 CX (ode_fa CX) i"
definition "euler_err_fas X0 h CX = map (euler_err_fas_nth X0 h CX) [0..<length X0]"
definition "euler_fas X0 h CX =
map (λi. (euler_incr_fas_nth X0 h X0 i + euler_err_fas_nth X0 h CX i)) [0..<length X0] @
euler_err_fas X0 h CX"
definition "rk2_fas_err_nth rkp x0 h cx s2 i =
((((h ^⇩e 3 / 6) *
(ode_d_fa_nth 1 cx (ode_fa cx) i +
ode_d_fa_nth 0 cx (ode_d_fa 0 cx (ode_fa cx)) i)))
- ((h ^⇩e 3 * rkp / 4) *
ode_d_fa_nth 1 (euler_incr_fas x0 (s2 * h * rkp) x0) (ode_fa x0) i))"
definition "rk2_fas_err rkp x0 h cx s2 = map (rk2_fas_err_nth rkp x0 h cx s2) [0..<length x0]"
definition "rk2_fas rkp x0 h cx s2 =
(map (λi.
((x0 ! i +
h * ((1 - (1 / (rkp * 2))) * ode_fa_nth x0 i +
(1 / (rkp * 2)) * ode_fa_nth (euler_incr_fas x0 (h * rkp) x0) i))
+ rk2_fas_err_nth rkp x0 h cx s2 i)) [0..<length x0]) @ rk2_fas_err rkp x0 h cx s2"
lemma ode_d_expr_nth: "i < length ode_e ⟹ ode_d_expr_nth N n i = ode_d_expr N n ! i "
by (auto simp: ode_d_expr_nth_def ode_d_expr_def
FDERIV_n_floatariths_nth)
lemma length_ode_d_expr[simp]: "length (ode_d_expr f n) = length ode_e"
by (induction n) (auto simp: ode_d_expr_def FDERIV_n_floatariths_def)
lemma ode_fa_nth: "i < length ode_e ⟹ ode_fa xs ! i = ode_fa_nth xs i"
by (auto simp: ode_fa_nth_def ode_fa_def)
lemma ode_d_fa_nth: "i < length ode_e ⟹ ode_d_fa_nth n xs ds i = ode_d_fa n xs ds ! i"
by (auto simp: ode_d_fa_def ode_d_fa_nth_def ode_d_expr_nth)
lemma length_ode_d_fa[simp]: "length (ode_d_fa n xs ds) = length ode_e"
by (auto simp: ode_d_fa_def FDERIV_n_floatariths_def)
lemma length_rk2_fas_err[simp]: "length (rk2_fas_err rkp x0 h cx s2) = length x0"
by (simp add: rk2_fas_err_def)
lemma length_euler_incr_fas[simp]: "length (euler_incr_fas X0 h CX) = length X0"
by (auto simp: euler_incr_fas_def)
lemma length_euler_err_fas[simp]: "length (euler_err_fas X0 h CX) = length X0"
by (auto simp: euler_err_fas_def)
lemma length_euler_floatarith[simp]: "length (euler_fas X0 h CX) = 2 * length X0"
by (auto simp: euler_fas_def)
lemma length_rk2_fas[simp]: "length (rk2_fas rkp x0 h cx s2) = 2 * length x0"
by (simp add: rk2_fas_def)
lemma open_safe: "open Csafe"
proof -
have leq: "list_updates [0..<DIM('a)] (list_of_eucl x) (replicate DIM('a) 0) = list_of_eucl x" for x::'a
by (auto intro!: nth_equalityI simp: list_updates_nth)
have "(Csafe::'a set) =
(if length ode_e = DIM('a) ∧ max_Var_floatariths ode_e ≤ DIM('a) ∧ max_Var_form safe_form ≤ DIM('a) ∧ open_form safe_form then
{x. interpret_form safe_form (list_of_eucl x)} ∩
{x. isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x)}
else {})"
by (auto simp: Csafe_def safe_def)
also have "open …"
apply (auto intro!: open_Int)
subgoal premises prems using open_form[OF prems(4), where 'a='a, of "[0..<DIM('a)]" "replicate (DIM('a)) 0"]
by (auto simp: leq)
subgoal
apply (rule isFDERIV_open)
apply (rule order_trans)
apply assumption
apply arith
done
done
finally show ?thesis .
qed
lemma safeD:
fixes x::"'a::executable_euclidean_space"
assumes "x ∈ Csafe"
shows "interpret_form safe_form (list_of_eucl x)"
and safe_isFDERIV: "isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x)"
using assms
by (auto simp: Csafe_def safe_def)
lemma
fixes x::"'a::executable_euclidean_space"
shows safe_max_Var: "x ∈ Csafe ⟹ max_Var_floatariths ode_e ≤ DIM('a)"
and safe_length: "x ∈ Csafe ⟹ length ode_e = DIM('a)"
and safe_max_Var_form: "x ∈ Csafe ⟹ max_Var_form safe_form ≤ DIM('a)"
by (auto simp: safe_def Csafe_def)
lemma safe_isFDERIV_append:
fixes x::"'a::executable_euclidean_space"
shows "x ∈ Csafe ⟹ isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x @ xs)"
apply (rule isFDERIV_max_Var_congI)
apply (rule safe_isFDERIV)
apply assumption
using safe_max_Var[of x]
by (auto simp: nth_append)
lemma ode_d_raw_0:
assumes "x ∈ Csafe"
shows "(ode has_derivative ode_d_raw 0 x d) (at x)"
using assms safe_max_Var[OF assms] safe_length[OF assms]
unfolding ode_def ode_d_raw_def ode_d_expr_def
apply (intro interpret_floatarith_FDERIV_floatariths[THEN has_derivative_eq_rhs])
apply (auto simp: isFDERIV_def FDERIV_n_floatariths_def safe_max_Var nth_append
max_Var_floatariths_Max Csafe_def safe_def
intro!: arg_cong[where f=eucl_of_list] ext interpret_floatariths_FDERIV_floatariths_cong
freshs_floatariths_max_Var_floatarithsI
max_Var_floatarith_le_max_Var_floatariths[le])
apply (rule interpret_floatariths_max_Var_cong)
apply (auto simp: max_Var_floatariths_Max Max_gr_iff nth_append
dest!: less_le_trans[OF _ max_Var_floatarith_DERIV_floatarith])
apply (drule max_Var_floatariths_lessI) apply simp
apply (auto dest!: less_le_trans[OF _ safe_max_Var[OF assms]])
apply (drule max_Var_floatariths_lessI) apply simp
apply (auto dest!: less_le_trans[OF _ safe_max_Var[OF assms]])
done
lemma not_fresh_odeD: "x ∈ Csafe ⟹ ¬fresh_floatariths ode_e i ⟹ i < DIM('a)" for x::"'a::executable_euclidean_space"
using fresh_floatariths_max_Var[of ode_e i] safe_max_Var[of x] by arith
lemma safe_isnFDERIV:
fixes x::"'a::executable_euclidean_space"
assumes "x ∈ Csafe"
shows "isnFDERIV DIM('a) ode_e [0..<DIM('a)] [DIM('a)..<2 * DIM('a)] (list_of_eucl x @ ys) n"
apply (rule isFDERIV_imp_isnFDERIV)
apply (rule isFDERIV_max_Var_congI)
apply (rule safe_isFDERIV[OF assms])
using safe_max_Var[OF assms] safe_length[OF assms]
by (auto simp: nth_append)
lemma safe_isnFDERIVI:
assumes "(eucl_of_list xs::'a::executable_euclidean_space) ∈ Csafe"
assumes [simp]: "length xs = DIM('a)" "length ds = DIM('a)"
shows "isnFDERIV DIM('a) ode_e [0..<DIM('a)] [DIM('a)..<2 * DIM('a)] (xs@ds) n"
proof -
have "isnFDERIV DIM('a) ode_e [0..<DIM('a)] [DIM('a)..<2 * DIM('a)] (list_of_eucl (eucl_of_list xs::'a)@ds) n"
by (rule safe_isnFDERIV; fact)
also
have "list_of_eucl (eucl_of_list xs::'a) = xs"
by (auto intro!: nth_equalityI)
finally show ?thesis .
qed
lemma dest_Num_eq_Some_iff[simp]: "dest_Num_fa fa = (Some x) ⟷ fa = floatarith.Num x"
by (cases fa) auto
lemma ode_d_raw_Suc:
includes floatarith_notation
assumes "x ∈ Csafe"
shows "((λx. ode_d_raw n x d d) has_derivative ode_d_raw (Suc n) x d) (at x)"
proof -
let ?shift = "λx. floatarith.Var (if 2 * DIM('a) ≤ x ∧ x < 3 * DIM('a) then x - DIM('a) else x)"
have subst_ode_e[simp]: "map (subst_floatarith ?shift) ode_e = ode_e"
apply (auto intro!: nth_equalityI)
apply (rule subst_floatarith_Var_max_Var_floatarith)
by (auto dest!: max_Var_floatariths_lessI
less_le_trans[OF _ safe_max_Var[OF assms]])
have map_shift[simp]:
"(map ?shift [DIM('a)..<2 * DIM('a)]) = (map floatarith.Var [DIM('a)..<2 * DIM('a)])"
"(map ?shift [2 * DIM('a)..<3 * DIM('a)]) =
(map floatarith.Var [DIM('a)..<2 * DIM('a)])"
by (auto intro!: nth_equalityI)
show ?thesis
unfolding ode_def ode_d_raw_def ode_d_expr_def
apply (rule interpret_floatarith_FDERIV_floatariths_append[THEN has_derivative_eq_rhs])
subgoal
proof -
let ?shift = "λx. if 2 * DIM('a) ≤ x ∧ x < 3 * DIM('a) then x - DIM('a) else x"
have mv: "max_Var_floatariths
(FDERIV_floatariths (FDERIV_n_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]) n)
[0..<DIM('a)] (map floatarith.Var [2 * DIM('a)..<3 * DIM('a)])) ≤ 3 * DIM('a)"
and mv2: "max_Var_floatariths
(FDERIV_floatariths (FDERIV_n_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]) n)
[0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)])) ≤ 2 * DIM('a)"
by (auto intro!:
max_Var_floatarith_FDERIV_floatariths[le]
max_Var_floatarith_FDERIV_n_floatariths[le]
safe_max_Var[OF assms, le])
have eq: "(map (subst_floatarith (λi. floatarith.Var (if 2 * DIM('a) ≤ i ∧ i < 3 * DIM('a) then i - DIM('a) else i)))
((FDERIV_floatariths (FDERIV_n_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]) n)
[0..<DIM('a)] (map floatarith.Var [2 * DIM('a)..<3 * DIM('a)])))) =
(FDERIV_floatariths (FDERIV_n_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]) n)
[0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]))"
apply (rule nth_equalityI)
apply auto defer
apply (subst subst_floatarith_Var_FDERIV_floatarith[where 'a='a], force, force, force)
apply (subst subst_floatarith_Var_FDERIV_n_nth[where 'a='a], force, force, force, force)
by (simp add: o_def)
show ?thesis
apply (subst isFDERIV_subst_Var_floatarith[symmetric, where s="?shift"])
subgoal by (auto intro!: mv[le] max_Var_floatariths_fold_const_fa[le])
subgoal by (auto simp: nth_append)
subgoal by (auto intro!: mv[le])
subgoal
proof -
have "isnFDERIV DIM('a) ode_e [0..<DIM('a)] [DIM('a)..<2*DIM('a)] (list_of_eucl x @ list_of_eucl d) (Suc (Suc n))"
apply (rule safe_isnFDERIVI)
using assms
by auto
from this[simplified, THEN conjunct1]
show ?thesis
unfolding eq isnFDERIV.simps
apply (rule isFDERIV_max_Var_congI)
apply (frule less_le_trans[OF _ mv2])
apply (auto simp: nth_append)
done
qed
done
qed
subgoal
by (auto intro!: safe_max_Var[OF assms, le]
max_Var_floatarith_FDERIV_floatariths[le]
max_Var_floatarith_FDERIV_n_floatariths[le])
subgoal using safe_length assms by simp
subgoal
apply (auto simp add: nth_append
intro!: ext arg_cong[where f=eucl_of_list] interpret_floatariths_FDERIV_floatariths_cong
freshs_floatariths_max_Var_floatarithsI
safe_max_Var[OF assms, le]
max_Var_floatarith_FDERIV_floatariths[le]
max_Var_floatarith_FDERIV_n_floatariths[le])
apply (rule nth_equalityI)
apply auto
subgoal premises prems for h i j
proof -
have *: "(list_of_eucl x @ list_of_eucl d @ list_of_eucl d @ list_of_eucl h) =
(map (λi. interpret_floatarith (?shift i)
(list_of_eucl x @ list_of_eucl d @ list_of_eucl d @ list_of_eucl h)) [0..<4 * DIM('a)])"
by (auto intro!: nth_equalityI simp: nth_append)
have mv_le: "max_Var_floatarith
(DERIV_floatarith j
(FDERIV_floatarith
(FDERIV_n_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]) n ! i)
[0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]))) ≤
2 * DIM('a)"
"max_Var_floatarith
(DERIV_floatarith j
(FDERIV_floatarith (FDERIV_n_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [DIM('a)..<2 * DIM('a)]) n ! i)
[0..<DIM('a)] (map floatarith.Var [2 * DIM('a)..<3 * DIM('a)])))
≤ 3 * DIM('a)"
by (auto intro!: prems
safe_max_Var[OF assms, le]
max_Var_floatarith_le_max_Var_floatariths_nth[le]
max_Var_floatarith_DERIV_floatarith[le]
max_Var_floatarith_FDERIV_floatarith[le]
max_Var_floatarith_FDERIV_n_floatariths[le])
show ?thesis
apply (subst *)
apply (subst interpret_floatarith_subst_floatarith[symmetric])
apply (auto intro!: prems mv_le[le])
apply (subst subst_floatarith_Var_DERIV_floatarith, use prems in force)
apply (subst subst_floatarith_Var_FDERIV_floatarith[where 'a='a], force, force, force)
apply (subst subst_floatarith_Var_FDERIV_n_nth[where 'a='a], force, force, force, use prems in force)
apply (auto simp: o_def prems nth_append intro!: interpret_floatarith_max_Var_cong
dest!: less_le_trans[OF _ mv_le(1)])
done
qed
done
done
qed
lift_definition ode_d::"nat ⇒ 'a::executable_euclidean_space ⇒ 'a ⇒ 'a ⇒⇩L 'a" is
"λn x dn d. if x ∈ Csafe
then ode_d_raw n x dn d
else 0"
subgoal for n x dn
apply (cases n)
subgoal
by (cases "x ∈ Csafe")
(auto intro: has_derivative_bounded_linear[OF ode_d_raw_0])
subgoal for n'
apply (cases "x ∈ Csafe")
subgoal
apply (simp del: isnFDERIV.simps)
apply (rule has_derivative_bounded_linear)
apply (rule ode_d_raw_Suc)
apply assumption
done
subgoal by (simp del: isnFDERIV.simps)
done
done
done
definition "ode_d1 x = ode_d 0 x 0"
lemma ode_has_derivative:
assumes "isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x)"
assumes "(x::'a::executable_euclidean_space) ∈ Csafe"
shows "(ode has_derivative ode_d1 x) (at x)"
proof -
from assms(1) have *: "x ∈ Csafe ⟹ isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x @ list_of_eucl (0::'a))"
apply (rule isFDERIV_max_Var_congI)
using safe_max_Var[of x]
by (auto simp: nth_append)
show ?thesis
unfolding ode_d1_def
apply (transfer fixing: x)
apply (rule ode_d_raw_0[THEN has_derivative_eq_rhs])
by (auto intro!: * assms)
qed
lemma ode_has_derivative_safeI:
assumes "x ∈ Csafe"
shows "(ode has_derivative ode_d1 x) (at x)"
using assms
by (auto simp: safe_def Csafe_def intro!: ode_has_derivative)
lemma ode_d1_eq: "ode_d1 x = ode_d 0 x j"
unfolding ode_d1_def
proof (transfer fixing: x j, rule ext, cases "x ∈ Csafe", clarsimp_all, goal_cases)
case (1 d)
have "isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x @ list_of_eucl (0::'a)) =
isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x @ list_of_eucl j)"
by (rule isFDERIV_max_Var_cong)
(auto dest!: less_le_trans[OF _ safe_max_Var[OF 1]] simp: nth_append)
moreover
have "interpret_floatariths (FDERIV_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [2 * DIM('a)..<3 * DIM('a)]))
(list_of_eucl x @ list_of_eucl (0::'a) @ list_of_eucl d) =
interpret_floatariths (FDERIV_floatariths ode_e [0..<DIM('a)] (map floatarith.Var [2 * DIM('a)..<3 * DIM('a)]))
(list_of_eucl x @ list_of_eucl j @ list_of_eucl d)"
using 1
by (intro interpret_floatariths_fresh_cong)
(auto dest!: not_fresh_FDERIV_floatariths not_fresh_odeD
simp: nth_append)
ultimately show ?case
by (auto simp: ode_d_raw_def ode_d_expr_def)
qed
lemma eventually_Collect_open:
assumes "P x" "open (Collect P)"
shows "eventually P (at x)"
using assms(1) assms(2) eventually_at_topological by blast
lemma ode_d_has_derivative:
assumes "x ∈ Csafe"
shows "((λx. ode_d n x d d) has_derivative ode_d (Suc n) x d) (at x)"
apply (transfer fixing: n d x)
using assms
apply (simp del: isnFDERIV.simps)
apply (rule if_eventually_has_derivative)
subgoal by (rule ode_d_raw_Suc)
subgoal
by (rule eventually_Collect_open)
(auto simp: safe_max_Var[OF assms] open_safe intro!: safe_max_Var[OF assms, le])
subgoal by (simp add: isnFDERIV.simps)
subgoal by simp
done
lemma ode_d1_has_derivative:
assumes "x ∈ Csafe"
shows "(ode_d1 has_derivative ode_d (Suc 0) x) (at x)"
proof (rule blinfun_has_derivative_componentwiseI[THEN has_derivative_eq_rhs])
fix i::'a assume "i ∈ Basis"
show "((λx. blinfun_apply (ode_d1 x) i) has_derivative ode_d (Suc 0) x i) (at x)"
unfolding ode_d1_eq[of _ i]
apply (rule ode_d_has_derivative)
apply fact
done
next
show "(λxa. ∑i∈Basis. blinfun_scaleR (blinfun_inner_left i) (blinfun_apply (ode_d (Suc 0) x i) xa)) = ode_d (Suc 0) x"
apply (rule ext)
apply (auto intro!: ext euclidean_eqI[where 'a='a] blinfun_euclidean_eqI
simp: blinfun.bilinear_simps inner_sum_left inner_Basis if_distrib if_distribR
sum.delta' cong: if_cong)
apply (rule arg_cong[where f="λx. x ∙ b" for b])
proof goal_cases
case (1 j i b)
from eventually_isFDERIV[where params=Nil, simplified, OF safe_isFDERIV[OF assms] order_trans[OF safe_max_Var[of x]]]
have "∀⇩F x in at x. isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl x)"
by (auto simp: assms)
then obtain S where S: "x ∈ S" "open S" "S ⊆ Csafe"
and "⋀xa. xa ∈ S ⟹ xa ≠ x ⟹ isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl xa)"
using assms open_safe safe_isFDERIV by auto
then have S_FDERIV: "⋀s. s ∈ S ⟹
isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl s)"
using safe_isFDERIV[OF assms]
by auto
interpret second_derivative_on_open "ode" ode_d1 "ode_d (Suc 0) x" x S
proof standard
fix a assume "a ∈ S"
with S have "a ∈ Csafe" by auto
from S_FDERIV[OF ‹a ∈ S›]
have "isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl a)" by simp
then have "isFDERIV DIM('a) [0..<DIM('a)] ode_e (list_of_eucl a)"
apply (rule isFDERIV_max_Var_congI)
using safe_max_Var[of x]
by (auto simp: nth_append)
then show "(ode has_derivative blinfun_apply (ode_d1 a)) (at a)"
using ‹a ∈ Csafe›
by (rule ode_has_derivative)
next
fix i
interpret linear "ode_d (Suc 0) x"
proof
fix y z
have 1: "((λx. ode_d 0 x (y + z) (y + z)) has_derivative ode_d (Suc 0) x (y + z)) (at x)"
apply (rule ode_d_has_derivative)
apply (rule assms)
done
have *: "ode_d 0 x (y + z) (y + z) = ode_d 0 x y y + ode_d 0 x z z" for x
by (auto simp: blinfun.bilinear_simps ode_d1_eq[symmetric])
have 2: "((λx. ode_d 0 x (y + z) (y + z)) has_derivative
ode_d (Suc 0) x y + ode_d (Suc 0) x z) (at x)"
apply (subst *)
apply (rule derivative_eq_intros)
apply (rule ode_d_has_derivative)
apply fact
apply (rule ode_d_has_derivative)
apply fact
apply (auto simp: blinfun.bilinear_simps)
done
from has_derivative_unique[OF 1 2]
show "ode_d (Suc 0) x (y + z) = ode_d (Suc 0) x y + ode_d (Suc 0) x z"
by (auto intro!: blinfun_eqI)
next
fix r y
have 1: "((λx. ode_d 0 x (r *⇩R y) (r *⇩R y)) has_derivative ode_d (Suc 0) x (r *⇩R y)) (at x)"
by (rule ode_d_has_derivative; fact)
have *: "ode_d 0 x (r *⇩R y) (r *⇩R y) = r *⇩R ode_d 0 x y y" for x
by (auto simp: blinfun.bilinear_simps ode_d1_eq[symmetric])
have 2: "((λx. ode_d 0 x (r *⇩R y) (r *⇩R y)) has_derivative
r *⇩R ode_d (Suc 0) x y) (at x)"
apply (subst *)
apply (rule derivative_eq_intros)
apply (rule ode_d_has_derivative; fact)
apply (auto simp: blinfun.bilinear_simps)
done
from has_derivative_unique[OF 1 2]
show "(ode_d (Suc 0) x (r *⇩R y)) = (r *⇩R ode_d (Suc 0) x y)"
by (auto intro!: blinfun_eqI)
qed
show "((λx. blinfun_apply (ode_d1 x) i) has_derivative blinfun_apply (ode_d (Suc 0) x i))
(at x)"
apply (subst euclidean_representation[of i, symmetric])
apply (subst (2) euclidean_representation[of i, symmetric])
apply (auto simp: blinfun.bilinear_simps)
apply (rule derivative_eq_intros)
apply (rule derivative_eq_intros)
apply (subst_tac j = i in ode_d1_eq)
apply (rule ode_d_has_derivative)
apply (rule assms)
apply force
apply (auto simp: blinfun.bilinear_simps[symmetric]
intro!: ext euclidean_eqI[where 'a='a] blinfun_euclidean_eqI)
apply (rule arg_cong[where f="λx. x ∙ b" for b])
by (auto simp: sum scaleR)
next
show "x ∈ S" "open S" by fact+
qed
show ?case
by (rule symmetric_second_derivative) fact
qed
qed
lemma ode_d1_has_derivative_safeI:
assumes "x ∈ Csafe"
shows "(ode_d1 has_derivative ode_d (Suc 0) x) (at x)"
apply (rule ode_d1_has_derivative)
using assms by (auto simp: safe_def)
sublocale c1_on_open_euclidean ode ode_d1 Csafe
by unfold_locales
(auto simp: continuous_on_eq_continuous_within at_within_open[OF _ open_safe]
intro!: derivative_eq_intros continuous_at_imp_continuous_on open_safe
ode_has_derivative_safeI continuous_blinfun_componentwiseI
has_derivative_continuous ode_d1_has_derivative_safeI)
definition ivlflows ::
"'a::executable_euclidean_space sctn set
⇒ (('a × 'a ⇒⇩L 'a) set
⇒ ('a × 'a ⇒⇩L 'a) set × ('a × 'a ⇒⇩L 'a) set)
⇒ ('a × 'a ⇒⇩L 'a) set ⇒ 'a sctn ⇒ bool"
where "ivlflows stops stopcont trap rsctn =
(∀ivl. ivl ⊆ ⋃(plane_of ` stops) × UNIV ⟶
ivl ⊆ (snd (stopcont ivl)) ∧
fst (stopcont ivl) ⊆ snd (stopcont ivl) ∧
(fst (stopcont ivl)) ⊆ sbelow_halfspace rsctn × UNIV ∧
(snd (stopcont ivl)) ⊆ sbelow_halfspace rsctn × UNIV ∧
flowsto (ivl) {0..} ((snd (stopcont ivl))) ((fst (stopcont ivl)) ∪ trap))"
lift_definition ode_d2::"'a::executable_euclidean_space ⇒ 'a ⇒⇩L 'a ⇒⇩L 'a" is "λx.
if x ∈ Csafe then ode_d 1 x else (λ_. 0)"
by (auto intro!: has_derivative_bounded_linear ode_d1_has_derivative)
definition ode_na::"real × _ ⇒ _" where "ode_na = (λa. ode (snd a))"
definition ode_d_na::"real × _ ⇒ (real × _) ⇒⇩L _" where "ode_d_na = (λtx. ode_d1 (snd tx) o⇩L snd_blinfun)"
definition ode_d2_na::"real × _ ⇒ (real × _) ⇒⇩L (real × _) ⇒⇩L _" where
"ode_d2_na = (λtx. flip_blinfun (flip_blinfun (ode_d2 (snd tx) o⇩L snd_blinfun) o⇩L snd_blinfun))"
definition "euler_incr_fas' D = (map fold_const_fa (euler_incr_fas (map floatarith.Var [0..<D]) (floatarith.Var (D))
(map floatarith.Var [Suc D..<Suc (2*D)])))"
definition "euler_fas' D = (map fold_const_fa (euler_fas (map floatarith.Var [0..<D])
(floatarith.Var (2*D)) (map floatarith.Var [D..<2*D])))"
definition "rk2_fas' D = (map fold_const_fa (rk2_fas
(floatarith.Var (2*D))
(map floatarith.Var [0..<D])
(floatarith.Var (2*D+1))
(map floatarith.Var [D..<2*D])
(floatarith.Var (2*D+2))))"
lemma [autoref_rules]: "(euler_incr_fas', euler_incr_fas') ∈ nat_rel → fas_rel"
"(euler_fas', euler_fas') ∈ nat_rel → fas_rel"
"(rk2_fas', rk2_fas') ∈ nat_rel → fas_rel"
by auto
definition "solve_poincare_fas n =
(let D = length ode_e in
map floatarith.Var [0..<D] @ concat (map (λi . map (λj .
(if i ≠ n then floatarith.Var (D + i * D + j) - (floatarith.Var(D + n * D + j) * (ode_e ! i) / (ode_e ! n))
else 0)
) [0..<D]) [0..<D]))"
end
definition "nonempty X ⟷ X ≠ {}"
definition pad_zeroes :: "nat ⇒ real list set ⇒ real list set"
where [simp]: "pad_zeroes n X = (λxs. xs @ replicate n (0::real)) ` X"
locale approximate_sets_ode = approximate_sets where ops = ops + ode_operations
where ode_ops = ode_ops
for ops:: "'b approximate_set_ops"
and ode_ops::"ode_ops"
begin
definition "D = (length ode_e)"
definition "ode_slp = slp_of_fas ode_e"
definition "euler_slp = slp_of_fas (euler_fas' D)"
definition "euler_incr_slp = slp_of_fas (euler_incr_fas' D)"
definition "rk2_slp = slp_of_fas (rk2_fas' D)"
definition "solve_poincare_slp = map (λi. slp_of_fas (map fold_const_fa (solve_poincare_fas i))) [0..<D]"
definition safe_set
where "safe_set (X::'a::executable_euclidean_space set) = do {
b1 ← approx_form_spec safe_form (list_of_eucl ` X);
b2 ← isFDERIV_spec D [0..<D] ode_e (list_of_eucl ` X);
RETURN (b1 ∧ b2)
}"
definition "wd TYPE('a::executable_euclidean_space) ⟷ length ode_e = DIM('a)"
lemma open_safe_form[intro, simp]: "open_form safe_form"
by (auto simp: safe_form_def)
lemma max_Var_floatariths_ode_e_le: "max_Var_floatariths ode_e ≤ D"
and max_Var_form_safe_form_le: "max_Var_form safe_form ≤ D"
using max_Var_ode_expression[of ode_ops] max_Var_form_safe_form_expr[of ode_ops]
by (auto simp: ode_e_def safe_form_def D_def)
lemma wdD:
assumes "wd TYPE('a::executable_euclidean_space)"
shows "length ode_e = DIM('a)" "max_Var_floatariths ode_e ≤ DIM('a)"
"max_Var_form safe_form ≤ DIM('a)"
"ode_e ≠ []" "D = DIM('a)"
using assms max_Var_floatariths_ode_e_le max_Var_form_safe_form_le
by (auto simp: wd_def D_def safe_form_def ode_e_def)
definition "mk_safe (X::'a::executable_euclidean_space set) = do {
ASSERT (wd TYPE('a));
s ← safe_set (X:::appr_rel::'a set);
if s then RETURN (X:::appr_rel) else SUCCEED
}"
definition
"mk_safe_coll X = do {
XS ← (sets_of_coll X);
FORWEAK XS (RETURN op_empty_coll)
(λx. do {
s ← mk_safe (x);
RETURN (mk_coll s)
})
(λb c. RETURN (b ∪ c))
}"
definition ode_set::"'a::executable_euclidean_space set ⇒ 'a set nres" where "ode_set X = do {
_ ← mk_safe X;
approx_slp_appr ode_e ode_slp (list_of_eucl ` (X))
}"
definition
"Picard_step X0 t0 h X = SPEC (λR.
case R of
Some R ⇒
nonempty R ∧ compact R ∧ (R ⊆ Csafe) ∧
(∀x0 ∈ X0. ∀h'∈{t0 .. t0 + h}. ∀phi∈cfuncset t0 h' X.
x0 + integral {t0 .. h'} (λt. ode (phi t)) ∈ R)
| None ⇒ True)"
lemmas [refine_vcg_def] = approx_form_spec_def isFDERIV_spec_def
lemma safe_set_spec[THEN order.trans, refine_vcg]:
assumes "wd TYPE('a::executable_euclidean_space)"
shows "safe_set X ≤ SPEC (λr. r ⟶ (X::'a set) ⊆ Csafe)"
unfolding safe_set_def
by (refine_vcg) (auto simp del: isnFDERIV.simps simp add: Csafe_def safe_def replicate_eq_list_of_eucl_zero wdD[OF ‹wd _›])
definition Picard_step_ivl :: "'a::executable_euclidean_space set ⇒ real ⇒ real ⇒ 'a set ⇒ 'a set option nres" where
"Picard_step_ivl X0 t0 h X = do {
ASSERT (0 ≤ h);
ASSERT (wd TYPE('a));
let H = lv_ivl [0] [h];
let D = DIM('a);
let env = concat ` listset [list_of_eucl ` X0, H, list_of_eucl ` X];
env ← approx_slp_spec (euler_incr_fas' D) D euler_incr_slp env;
(case env of
Some env ⇒
do {
(l, u) ← op_ivl_rep_of_set ((eucl_of_list ` env::'a set));
ASSERT (l ≤ u);
s ← safe_set {l .. u};
if s then do {
r ← mk_safe ({l .. u}:::appr_rel);
RETURN (Some (r:::appr_rel))
} else RETURN None
}
| None ⇒ RETURN None)
}"
definition "do_widening_spec (i::nat) = SPEC (λb::bool. True)"
primrec P_iter::"'a::executable_euclidean_space set ⇒ real ⇒ nat ⇒ ('a) set ⇒ ('a) set option nres" where
"P_iter X0 h 0 X = do {
let _ = trace_set (ST ''P_iter failed (0)'') (Some (X));
RETURN None
}"
| "P_iter X0 h (Suc i) X = do {
ASSERT (0 ≤ h);
(l, u) ← op_ivl_rep_of_set (X);
ASSERT (l ≤ u);
ivl ← mk_safe ({l .. u}:::appr_rel);
X' ← Picard_step_ivl X0 0 h ivl;
(case X' of
Some X' ⇒ do {
(l', u') ← op_ivl_rep_of_set (X');
do_widening ← do_widening_spec i;
let l' = inf l' l - (if do_widening then abs (l' - l) else 0);
let u' = sup u' u + (if do_widening then abs (u' - u) else 0);
ASSERT (l' ≤ u');
s ← safe_set {l' .. u'};
if s
then do {
ivl' ← mk_safe {l' .. u'};
if (l ≤ l' ∧ u' ≤ u) then RETURN (Some ivl)
else P_iter X0 h i ivl'
} else do {
let _ = trace_set (ST ''P_iter failed (unsafe interval after widening)'') (Some (X));
RETURN None
}
}
| None ⇒ do {
let _ = trace_set (ST ''P_iter failed (Picard_step)'') (Some (X));
RETURN None
}
)
}"
context fixes m::"('a::executable_euclidean_space set ⇒ real ⇒ real ⇒ 'a set ⇒ ('a set × 'c) option nres)"
begin
primrec cert_stepsize::
"'a set ⇒ real ⇒ nat ⇒ nat ⇒ (real × 'a set × 'a set × 'c) nres"
where
"cert_stepsize X0 h n 0 = do { let _ = trace_set (ST ''cert_stepsize failed'') (Some (X0)); SUCCEED}"
| "cert_stepsize X0 h n (Suc i) = do {
(l, u) ← op_ivl_rep_of_set (X0);
ASSERT (0 ≤ h);
ASSERT (l ≤ u);
ivl ← mk_safe {l .. u};
ASSERT (ivl ≠ {});
X' ← P_iter X0 h n ivl;
case X' of Some X' ⇒
do {
r1 ← m X0 h h X';
r2 ← m X0 0 h X';
(case (r1, r2) of
(Some (res, err), Some (res_ivl, _)) ⇒
do {
_ ← mk_safe res;
_ ← mk_safe res_ivl;
RETURN (h, res, res_ivl, err)
}
| _ ⇒
do {
let _ = trace_set (ST ''cert_stepsize method failed'') (Some (X'));
cert_stepsize X0 (h / 2) n i
}
)
}
| None ⇒ cert_stepsize X0 (h / 2) n i
}"
end
definition "one_step_method m ⟷ (∀X0 CX hl hu. m X0 hl hu CX ≤
SPEC (λr. case r of None ⇒ True | Some (res, err) ⇒ nonempty res ∧
(∀x0 ∈ X0. ∀h ∈ {hl .. hu}. x0 ∈ Csafe ⟶ h ≥ 0 ⟶ h ∈ existence_ivl0 x0 ⟶
(∀h' ∈ {0 .. h}. flow0 x0 h' ∈ CX) ⟶ x0 + h *⇩R ode x0 ∈ CX ⟶ flow0 x0 h ∈ res)))"
definition "one_step X0 h m = do {
CHECKs ''one_step nonneg'' (0 < h);
its ← iterations_spec;
halvs ← halve_stepsizes_spec;
(h, res, res_ivl, err) ← cert_stepsize m X0 h its halvs;
ASSERT (0 < h);
RETURN (h, err, res_ivl, res)
}"
definition [refine_vcg_def]: "default_reduce_argument_spec = SPEC (λx::unit. True)"
definition "euler_step X0 h = one_step X0 h (λX0 hl hu CX.
do {
let H = lv_ivl [min hl hu] [max hl hu];
_ ← mk_safe CX;
let env = concat ` listset [list_of_eucl ` X0, list_of_eucl ` CX, H];
env ← approx_slp_spec (euler_fas' DIM('a)) (2 * DIM('a)) euler_slp env;
case env of None ⇒ RETURN None
| Some env ⇒ do {
let res' = take DIM('a) ` env;
ASSERT (env_len res' DIM('a));
let res = (eucl_of_list ` res');
ASSUME (ncc res);
let err' = drop DIM('a) ` take (DIM('a) * 2) ` env;
ASSERT (env_len err' DIM('a));
let err = (eucl_of_list ` err'::'a::executable_euclidean_space set);
ra ← default_reduce_argument_spec;
res ← reduce_spec ra res;
ASSUME (ncc res);
s ← safe_set res;
if s then
do {
res ← mk_safe res;
RETURN (Some (res::'a set, err))
} else RETURN None
}
})"
definition "rk2_step X0 h = one_step X0 h (λX0 hl hu CX.
do {
let H = lv_ivl [min hl hu] [max hl hu];
rps ← rk2_param_spec;
let rkp = lv_ivl [rps] [rps];
let s2 = lv_ivl [0] [1];
_ ← mk_safe CX;
ASSUME (ncc CX);
let env = concat ` listset [list_of_eucl ` X0, list_of_eucl ` CX, rkp, H, s2];
env ← approx_slp_spec (rk2_fas' DIM('a)) (2 * DIM('a)) rk2_slp env;
case env of None ⇒ RETURN None
| Some env ⇒ do {
let res' = take DIM('a) ` env;
ASSERT (env_len res' DIM('a));
let res = (eucl_of_list ` res'::'a::executable_euclidean_space set);
ASSUME (ncc res);
let err' = drop DIM('a) ` take (DIM('a) * 2) ` env;
ASSERT (env_len err' DIM('a));
let err = (eucl_of_list ` err'::'a set);
ra ← default_reduce_argument_spec;
res ← reduce_spec ra res;
ASSUME (ncc res);
s ← safe_set res;
if s then
do {
res ← mk_safe res;
RETURN (Some (res, err))
} else RETURN None
}
})"
definition "choose_step X0 h = do {
mid ← method_spec;
(if mid = 2 then rk2_step X0 h else euler_step X0 h)
}"
definition "ode_e' = (ode_e @
mmult_fa D D D (concat (map (λj. map (λi.
(FDERIV_floatarith (ode_e ! j) [0..<D] ((replicate D 0)[i := 1]))) [0..<D]) [0..<D]))
(map floatarith.Var [D..<D + D*D]))"
definition "transversal_directions f =
do {
(I, S) ← op_ivl_rep_of_set f;
RETURN (sum_list (map (λb. (if I ∙ b ≤ 0 then if S ∙ b ≤ 0 then S ∙ b else 0 else if S ∙ b ≥ 0 then I ∙ b else 0) *⇩R b)
(Basis_list::'a::executable_euclidean_space list)))
}"
definition "intersects_sctns X' sctns = do {
ASSUME (finite sctns);
FOREACH⇗λsctns' b. ¬b ⟶ X' ∩ ⋃(plane_of ` (sctns - sctns')) = {}⇖ sctns
(λsctn b. do {b' ← op_intersects ( X') sctn; RETURN (b ∨ b')}) False
}"
definition "trace_sets s X = do {
XS ← sets_of_coll (X:::clw_rel (appr_rel)); FORWEAK XS (RETURN ()) (λX. RETURN (trace_set s (Some X))) (λ_ _. RETURN ())
}"
definition "print_sets s X = do {
XS ← sets_of_coll (X:::clw_rel (appr_rel)); FORWEAK XS (RETURN ()) (λX. RETURN (print_set s (X))) (λ_ _. RETURN ())
}"
definition "intersects_sctns_spec_clw R sctns = do {
Rs ← sets_of_coll ((R:::clw_rel appr_rel):::clw_rel(appr_rel));
FORWEAK Rs (RETURN False) (λR. intersects_sctns R sctns) (λa b. RETURN (a ∨ b))
}"
definition [simp]: "nonneg_reals = ({0..}::real set)"
definition [simp]: "pos_reals = ({0<..}::real set)"
definition "nonzero_component s X n = do {
I ← Inf_inner X n;
S ← Sup_inner X n;
CHECKs s (I > 0 ∨ S < 0)
}"
definition "disjoints_spec X Y = do {
Xi ← ivls_of_sets X;
IS ← inter_overappr (Xi:::clw_rel lvivl_rel) (Y:::clw_rel lvivl_rel);
RETURN (is_empty IS)
}"
definition subset_spec_plane :: "'a::executable_euclidean_space set ⇒ 'a sctn ⇒ bool nres" where
"subset_spec_plane X sctn = do {
CHECKs ''subset_spec_plane: not in Basis'' (abs (normal sctn) ∈ set Basis_list);
(i, s) ← ivl_rep X;
RETURN (i ∙ normal sctn = pstn sctn ∧ s ∙ normal sctn = pstn sctn)
}"
definition "op_eventually_within_sctn X sctn S = do {
(l, u) ← ivl_rep S;
(xl, xu) ← op_ivl_rep_of_set X;
CHECKs (ST ''op_eventually_within_sctn: empty ivl'') (l ≤ u);
CHECKs (ST ''op_eventually_within_sctn: not in Basis'') (abs (normal sctn) ∈ set Basis_list);
b ← subset_spec_plane S sctn;
CHECKs (ST ''op_eventually_within_sctn: subset_spec_plane 1'') b;
b ← subset_spec_plane ({xl .. xu}:::lvivl_rel) sctn;
CHECKs (ST ''op_eventually_within_sctn: subset_spec_plane 2'') b;
RETURN (b ∧ (∀i ∈ set Basis_list - {abs (normal sctn)}. l ∙ i < xl ∙ i ∧ xu ∙ i < u ∙ i))
}"
definition [simp]: "uninfo X = X"
definition [simp]: "op_subset_ivl a b ⟷ a ⊆ b"
definition [simp]: "op_eq_ivl a b ⟷ a = b"
abbreviation "iplane_rel ≡ λA. ⟨A, ⟨lv_rel⟩plane_rel⟩inter_rel"
abbreviation "isbelow_rel ≡ λA. ⟨A, ⟨lv_rel⟩sbelow_rel⟩inter_rel"
abbreviation "isbelows_rel ≡ λA. ⟨A, ⟨lv_rel⟩sbelows_rel⟩inter_rel"
definition [refine_vcg_def]: "get_plane X = SPEC (λsctn. X = plane_of sctn)"
definition "tolerate_error Y E =
do {
(ei, es) ← op_ivl_rep_of_set (E);
(yi, ys) ← op_ivl_rep_of_set (Y);
let ea = sup (abs ei) (abs es);
let ya = sup (abs yi) (abs ys);
rtol ← adaptive_rtol_spec;
atol ← adaptive_atol_spec;
let errtol = sup (rtol *⇩R ya) (atol *⇩R sum_list Basis_list);
RETURN (ea ≤ errtol, infnorm ea)
}"
definition "adapt_stepsize_fa rtol m e h' =
floatarith.Num (float_of h') *
floatarith.Powr (floatarith.Num (float_of (rtol)) / floatarith.Num (float_of e))
(inverse (floatarith.Num (float_of (real_of_nat m) + 1)))"
end
text ‹With ODE operations for variational equation›
locale approximate_sets_ode' = approximate_sets_ode
where ops = ops
and ode_ops = ode_ops
for ops :: "'b approximate_set_ops"
and ode_ops
begin
lift_definition var_ode_ops::ode_ops is "(ode_e', safe_form)"
using max_Var_floatariths_ode_e_le max_Var_form_safe_form_le
by (auto simp: ode_e'_def D_def length_concat o_def sum_list_triv
intro!: max_Var_floatariths_mmult_fa[le] max_Var_floatariths_concat max_Var_floatariths_mapI
max_Var_floatarith_FDERIV_floatarith[le] max_Var_floatariths_list_update
max_Var_floatariths_replicateI
max_Var_floatarith_le_max_Var_floatariths_nth[le])
sublocale var: approximate_sets_ode where ode_ops = var_ode_ops
by unfold_locales
end
lifting_update ode_ops.lifting
lifting_forget ode_ops.lifting
end