section{* Runge-Kutta methods*}
theory Runge_Kutta
imports
  One_Step_Method
  "~~/src/HOL/Library/Float"
  "../../Affine_Arithmetic/Executable_Euclidean_Space"
  "../Library/Multivariate_Taylor"
begin

lemma cball_trans: "y \<in> cball z b \<Longrightarrow> x \<in> cball y a \<Longrightarrow> x \<in> cball z (b + a)"
  apply (auto simp: dist_triangle3)
  by (smt dist_triangle)

text{*\label{sec:rk}*}

subsection {* Definitions *}
text{*\label{sec:rk-definition}*}

declare setsum.cong[fundef_cong]
fun rk_eval :: "(nat\<Rightarrow>nat\<Rightarrow>real) \<Rightarrow> (nat\<Rightarrow>real) \<Rightarrow> (real\<times>'a::real_vector \<Rightarrow> 'a) \<Rightarrow> real \<Rightarrow> real \<Rightarrow> 'a \<Rightarrow> nat \<Rightarrow> 'a" where
  "rk_eval A c f t h x j =
  f (t + h * c j, x + h *\<^sub>R (\<Sum>l=1 ..< j. A j l *\<^sub>R rk_eval A c f t h x l))"

primrec rk_eval_dynamic :: "(nat\<Rightarrow>nat\<Rightarrow>real) \<Rightarrow> (nat\<Rightarrow>real) \<Rightarrow> (real\<times>'a::{comm_monoid_add, scaleR} \<Rightarrow> 'a) \<Rightarrow> real \<Rightarrow> real \<Rightarrow> 'a \<Rightarrow> nat \<Rightarrow> (nat \<Rightarrow> 'a)" where
  "rk_eval_dynamic A c f t h x 0 = (\<lambda>i. 0)"
| "rk_eval_dynamic A c f t h x (Suc j) =
  (let K = rk_eval_dynamic A c f t h x j in
  K(Suc j:=f (t + h * c (Suc j), x + h *\<^sub>R (\<Sum>l=1..j. A (Suc j) l *\<^sub>R K l))))"

definition rk_increment where
  "rk_increment f s A b c h t x = (\<Sum>j=1..s. b j *\<^sub>R rk_eval A c f t h x j)"

definition rk_increment' where
  "rk_increment' error f s A b c h t x =
  eucl_down error (\<Sum>j=1..s. b j *\<^sub>R rk_eval A c f t h x j)"

definition euler_increment where
  "euler_increment f = rk_increment f 1 (\<lambda>i j. 0) (\<lambda>i. 1) (\<lambda>i. 0)"

definition euler where
  "euler f = grid_function (discrete_evolution (euler_increment f))"

definition euler_increment' where
"euler_increment' e f = rk_increment' e f 1 (\<lambda>i j. 0) (\<lambda>i. 1) (\<lambda>i. 0)"

definition euler' where
  "euler' e f = grid_function (discrete_evolution (euler_increment' e f))"

definition rk2_increment where
  "rk2_increment x f = rk_increment f 2 (\<lambda>i j. if i = 2 \<and> j = 1 then x else 0)
    (\<lambda>i. if i = 1 then 1 - 1 / (2 * x) else 1 / (2 * x)) (\<lambda>i. if i = 2 then x else 0)"

definition rk2 where
  "rk2 x f = grid_function (discrete_evolution (rk2_increment x f))"

subsection {* Euler method is consistent *}
text{*\label{sec:rk-euler-cons}*}

lemma euler_increment:
  fixes f::"_ \<Rightarrow> 'a::real_vector"
  shows "euler_increment f h t x = f (t, x)"
  unfolding euler_increment_def rk_increment_def
  by (subst rk_eval.simps) (simp del: rk_eval.simps)

lemma euler_float_increment:
  fixes f::"_ \<Rightarrow> 'a::executable_euclidean_space"
  shows "euler_increment' e f h t x = eucl_down e (f (t, x))"
  unfolding euler_increment'_def rk_increment'_def
  by (subst rk_eval.simps) (simp del: rk_eval.simps)

lemma euler_lipschitz:
  fixes x::"real \<Rightarrow> real"
  fixes f::"_ \<Rightarrow> 'a::real_normed_vector"
  assumes t: "t \<in> {t0..T}"
  assumes lipschitz: "\<forall>t\<in>{t0..T}. lipschitz D' (\<lambda>x. f (t, x)) L"
  shows "lipschitz D' (euler_increment f h t) L"
  using t lipschitz
  by (simp add: lipschitz_def euler_increment del: One_nat_def)

lemma rk2_increment:
  fixes f::"_ \<Rightarrow> 'a::real_vector"
  shows "rk2_increment p f h t x =
    (1 - 1 / (p * 2)) *\<^sub>R f (t, x) +
    (1 / (p * 2)) *\<^sub>R f (t + h * p, x + (h * p) *\<^sub>R f (t, x))"
  unfolding rk2_increment_def rk_increment_def
  apply (subst rk_eval.simps)
  apply (simp del: rk_eval.simps add: numeral_2_eq_2)
  apply (subst rk_eval.simps)
  apply (simp del: rk_eval.simps add:  field_simps)
  done

subsection {* Set-Based Consistency of Euler Method *}
text {*\label{sec:setconsistent}*}

locale derivative_set_bounded =
  derivative_on_convex +
  fixes F F'
  assumes f_set_bounded: "bounded F" "\<And>t x. t\<in>T \<Longrightarrow> x\<in>X \<Longrightarrow> (x, f (t, x)) \<in> F"
  assumes f'_ivl_bounded: "is_interval F'" "bounded F'" "\<And>t x d. t\<in>T \<Longrightarrow> (x, d) \<in> F \<Longrightarrow>
    f' (t,x) (1, d) \<in> F'"
begin

lemma F_nonempty: "F \<noteq> {}"
  and F'_nonempty: "F' \<noteq> {}"
  using nonempty
  unfolding ex_in_conv[symmetric]
  by (auto intro!: f_set_bounded f'_ivl_bounded)

lemma euler_consistent_traj_set:
  fixes t
  assumes ht: "0 \<le> h" "t + h \<le> u"
  assumes T: "{t..u} \<subseteq> T"
  assumes x': "\<And>s. s \<in> {t..u} \<Longrightarrow> (x has_vector_derivative f (s, x s)) (at s within {t..u})"
  assumes x: "\<And>s. s \<in> {t..u} \<Longrightarrow> x s \<in> X"
  shows "x (t + h) - discrete_evolution (euler_increment f) (t + h) t (x t) \<in> op *\<^sub>R (h\<^sup>2 / 2) ` F'"
proof cases
  assume "h = 0"
  from F'_nonempty obtain f' where "f' \<in> F'" by auto
  from this `h = 0` show ?thesis
    by (auto simp: discrete_evolution_def)
next
  assume "h \<noteq> 0"
  from this ht have "t < u" by simp
  from ht have line_subset: "(\<lambda>ta. t + ta * h) ` {0..1} \<subseteq> {t..u}"
    by (auto intro!: order_trans[OF add_left_mono[OF mult_left_le_one_le]])
  hence line_in: "\<And>s. 0 \<le> s \<Longrightarrow> s \<le> 1 \<Longrightarrow> t + s * h \<in> {t..u}"
    by (rule set_mp) auto
  from ht have subset: "{t .. t + h} \<subseteq> {t .. u}" by simp
  let ?T = "{t..u}"
  from ht have subset: "{t .. t + h} \<subseteq> {t .. u}" by simp
  from `t < u` have "t \<in> ?T" by auto
  from `t < u` have tx: "t \<in> T" "x t \<in> X" using assms by auto
  from tx assms have "0 \<le> norm (f (t, x t))" by simp
  have x_diff: "\<And>s. s \<in> ?T \<Longrightarrow> x differentiable at s within ?T"
    by (rule differentiableI, rule x'[simplified has_vector_derivative_def])
  have f': "\<And>t x. t \<in> ?T \<Longrightarrow> x \<in> X \<Longrightarrow> (f has_derivative f' (t, x)) (at (t, x) within (?T \<times> X))"
    using T by (intro has_derivative_subset[OF f']) auto
  let ?p = "(\<lambda>t. f' (t, x t) (1, f (t, x t)))"
  def diff \<equiv> "\<lambda>n::nat. if n = 0 then x else if n = 1 then \<lambda>t. f (t, x t) else ?p"
  have diff_0[simp]: "diff 0 = x" by (simp add: diff_def)
  {
    fix m::nat and ta::real
    assume mta: "m < 2" "t \<le> ta" "ta \<le> t + h"
    have image_subset: "(\<lambda>xa. (xa, x xa)) ` {t..u} \<subseteq> {t..u} \<times> X"
      using assms by auto
    note has_derivative_in_compose[where f="(\<lambda>xa. (xa, x xa))" and g = f, derivative_intros]
    note has_derivative_subset[OF _ image_subset, derivative_intros]
    note f'[derivative_intros]
    note x'[simplified has_vector_derivative_def, derivative_intros]
    have [simp]: "\<And>c x'. c *\<^sub>R f' (ta, x ta) x' = f' (ta, x ta) (c *\<^sub>R x')"
      using mta ht assms by (auto intro!: f' linear_cmul[symmetric] has_derivative_linear)
    have "((\<lambda>t. f (t, x t)) has_vector_derivative f' (ta, x ta) (1, f (ta, x ta))) (at ta within {t..u})"
      unfolding has_vector_derivative_def
      using assms ht mta by (auto intro!: derivative_eq_intros)
    hence "(diff m has_vector_derivative diff (Suc m) ta) (at ta within {t..t + h})"
      using mta ht
      by (auto simp: diff_def intro!: has_vector_derivative_within_subset[OF _ subset] x')
  } note diff = this
  from taylor_up_within_vector[of 2 t "t + h" diff x t, OF _ _ diff] ht `h\<noteq>0`
  obtain tt where tt: "tt \<in> Basis \<rightarrow> {t <..< t + h}"
    and taylor: "x (t + h) = x t + h *\<^sub>R f (t, x t) + (h\<^sup>2/2) *\<^sub>R (\<Sum>i\<in>Basis. ((?p (tt i)) \<bullet> i) *\<^sub>R i)"
    by (auto simp: Pi_iff Bex_def numeral_2_eq_2 diff_def real_of_nat_Suc
      scaleR_setsum_right ac_simps power2_eq_square inverse_eq_divide setprod_numseg)
  hence "x (t + h) - discrete_evolution (euler_increment f) (t + h) t (x t) =
    (h\<^sup>2 / 2) *\<^sub>R (\<Sum>i\<in>Basis. ((?p (tt i)) \<bullet> i) *\<^sub>R i)" (is "?d = _")
    by (auto simp: discrete_evolution_def euler_increment)
  also have "\<dots> \<in> op *\<^sub>R (h\<^sup>2 / 2) ` F'"
    apply (rule image_eqI[OF refl mem_box_componentwiseI[OF f'_ivl_bounded(1)]])
    using T tt ht
    apply (force simp: Pi_iff intro: f_set_bounded x image_eqI[OF refl f'_ivl_bounded(3)])
    done
  finally show ?thesis .
qed

lemma euler_consistent_traj_set2:
  fixes t
  assumes ht: "0 \<le> h" "t1 \<le> u"
  assumes T: "{t..u} \<subseteq> T"
  assumes x': "\<And>s. s \<in> {t..u} \<Longrightarrow> (x has_vector_derivative f (s, x s)) (at s within {t..u})"
  assumes x: "\<And>s. s \<in> {t..u} \<Longrightarrow> x s \<in> X"
  assumes *: "t1 = t + h"
  shows "x t1 - discrete_evolution (euler_increment f) t1 t (x t) \<in> op *\<^sub>R (h\<^sup>2 / 2) ` F'"
  using ht T x' x
  unfolding *
  by (rule euler_consistent_traj_set)

end

lemma bilinear_lsetsum: "bilinear f \<Longrightarrow> f (setsum g S) x = setsum ((\<lambda>i. f (g i) x)) S"
  apply (auto simp add: bilinear_def linear_setsum)
  apply (subst linear_setsum) back back back
  apply (auto simp: )
  done

lemma bilinear_rsetsum: "bilinear f \<Longrightarrow> finite S \<Longrightarrow> f x (setsum g S) = setsum ((\<lambda>i. (f x (g i)))) S"
  by (auto simp add: bilinear_def linear_setsum)

lemma mult_le_oneI: "0 \<le> (a::real) \<and> a \<le> 1 \<and> b \<le> 1 \<Longrightarrow> a * b \<le> 1"
  by (metis ODE_Auxiliarities.linordered_idom_class.mult_right_le_one_le order_trans)

text {* TODO: use @{term second_derivative} here *}
locale derivative2_set_bounded =
  fixes f::"real*'a \<Rightarrow> 'a::euclidean_space" and f' f''
  assumes f': "\<And>t x. (f has_derivative f' (t, x)) (at (t, x))"
  assumes f'': "\<And>t x dt dx. ((\<lambda>(t, x). f' (t, x) (dt, dx)) has_derivative f'' (t, x) (dt, dx)) (at (t, x))"
begin

(*   let ?d = "\<lambda>tq tr. (h^3) *\<^sub>R ((1/6)*\<^sub>R ?q tq - (p / 4) *\<^sub>R ?r tr)" *)
(*
  let ?r = "\<lambda>a. f'' (t + a, x t + a *\<^sub>R f (t, x t)) (1, f (t, x t))
      (1, f (t, x t))"
  let ?q = "\<lambda>s. f'' (s, x s) (1, f (s, x s)) (1, f (s, x s)) +
     f' (s, x s) (0, f' (s, x s) (1, f (s, x s)))"

*)
interpretation local_rect UNIV
  by (metis local_rect_openI open_UNIV)

lemma numeral_6_eq_6: "6 = Suc (Suc (Suc (Suc (Suc (Suc 0)))))"
  by (simp add: eval_nat_numeral)

lemma scale_back: "(r, r *\<^sub>R x) = r *\<^sub>R (1, x)"
  by simp

lemma rk2_consistent_traj_set:
  fixes t
  assumes ht: "0 \<le> h" "t + h \<le> u"
  assumes T: "{t..u} \<subseteq> T" and X0_nonempty: "X0 \<noteq> {}" and X_nonempty: "X \<noteq> {}" and convex_X: "convex X"
  assumes x': "\<And>s. s \<in> {t..u} \<Longrightarrow> (x has_vector_derivative f (s, x s)) (at s within {t..u})"
  assumes x: "\<And>s. s \<in> {t..u} \<Longrightarrow> x s \<in> X"
  assumes f_set_bounded: "bounded F" "\<And>t x x0. t\<in>T \<Longrightarrow> x0 \<in> X0 \<Longrightarrow> x\<in>X \<Longrightarrow> (x0, x, f (t, x)) \<in> F"
  assumes p: "0 < p" "p \<le> 1"
  assumes in_X0: "x t \<in> X0"
  assumes step_in: "x t + (h * p) *\<^sub>R f (t, x t) \<in> X"
  assumes R_ivl: "is_interval R"
  assumes heun_remainder_bounded:
      "\<And>tt xttt fxttt ttt t xt. tt \<in> {0..1} \<Longrightarrow> (xt, xttt, fxttt) \<in> F \<Longrightarrow> h ^ 3 *\<^sub>R
    ((1 / 6) *\<^sub>R (f'' (ttt, xttt) (1, fxttt) (1, fxttt) + f' (ttt, xttt) (0, f' (ttt, xttt) (1, fxttt))) -
     (p / 4) *\<^sub>R f'' (t + h * (p * tt), xt + (h * (p * tt)) *\<^sub>R f (t, xt)) (1, f (t, xt)) (1, f (t, xt))) \<in> R"
  shows "x (t + h) - discrete_evolution (rk2_increment p f) (t + h) t (x t) \<in> R"
proof cases
  assume "h = 0"
  from T ht have "t \<in> T" by auto
  hence F_nonempty: "F \<noteq> {}" using X_nonempty X0_nonempty
    unfolding ex_in_conv[symmetric]
    by (force intro!: f_set_bounded)
  with F_nonempty X_nonempty heun_remainder_bounded[of 0] `h = 0` have "0 \<in> R" by auto
  from this `h = 0` show ?thesis
    by (auto simp: discrete_evolution_def)
next
  assume "h \<noteq> 0"
  from this ht have "t < u" by simp
  have [simp]: "p \<noteq> 0" using p by simp

  let ?r = "\<lambda>a. f'' (t + a, x t + a *\<^sub>R f (t, x t)) (1, f (t, x t))
      (1, f (t, x t))"
  let ?q = "\<lambda>s. f'' (s, x s) (1, f (s, x s)) (1, f (s, x s)) +
     f' (s, x s) (0, f' (s, x s) (1, f (s, x s)))"

  let ?d = "\<lambda>tq tr. (h^3) *\<^sub>R ((1/6)*\<^sub>R ?q tq - (p / 4) *\<^sub>R ?r tr)"
  term ?d
  from ht have line_subset: "(\<lambda>ta. t + ta * h) ` {0..1} \<subseteq> {t..u}"
    by (auto intro!: order_trans[OF add_left_mono[OF mult_left_le_one_le]])
  hence line_in: "\<And>s. 0 \<le> s \<Longrightarrow> s \<le> 1 \<Longrightarrow> t + s * h \<in> {t..u}"
    by (rule set_mp) auto
  from ht have subset: "{t .. t + h} \<subseteq> {t .. u}" by simp
  let ?T = "{t..u}"
  from ht have subset: "{t .. t + h} \<subseteq> {t .. u}" by simp
  from `t < u` have "t \<in> ?T" by auto
  from `t < u` have tx: "t \<in> T" "x t \<in> X" using T ht x by auto
  from tx assms have "0 \<le> norm (f (t, x t))" by simp
  have x_diff: "\<And>s. s \<in> ?T \<Longrightarrow> x differentiable at s within ?T"
    by (rule differentiableI, rule x'[simplified has_vector_derivative_def])
  let ?p = "(\<lambda>t. f' (t, x t) (1, f (t, x t)))"
  {
    fix s assume "s \<in> ?T"
    then interpret bounded_linear "f' (s, x s)"
      apply (intro has_derivative_bounded_linear)
      apply (rule f')
      done
    have "?p s = f' (s, x s) (1, 0) + (\<Sum>i\<in>Basis. (f (s, x s) \<bullet> i) *\<^sub>R f' (s, x s) (0, i))"
      by (auto simp: plus_prod_def fst_setsum snd_setsum euclidean_representation
        add[symmetric] setsum[symmetric] scaleR[symmetric])
  } note * = this
  note f'[derivative_intros]
  note f''[derivative_intros]
  note x'[simplified has_vector_derivative_def, derivative_intros]
  interpret second_derivative UNIV f f' f''
    using f' f''
    by unfold_locales (auto intro!: f' `t < u`  simp: split_beta')
  have p': "\<And>s. s \<in> {t .. u} \<Longrightarrow> (?p
      has_vector_derivative
      (?q s)) (at s within ?T)"
    apply (rule has_vector_derivative_imp)
    apply assumption
    apply (rule *[symmetric])
    apply (auto intro!: derivative_eq_intros simp: has_vector_derivative_def)
    apply (rule diff_chain_within[simplified o_def]) back
    apply (auto intro!: derivative_eq_intros simp: has_vector_derivative_def)
    apply (rule has_derivative_subset)
    apply (rule derivative_intros) apply simp apply (simp add: x)
    apply (rule diff_chain_within[simplified o_def]) back
    apply (auto intro!: derivative_eq_intros simp: has_vector_derivative_def)
    apply (rule has_derivative_subset)
    apply (rule derivative_intros) apply simp
    apply (rule diff_chain_within[simplified o_def]) back
    apply (auto intro!: derivative_eq_intros simp: has_vector_derivative_def)
    apply (rule has_derivative_subset)
    apply (rule derivative_intros) apply simp apply (simp add: x)
  proof
    case (goal1 s sa)
    have b: "bilinear (f'' (s, x s))" using bilinear_second by (force simp: x goal1)
    interpret linear "f' (s, x s)" using linear_first by (force simp: x goal1)
    show "f'' (s, x s) (1, 0) (sa, sa *\<^sub>R f (s, x s)) +
       (\<Sum>i\<in>Basis. (f (s, x s) \<bullet> i) *\<^sub>R f'' (s, x s) (0, i) (sa, sa *\<^sub>R f (s, x s)) +
                   (f' (s, x s) (sa, sa *\<^sub>R f (s, x s)) \<bullet> i) *\<^sub>R f' (s, x s) (0, i)) =
      sa *\<^sub>R (f'' (s, x s) (1, f (s, x s)) (1, f (s, x s)) + f' (s, x s) (0, f' (s, x s) (1, f (s, x s))))"
      apply (simp add: algebra_simps setsum.distrib)
      apply (simp add: bilinear_lmul[OF b, symmetric] scaleR[symmetric] setsum[symmetric] setsum_prod euclidean_representation)
      apply (subst bilinear_lsetsum[OF b, symmetric])
      apply (simp add: euclidean_representation setsum_prod bilinear_ladd[OF b, symmetric])
    proof -
      case goal1
      have "f'' (s, x s) (1, f (s, x s)) (sa, sa *\<^sub>R f (s, x s)) =
         sa *\<^sub>R f'' (s, x s) (1, f (s, x s)) (1, f (s, x s))"
        by (subst bilinear_rmul[OF b, symmetric]) simp
      also have "\<dots> =  f'' (s, x s) (sa, sa *\<^sub>R f (s, x s)) (1, f (s, x s))"
        by (subst bilinear_lmul[OF b, symmetric]) simp
      finally show ?case .
    qed
  qed

  def diff \<equiv> "\<lambda>n::nat. if n = 0 then x else if n = 1 then \<lambda>t. f (t, x t) else if n = 2 then ?p
    else ?q"
  have diff_0[simp]: "diff 0 = x" by (simp add: diff_def)
  {
    fix m::nat and ta::real
    assume mta: "m < 3" "t \<le> ta" "ta \<le> t + h"
    have image_subset: "(\<lambda>xa. (xa, x xa)) ` {t..u} \<subseteq> {t..u} \<times> X"
      using assms by auto
    note has_derivative_in_compose[where f="(\<lambda>xa. (xa, x xa))" and g = f, derivative_intros]
    note has_derivative_subset[OF _ image_subset, derivative_intros]
    note f'[derivative_intros]
    note x'[simplified has_vector_derivative_def, derivative_intros]
    have [simp]: "\<And>c x'. c *\<^sub>R f' (ta, x ta) x' = f' (ta, x ta) (c *\<^sub>R x')"
      using mta ht assms by (auto intro!: f' linear_cmul[symmetric] has_derivative_linear)
    have "((\<lambda>t. f (t, x t)) has_vector_derivative f' (ta, x ta) (1, f (ta, x ta))) (at ta within {t..u})"
      unfolding has_vector_derivative_def
      using assms ht mta by (auto intro!: derivative_eq_intros has_derivative_within_subset[OF f'])
    hence "(diff m has_vector_derivative diff (Suc m) ta) (at ta within {t..t + h})"
      using mta ht
      by (auto simp: diff_def intro!: has_vector_derivative_within_subset[OF _ subset] x' p')
  } note diff = this
  obtain ttt where ttt: "ttt \<in> Basis \<rightarrow> {t <..< t + h}"
    and taylor: "x (t + h) = x t +
         h *\<^sub>R (f (t, x t)) +
         (h ^ 2 / (2::nat)) *\<^sub>R (?p t) +
         (h ^ 3 / (6::nat)) *\<^sub>R (\<Sum>x\<in>Basis. ((?q (ttt x)) \<bullet> x) *\<^sub>R x)"
    using taylor_up_within_vector[of 3 t "t + h" diff x t, OF _ _ diff]
    using ht `h\<noteq>0`
    by (auto simp: field_simps Pi_iff numeral_2_eq_2 numeral_3_eq_3
      numeral_6_eq_6 power2_eq_square diff_def scaleR_setsum_right)
  hence taylor: "x (t + h) = x t + h *\<^sub>R f (t, x t) + (h\<^sup>2/2) *\<^sub>R ?p t +
      (h^3/6) *\<^sub>R (\<Sum>i\<in>Basis. ((?q (ttt i)) \<bullet> i) *\<^sub>R i)"
    by simp

  interpret higher_derivative_lrect UNIV f 1 "\<lambda>xs. case xs of [] \<Rightarrow> f' | [d] \<Rightarrow> (\<lambda>x. f'' x d)"
    using f' f'' by unfold_locales (auto simp: split_beta')
  {
    fix x h::"real*'a"
    assume line_in: "(\<lambda>s. x + s *\<^sub>R h) ` {0..1} \<subseteq> UNIV"
    obtain i::'a where "i \<in> Basis" by (metis SOME_Basis)
    have *: "(\<lambda>(xs, x). x # xs) ` ({[]} \<times> Basis) = {[i]|i. i \<in> Basis}"
      by (auto intro!: image_eqI[where x="([], i)"])
    from multivariate_taylor[OF line_in, simplified] obtain tt
    where tt: "tt\<in>Basis \<rightarrow> {0<..<1}"
    and "f (x + h) = f x + (\<Sum>d\<in>Basis. (h \<bullet> d) *\<^sub>R f' x d) +
        (\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) *
          ((\<Sum>d\<in>Basis. \<Sum>ds\<in>{[i] |i. i \<in> Basis}.
            fold (\<lambda>d. op * (h \<bullet> d)) ds (h \<bullet> d) *\<^sub>R (case ds of [] \<Rightarrow> f' | [d] \<Rightarrow> \<lambda>y. f'' y d)
            (x + tt i *\<^sub>R h) d) \<bullet> i)) *\<^sub>R i)"
      (is "_ = _ + _ + (\<Sum>i\<in>Basis. (_ * ?s i) *\<^sub>R _)")
      by (auto simp add: *)
    note this(2)
    also have
      "(\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) * ((\<Sum>d\<in>Basis. \<Sum>ds\<in>{[i] |i. i \<in> Basis}.
            fold (\<lambda>d. op * (h \<bullet> d)) ds (h \<bullet> d) *\<^sub>R (case ds of [] \<Rightarrow> f' | [d] \<Rightarrow> \<lambda>y. f'' y d)
            (x + tt i *\<^sub>R h) d) \<bullet> i)) *\<^sub>R i) =
       (\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) * ((\<Sum>d\<in>Basis. \<Sum>d' \<in> Basis.
         fold (\<lambda>d. op * (h \<bullet> d)) [d'] (h \<bullet> d) *\<^sub>R (\<lambda>y. f'' y d') (x + tt i *\<^sub>R h) d) \<bullet> i)) *\<^sub>R i)"
      apply (rule setsum.cong, simp)
      apply (rule arg_cong) back back back
      apply (rule setsum.cong, simp)
      apply (rule setsum.reindex_cong[where l="\<lambda>x. [x]"])
      apply (auto )
      by (metis (erased, lifting) inj_on_def list.inject)
    also have "\<dots> = (\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) * ((\<Sum>d\<in>Basis. \<Sum>d'\<in>Basis. (h \<bullet> d' * (h \<bullet> d)) *\<^sub>R f'' (x + tt i *\<^sub>R h) d' d) \<bullet> i)) *\<^sub>R i)"
      by simp
    also
    interpret linear "f' x"
      using tt `i \<in> Basis`
      by (auto simp: Pi_iff zero_prod_def intro!:  linear_first)
    have "(\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) * ((\<Sum>d\<in>Basis. \<Sum>d'\<in>Basis. (h \<bullet> d' * (h \<bullet> d)) *\<^sub>R f'' (x + tt i *\<^sub>R h) d' d) \<bullet> i)) *\<^sub>R i)=
      (\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) * ((\<Sum>d\<in>Basis. \<Sum>d'\<in>Basis. f'' (x + tt i *\<^sub>R h) ((h \<bullet> d') *\<^sub>R d') ((h \<bullet> d) *\<^sub>R d)) \<bullet> i)) *\<^sub>R i)"
      apply (intro setsum.cong, rule refl)
      apply (rule arg_cong) back back back
      apply (intro setsum.cong, rule refl, rule refl)
      apply (auto simp: in_Basis_prod_iff)
      done
    also
    have "\<dots> = (\<Sum>i\<in>Basis. (inverse (real (Suc (Suc 0))) * (f'' (x + tt i *\<^sub>R h) h h \<bullet> i)) *\<^sub>R i)"
      apply (rule setsum.cong, rule refl)
    proof -
      case (goal1 i)
      have bi: "bilinear (f'' (x + tt i *\<^sub>R h))"
        using tt `i \<in> Basis`
        by (force simp: Pi_iff intro!: bilinear_second set_rev_mp[OF _ line_in])
      thus ?case
        apply (subst bilinear_lsetsum[OF bi, symmetric])
        apply (subst bilinear_rsetsum[OF bi, symmetric])
        apply simp
        apply (simp add: euclidean_representation)
        done
    qed
    also have "\<dots> = (inverse (real (Suc (Suc 0)))) *\<^sub>R (\<Sum>i\<in>Basis. ((f'' (x + tt i *\<^sub>R h) h h \<bullet> i)) *\<^sub>R i)"
      by (subst scaleR_setsum_right) simp
    also have "(\<Sum>d\<in>Basis. (h \<bullet> d) *\<^sub>R f' x d) = f' x h"
      apply (subst scaleR[symmetric])
      apply (subst setsum[symmetric])
      apply (subst euclidean_representation)
      ..
    finally have "\<exists>tt \<in> Basis \<rightarrow> {0<..<1}. f (x + h) = f x + f' x h +
      (\<Sum>i\<in>Basis. (f'' (x + tt i *\<^sub>R h) h h \<bullet> i) *\<^sub>R i) /\<^sub>R 2"
      using tt by (auto simp add: scaleR_setsum_right inverse_eq_divide real_of_nat_Suc)
  } note f_taylor = this

  let ?k = "\<lambda>t. f ((t, x t) + (h * p) *\<^sub>R (1, f (t, x t)))"
  have line_in: "(\<lambda>s. (t, x t) + s *\<^sub>R ((h * p) *\<^sub>R (1, f (t, x t)))) ` {0..1} \<subseteq> UNIV"
    by simp
(*    apply (auto intro!: x mult_nonneg_nonneg ht less_imp_le[OF p(1)])
    apply (subst (2) mult.commute)
    apply (subst mult.assoc[symmetric])
    apply (rule line_in[simplified, THEN conjunct2])
    apply (rule mult_nonneg_nonneg) apply simp apply (simp add: less_imp_le p)
    apply (rule mult_le_oneI)
    apply (simp add: p)
  proof -
    case goal1
    have "(1 - s) *\<^sub>R x t + s *\<^sub>R (x t + (h * p) *\<^sub>R f (t, x t)) \<in> X"
      apply (rule convexD)
      using `t < u` goal1
      by (auto intro!: convex_X x simp: algebra_simps step_in)
    thus ?case by (simp add: algebra_simps)
  qed *)
  from f_taylor[OF line_in, simplified]
  obtain tt where tt: "tt\<in>Basis \<rightarrow> {0<..<1}"
    and k: "?k t = f (t, x t) + f' (t, x t) (h * p, (h * p) *\<^sub>R f (t, x t)) +
     (1 / 2) *\<^sub>R
     (\<Sum>i\<in>Basis. (f'' (t + tt i * (h * p),
                       x t + (tt i * (h * p)) *\<^sub>R f (t, x t))
                   (h * p, (h * p) *\<^sub>R f (t, x t))
                   (h * p, (h * p) *\<^sub>R f (t, x t)) \<bullet>
                  i) *\<^sub>R
                 i)"
    by auto
(*
  from line_in have "(\<lambda>s. x t + s *\<^sub>R ((h * p) *\<^sub>R f (t, x t))) ` {0..1} \<subseteq> X"
    by auto
  hence line_in': "\<And>i. i \<in> Basis \<Longrightarrow> x t + (h * (p * tt i)) *\<^sub>R f (t, x t) \<in> X"
    apply (rule set_mp)
    using tt
    apply (force simp: Pi_iff)
    done *)
  have line_in'': "\<And>i. i \<in> Basis \<Longrightarrow> 0 \<le> h * (p * tt i)"
    using tt
    by (auto intro!: mult_nonneg_nonneg `0 \<le> h` less_imp_le[OF `0 < p`] simp: Pi_iff)
(*  from line_in have "(\<lambda>s. t + s *\<^sub>R (h * p)) ` {0..1} \<subseteq> {t .. u}"
    by auto
  hence "\<And>i. i \<in> Basis \<Longrightarrow> t + h * (p * tt i) \<in> {t .. u}"
    apply (rule set_mp)
    using tt
    apply (force simp: Pi_iff)
    done
  hence line_in''': "\<And>i. i \<in> Basis \<Longrightarrow> t + h * (p * tt i) \<le> u"
    by auto *)
  have rk2: "discrete_evolution (rk2_increment p f) (t + h) t (x t) =
    x t + h *\<^sub>R f (t, x t) +
    h *\<^sub>R (( - 1 / (2 * p)) *\<^sub>R f (t, x t)) +
    h *\<^sub>R ((1 / (2 * p)) *\<^sub>R ?k t)"
    (is "_ = _ ?rk2 t")
    apply (simp add: rk2_increment_def discrete_evolution_def )
    apply (subst rk_increment_def)
    apply (subst rk_eval.simps)
    apply (simp del: rk_eval.simps add:  numeral_2_eq_2)
    apply (subst rk_eval.simps)
    apply (simp del: rk_eval.simps add: algebra_simps)
    done
  also have "\<dots> = x t + h *\<^sub>R f (t, x t) + (h / (2 * p)) *\<^sub>R (f' (t, x t) ((h * p), (h * p) *\<^sub>R f (t, x t))) +
     (h^3 * p / 4) *\<^sub>R (\<Sum>i\<in>Basis. (f'' (t + tt i * (h * p),
                       x t + (tt i * (h * p)) *\<^sub>R f (t, x t))
                   (1, f (t, x t))
                   (1, f (t, x t)) \<bullet>
                  i) *\<^sub>R
                 i)"
    using k
    apply (auto simp: algebra_simps)
    apply (subst setsum.cong, rule refl)
    apply (subst scale_back)+
    apply (subst bilinear_lmul[of "f'' (t + h * (p * tt xa), x t + (h * (p * tt xa)) *\<^sub>R f (t, x t))" for xa])
    apply (rule bilinear_second)
    using `t < u`
    apply (auto intro!:  line_in'' )[1]
    apply (subst bilinear_rmul[of "f'' (t + h * (p * tt xa), x t + (h * (p * tt xa)) *\<^sub>R f (t, x t))" for xa])
    apply (rule bilinear_second)
    using `t < u`
    apply (auto intro!: line_in'')[1]
    apply (rule refl)
    apply (auto simp: scaleR_setsum_right numeral_3_eq_3 algebra_simps)[1]
    done
  also
  have "linear (f' (t, x t))"
    using linear_Ds[of "[]" "(t, x t)"] `t < u` x by (auto )
  have "(h / (2 * p)) *\<^sub>R f' (t, x t) (h * p, (h * p) *\<^sub>R f (t, x t)) = (h\<^sup>2 / 2) *\<^sub>R ?p t"
    apply (subst linear_cmul[of "f' (t, x t)", symmetric], fact)+
    apply (subst scaleR_prod_def)
    by (auto simp: power2_eq_square)
  finally
  have discrete_taylor: "discrete_evolution (rk2_increment p f) (t + h) t (x t) =
    x t + h *\<^sub>R f (t, x t) +
    (h\<^sup>2 / 2) *\<^sub>R f' (t, x t) (1, f (t, x t)) +
    (h ^ 3 * p / 4) *\<^sub>R (\<Sum>i\<in>Basis. ((?r (tt i * (h * p))) \<bullet> i) *\<^sub>R i)"
    .
  from taylor discrete_taylor
  have "x (t + h) - discrete_evolution (rk2_increment p f) (t + h) t (x t) =
    (h^3) *\<^sub>R ((1/6)*\<^sub>R (\<Sum>i\<in>Basis. ((?q (ttt i)) \<bullet> i) *\<^sub>R i) - (p / 4) *\<^sub>R (\<Sum>i\<in>Basis. ((?r (tt i * (h * p))) \<bullet> i) *\<^sub>R i))"
    by (auto simp: algebra_simps)
  also have "\<dots> \<in> R"
    apply (rule mem_box_componentwiseI[OF R_ivl])
    apply (simp add: algebra_simps )
    apply (subst scaleR_add_left[symmetric])
    apply simp
  proof -
    case goal1
    have " h ^ 3 *
      (f' (ttt i, x (ttt i)) (0, f' (ttt i, x (ttt i)) (1, f (ttt i, x (ttt i)))) \<bullet> i +
       f'' (ttt i, x (ttt i)) (1, f (ttt i, x (ttt i))) (1, f (ttt i, x (ttt i))) \<bullet> i) /
      6 -
      p * (h ^ 3 *
           (f'' (t + h * (p * tt i), x t + (h * (p * tt i)) *\<^sub>R f (t, x t))
             (1, f (t, x t)) (1, f (t, x t)) \<bullet>
            i)) /
      4 = (h ^ 3 *\<^sub>R
      ((1 / 6) *\<^sub>R
       (f'' (ttt i, x (ttt i)) (1, f (ttt i, x (ttt i))) (1, f (ttt i, x (ttt i))) +
        f' (ttt i, x (ttt i)) (0, f' (ttt i, x (ttt i)) (1, f (ttt i, x (ttt i))))) -
      (p / 4) *\<^sub>R f'' (t + h * (p * tt i), x t + (h * (p * tt i)) *\<^sub>R f (t, x t)) (1, f (t, x t)) (1, f (t, x t)))) \<bullet> i"
      by (auto simp: algebra_simps)
    also have "\<dots> \<in> (\<lambda>x. x \<bullet> i) ` R"
      apply (rule image_eqI[OF refl])
      apply (rule heun_remainder_bounded)
      using tt `i \<in> Basis` `t \<in> T` `t < u` ttt ht T in_X0
      apply (force simp add: Pi_iff intro!: f_set_bounded x)+
      done
    finally show ?case .
  qed
  finally show ?thesis .
qed

end

locale derivative_norm_bounded = derivative_on_convex T X f f' for T and X::"'a::euclidean_space set" and f f' +
  fixes B B'
  assumes X_bounded: "bounded X"
  assumes f_bounded: "\<And>t x. t\<in>T \<Longrightarrow> x\<in>X \<Longrightarrow> norm (f (t,x)) \<le> B"
  assumes f'_bounded: "\<And>t x. t\<in>T \<Longrightarrow> x\<in>X \<Longrightarrow> onorm (f' (t,x)) \<le> B'"
begin

lemma f_bound_nonneg: "0 \<le> B"
proof -
  from nonempty obtain t x where "t \<in> T" "x \<in> X" by auto
  have "0 \<le> norm (f (t, x))" by simp
  also have "\<dots> \<le> B" by (rule f_bounded) fact+
  finally show ?thesis .
qed

lemma f'_bound_nonneg: "0 \<le> B'"
proof -
  from nonempty f_bounded ex_norm_eq_1[where 'a="real*'a"]
    obtain t x and d::"real*'a" where tx: "t \<in> T" "x \<in> X" "norm d = 1" by auto
  have "0 \<le> norm (f' (t, x) d)" by simp
  also have "... \<le> B'"
    apply (rule order_trans)
     apply (rule onorm[OF has_derivative_bounded_linear[OF f']])
       using tx
       apply (auto intro!: f'_bounded f' has_derivative_linear)
    done
  finally show ?thesis .
qed

sublocale g: global_lipschitz _ _ _ B'
proof
  fix t assume "t \<in> T"
  show "lipschitz X (\<lambda>x. f (t, x)) B'"
  proof (rule lipschitzI)
    show "0 \<le> B'" using f'_bound_nonneg .
    fix x y
    let ?I = "T \<times> X"
    have "convex ?I" by (intro convex convex_Times)
    moreover have "\<forall>x\<in>?I. (f has_derivative f' x) (at x within ?I)" "\<forall>x\<in>?I. onorm (f' x) \<le> B'"
      using f' f'_bounded
      by (auto simp add: intro!: f'_bounded has_derivative_linear)
    moreover assume "x \<in> X" "y \<in> X"
    with `t \<in> T` have "(t, x) \<in> ?I" "(t, y) \<in> ?I" by simp_all
    ultimately have "norm (f (t, x) - f (t, y)) \<le> B' * norm ((t, x) - (t, y))"
      by (rule differentiable_bound)
    thus "dist (f (t, x)) (f (t, y)) \<le> B' * dist x y"
      by (simp add: dist_norm norm_Pair)
  qed
qed

definition euler_C::"real" where "euler_C = (sqrt DIM('a) * (B' * (B + 1) / 2))"

lemma euler_C_nonneg: "euler_C \<ge> 0"
 using f_bounded f_bound_nonneg f'_bound_nonneg
 by (simp add: euler_C_def)

sublocale derivative_set_bounded T X f f' "X \<times> cball 0 B"
    "cbox (- (B' * (B + 1)) *\<^sub>R One) ((B' * (B + 1)) *\<^sub>R One)"
proof
  show "bounded (X \<times> cball 0 B)" using X_bounded by (auto intro!: bounded_Times)
  show "is_interval (cbox (-(B' * (B + 1)) *\<^sub>R One) ((B' * (B + 1)) *\<^sub>R One::'a))"
    "bounded (cbox (-(B' * (B + 1)) *\<^sub>R One) ((B' * (B + 1)) *\<^sub>R One::'a))"
    by (auto intro!: is_interval_cbox bounded_cbox)
  fix t x assume "t \<in> T" "x \<in> X"
  thus "(x, f (t, x)) \<in> X \<times> cball 0 B"
    by (auto simp: dist_norm f_bounded)
next
  fix t and x d::'a assume "t \<in> T" "(x, d) \<in> X \<times> cball 0 B"
  hence "x \<in> X" "norm d \<le> B" by (auto simp: dist_norm)
  have "norm (f' (t, x) (1, d)) \<le> onorm (f' (t, x)) * norm (1::real, d)"
    by (auto intro!: onorm has_derivative_bounded_linear f' `t \<in> T` `x \<in> X`)
  also have "\<dots> \<le> B' * (B + 1)"
    by (auto intro!: mult_mono f'_bounded f_bounded `t \<in> T` `x \<in> X` f'_bound_nonneg
      order_trans[OF norm_Pair_le] `norm d \<le> B`)
  finally have "f' (t, x) (1, d) \<in> cball 0 (B' * (B + 1))"
    by (auto simp: dist_norm)
  also note cball_in_cbox
  finally show "f' (t, x) (1, d) \<in> cbox (- (B' * (B + 1)) *\<^sub>R One) ((B' * (B + 1)) *\<^sub>R One)"
    by simp
qed

lemma euler_consistent_traj:
  fixes t
  assumes T: "{t..u} \<subseteq> T"
  assumes x': "\<And>s. s \<in> {t..u} \<Longrightarrow> (x has_vector_derivative f (s, x s)) (at s within {t..u})"
  assumes x: "\<And>s. s \<in> {t..u} \<Longrightarrow> x s \<in> X"
  shows "consistent x t u euler_C 1 (euler_increment f)"
proof
  fix h::real
  assume ht: "0 < h" "t + h \<le> u" hence "t < u" "0 < h\<^sup>2 / 2" by simp_all
  from euler_consistent_traj_set ht T x' x
  have "x (t + h) - discrete_evolution (euler_increment f) (t + h) t (x t) \<in>
      op *\<^sub>R (h\<^sup>2 / 2) ` cbox (- (B' * (B + 1)) *\<^sub>R One) ((B' * (B + 1)) *\<^sub>R One)"
    by auto
  also have "\<dots> = cbox (- ((h\<^sup>2 / 2) * (B' * (B + 1))) *\<^sub>R One) (((h\<^sup>2 / 2) * (B' * (B + 1))) *\<^sub>R One)"
    using f_bound_nonneg f'_bound_nonneg
    apply (auto simp add: image_smult_cbox box_eq_empty)
    by (metis eucl_less_le_not_le gr_implies_not0 less_eq_real_def
      linordered_field_class.sign_simps(42) numeral_One one_power2 pos_add_strict power_0_left
      real_0_le_add_iff realpow_square_minus_le zero_less_numeral)
  also
  note centered_cbox_in_cball
  finally show "dist (x (t + h)) (discrete_evolution (euler_increment f) (t + h) t (x t))
      \<le> euler_C * h ^ (1 + 1)"
    by (auto simp: euler_C_def dist_norm algebra_simps norm_minus_commute power2_eq_square)
qed

end

locale grid_from = grid +
  fixes t0
  assumes grid_min: "t0 = t 0"

locale euler_consistent =
  has_solution i +
  derivative_norm_bounded T X' f B f' B'
  for i::"'a::euclidean_space ivp" and t X' B f' B' +
  fixes r e
  assumes domain_subset: "X \<subseteq> X'"
  assumes interval: "T = {t0 - e .. t0 + e}"
  assumes lipschitz_area: "\<And>t. t \<in> T \<Longrightarrow> cball (solution t) \<bar>r\<bar> \<subseteq> X'"
begin

lemma euler_consistent_solution:
  fixes t'
  assumes t': "t' \<in> {t0 .. t0 + e}"
  shows "consistent solution t' (t0 + e) euler_C 1 (euler_increment f)"
proof (rule euler_consistent_traj)
  show "{t'..t0 + e} \<subseteq> T" using t' interval by simp
  fix s
  assume "s \<in> {t'..t0 + e}" hence "s \<in> T" using `{t'..t0 + e} \<subseteq> T` by auto
  show "(solution has_vector_derivative f (s, solution s)) (at s within {t'..t0 + e})"
    by (rule has_vector_derivative_within_subset[OF _ `{t'..t0 + e} \<subseteq> T`]) (rule solution(2)[OF `s \<in> T`])
  have "solution s \<in> ivp_X i" by (rule solution(3)[OF `s \<in> T`])
  thus "solution s \<in> X'" using domain_subset ..
qed

end

sublocale euler_consistent \<subseteq>
  consistent_one_step t0 "t0 + e" solution "euler_increment f" 1 euler_C r B'
proof
  show "0 < (1::nat)" by simp
  show "0 \<le> euler_C" using euler_C_nonneg by simp
  show "0 \<le> B'" using lipschitz_nonneg[OF lipschitz] iv_defined by simp
  fix s x assume s: "s \<in> {t0 .. t0 + e}"
  show "consistent solution s (t0 + e) euler_C 1 (euler_increment f)"
    using interval s f_bounded f'_bounded f'
      strip
    by (intro euler_consistent_solution) auto
  fix h
  assume "h\<in>{0..t0 + e - s}"
  have "lipschitz X' (euler_increment f h s) B'"
    using s lipschitz interval strip
    by (auto intro!: euler_lipschitz)
  thus "lipschitz (cball (solution s) \<bar>r\<bar>) (euler_increment f h s) B'"
    using lipschitz_area s interval
    apply (auto intro: )
    apply (rule lipschitz_subset)
    apply assumption
    apply simp
    done
qed

subsection {* Euler method is convergent *}
text{*\label{sec:rk-euler-conv}*}

locale max_step1 = grid +
  fixes t1 L B r
  assumes max_step: "\<And>j. t j \<le> t1 \<Longrightarrow> max_stepsize j \<le> \<bar>r\<bar> * L / B / (exp (L * (t1 - t 0) + 1) - 1)"

sublocale max_step1 < max_step: max_step t t1 1 L B r
using max_step by unfold_locales simp_all

locale euler_convergent =
  euler_consistent + max_step1 t "t0 + e" B' euler_C r +
  assumes grid_from: "t0 = t 0"

sublocale euler_convergent \<subseteq>
  convergent_one_step t0 "t0 + e" solution "euler_increment f" 1 euler_C r B' t
  by unfold_locales (simp add: grid_from)

subsection {* Euler method on Rectangle is convergent *}
text{*\label{sec:rk-euler-conv-on-rect}*}

locale ivp_rectangle_bounded_derivative = solution_in_cylinder "i::'a::euclidean_space ivp" e b B +
  derivative_norm_bounded T "cbox (x0 - (b + \<bar>r\<bar>) *\<^sub>R One) (x0 + (b + \<bar>r\<bar>) *\<^sub>R One)" f f' B B' for i e b r B f' B'

sublocale ivp_rectangle_bounded_derivative \<subseteq> unique_on_cylinder i e b B B'
  "cbox (x0 - (b + \<bar>r\<bar>) *\<^sub>R One) (x0 + (b + \<bar>r\<bar>) *\<^sub>R One)"
  using b_pos cball_in_cbox[of x0 "b + abs r"]
  by unfold_locales (auto simp: cylinder intro!: scaleR_mono One_nonneg)

sublocale ivp_rectangle_bounded_derivative \<subseteq>
  euler_consistent i t "cbox (x0 - (b + \<bar>r\<bar>) *\<^sub>R One) (x0 + (b + \<bar>r\<bar>) *\<^sub>R One)" f' B B' r e
proof
  show "X \<subseteq> cbox (x0 - (b + \<bar>r\<bar>) *\<^sub>R One) (x0 + (b + \<bar>r\<bar>) *\<^sub>R One)" using lipschitz_on_domain .
  fix t assume "t \<in> T"
  have "cball (solution t) \<bar>r\<bar> \<subseteq> cball x0 (b + \<bar>r\<bar>)"
    using solution_in_D[of t] cylinder `t \<in> T`
    by (auto intro: cball_trans simp: interval)
  also note cball_in_cbox
  finally show "cball (solution t) \<bar>r\<bar> \<subseteq> cbox (x0 - (b + \<bar>r\<bar>) *\<^sub>R One) (x0 + (b + \<bar>r\<bar>) *\<^sub>R One)" .
qed (simp_all add: interval)

locale euler_on_rectangle =
  ivp_rectangle_bounded_derivative i e b r B f' B' +
  grid_from t t0 +
  max_step1 t "t0 + e" B' euler_C r
  for i::"'a::euclidean_space ivp" and t e b r B f' B'

sublocale euler_on_rectangle \<subseteq>
  convergent: euler_convergent i t "cbox (x0 - (b + \<bar>r\<bar>) *\<^sub>R One) (x0 + (b + \<bar>r\<bar>) *\<^sub>R One)" f' B B' r e
proof unfold_locales
qed (rule grid_min)

lemma "B \<ge> (0::real) \<Longrightarrow> 0 \<le> (exp (B + 1) - 1)" by (simp add: algebra_simps)

context euler_on_rectangle begin

lemma convergence:
  assumes "t j \<le> t0 + e"
  shows "dist (solution (t j)) (euler f x0 t j)
    \<le> sqrt DIM('a) * (B + 1) / 2 * (exp (B' * e + 1) - 1) * max_stepsize j"
proof -
  have "dist (solution (t j)) (euler f x0 t j)
    \<le> sqrt DIM('a) * (B + 1) / 2 * B' / B' * ((exp (B' * e + 1) - 1) * max_stepsize j)"
    using assms convergence[OF assms] f'_bound_nonneg
    unfolding euler_C_def
    by (simp add: euler_def grid_min[symmetric] solution_t0 ac_simps)
  also have "\<dots> \<le> sqrt DIM('a) * (B + 1) / 2 * ((exp (B' * e + 1) - 1) * max_stepsize j)"
    using f_bound_nonneg f'_bound_nonneg
    by (auto intro!: mult_right_mono mult_nonneg_nonneg max_stepsize_nonneg add_nonneg_nonneg
      simp: le_diff_eq)
  finally show ?thesis by simp
qed

end

subsection {* Stability and Convergence of Approximate Euler *}
text{*\label{sec:rk-euler-stable}*}

locale euler_rounded_on_rectangle =
  ivp_rectangle_bounded_derivative i e1' b r B f' B' +
  grid: grid_from t t0' +
  max_step_r_2: max_step1 t "t0 + e2'" B' euler_C "r/2"
  for i::"'a::executable_euclidean_space ivp" and t :: "nat \<Rightarrow> real" and t0' e1' e2'::real and x0' :: 'a
  and b r B f' B' +
  fixes g::"(real\<times>'a)\<Rightarrow>'a" and e::int
  assumes t0_float: "t0 = t0'"
  assumes ordered_bounds: "e1' \<le> e2'"
  assumes approx_f_e: "\<And>j x. t j \<le> t0 + e1' \<Longrightarrow> dist (f (t j, x)) ((g (t j, x))) \<le> sqrt (DIM('a)) * 2 powr -e"
  assumes initial_error: "dist x0 (x0') \<le> euler_C / B' * (exp 1 - 1) * stepsize 0"
  assumes rounding_error: "\<And>j. t j \<le> t0 + e1' \<Longrightarrow> sqrt (DIM('a)) * 2 powr -e \<le> euler_C / 2 * stepsize j"
begin

lemma approx_f: "t j \<le> t0 +  e1' \<Longrightarrow> dist (f (t j, x)) ((g (t j, x)))
    \<le> euler_C / 2 * stepsize j"
  using approx_f_e[of j x] rounding_error[of j] by auto

lemma t0_le: "t 0 \<le> t0 + e1'"
  unfolding grid_min[symmetric] t0_float[symmetric]
  by (metis atLeastAtMost_iff interval iv_defined(1))

end

sublocale euler_rounded_on_rectangle \<subseteq> grid': grid_from t t0'
  using grid t0_float grid_min by unfold_locales auto

sublocale euler_rounded_on_rectangle \<subseteq> max_step_r: max_step1 t "t0 + e2'" B' euler_C r
proof unfold_locales
  fix j
  assume "(t j) \<le> t0 + e2'"
  moreover with grid_mono[of 0 j] have "t 0 \<le> t0 + e2'" by (simp add: less_eq_float_def)
  ultimately show "max_stepsize j
        \<le> \<bar>r\<bar> * B' / euler_C / (exp (B' * (t0 + e2' - (t 0)) + 1) - 1)"
    using max_step_mono_r lipschitz B_nonneg f'_bound_nonneg
    by (auto simp: less_eq_float_def euler_C_def mult_nonneg_nonneg)
qed

lemma max_step1_mono:
  assumes "t 0 \<le> t1"
  assumes "t1 \<le> t2"
  assumes "0 \<le> a"
  assumes "0 \<le> b"
  assumes ms2: "max_step1 t t2 a b c"
  shows "max_step1 t t1 a b c"
proof -
  interpret t2: max_step1 t t2 a b c using ms2 .
  show ?thesis
  proof
    fix j
    assume "t j \<le> t1" hence "t j \<le> t2" using assms by simp
    hence "t2.max_stepsize j \<le> \<bar>c\<bar> * a / b / (exp (a * (t2 - t 0) + 1) - 1)" (is "_ \<le> ?x t2")
      by (rule t2.max_step)
    also have "\<dots> \<le> ?x t1"
      using assms
      by (cases "b = 0") (auto intro!: divide_left_mono mult_mono abs_ge_zero add_increasing
        mult_pos_pos add_strict_increasing2 simp: le_diff_eq less_diff_eq)
    finally show "t2.max_stepsize j \<le> ?x t1" .
  qed
qed

sublocale euler_rounded_on_rectangle \<subseteq> max_step_r1: max_step1 t "t0 + e1'" B' euler_C r
  by (rule max_step1_mono[of t, OF t0_le add_left_mono[OF ordered_bounds] f'_bound_nonneg euler_C_nonneg])
    unfold_locales

sublocale euler_rounded_on_rectangle \<subseteq> c: euler_on_rectangle i t e1' b r B f' B'
  using t0_float grid_min by unfold_locales simp

sublocale euler_rounded_on_rectangle \<subseteq>
  consistent_one_step "t 0" "t0 + e1'" solution "euler_increment f" 1 euler_C r B'
  using consistent_nonneg consistent lipschitz_nonneg lipschitz_incr t0_float grid_min
  by unfold_locales simp_all

sublocale euler_rounded_on_rectangle \<subseteq> max_step1 t "t0 + e1'" B' euler_C "r / 2"
  by (rule max_step1_mono[of t, OF t0_le add_left_mono[OF ordered_bounds] f'_bound_nonneg euler_C_nonneg])
    unfold_locales

sublocale euler_rounded_on_rectangle \<subseteq>
  one_step:
  rounded_one_step t "t0 + e1'" solution "euler_increment f" 1 euler_C r B' "euler_increment' e g" x0'
proof
  fix h j x assume "t j \<le> t0 + e1'"
  have "dist (euler_increment f (h) (t j) (x))
        ((euler_increment' e g h (t j) x)) =
    dist (f (t j, x)) ((eucl_down e (g (t j, x))))"
    by (simp add: euler_increment euler_float_increment)
  also
  have "... \<le>
    dist (f (t j, x)) ((g (t j, x))) +
    dist ((g (t j, x))) ((eucl_down e (g (t j, x))))"
    by (rule dist_triangle)
  also
  from approx_f[OF `t j \<le> t0 + e1'`]
  have "dist (f (t j, x)) ((g (t j, x))) \<le>
    euler_C / 2 * stepsize j" .
  also
  from eucl_truncate_down_correct[of "g (t j, x)" e]
  have "dist ((g (t j, x))) ((eucl_down e (g (t j, x)))) \<le> sqrt (DIM('a)) * 2 powr - e" by simp
  also
  have "sqrt (DIM('a)) * 2 powr -e \<le> euler_C / 2 * stepsize j"
    using rounding_error `t j \<le> t0 + e1'` .
  finally
  have "dist (euler_increment f (h) (t j) (x)) ((euler_increment' e g h (t j) x)) \<le> euler_C * stepsize j"
    by arith
  thus "dist (euler_increment f h (t j) (x)) ((euler_increment' e g h (t j) x)) \<le> euler_C * stepsize j ^ 1"
    by simp
qed (insert initial_error grid_min solution_t0, simp_all)

context euler_rounded_on_rectangle begin

lemma stability:
  assumes "t j \<le> t0 + e1'"
  shows "dist (euler' e g x0' t j) (euler f x0 t j) \<le>
    sqrt DIM('a) * (B + 1) / 2 * (exp (B' * e1' + 1) - 1) * max_stepsize j"
proof -
  have "dist ((euler' e g x0' t j)) (euler f x0 t j) \<le>
    sqrt DIM('a) * (B + 1) / 2 * B' / B' * (exp (B' * e1' + 1) - 1) * max_stepsize j"
    using assms stability[OF assms]
    unfolding grid_min[symmetric] solution_t0 euler_C_def
    by (auto simp add: euler_def euler'_def t0_float)
  also have "\<dots> \<le> sqrt DIM('a) * (B + 1) / 2 * ((exp (B' * e1' + 1) - 1) * max_stepsize j)"
    using f_bound_nonneg f'_bound_nonneg
    by (auto intro!: mult_right_mono mult_nonneg_nonneg max_stepsize_nonneg add_nonneg_nonneg
      simp: le_diff_eq)
  finally show ?thesis by simp
qed

lemma convergence_float:
  assumes "t j \<le> t0 + e1'"
  shows "dist (solution (t j)) (euler' e g x0' t j) \<le>
    sqrt DIM('a) * (B + 1) * (exp (B' * e1' + 1) - 1) * max_stepsize j"
proof -
  have "dist (solution ((t j))) ((euler' e g x0' t j)) \<le>
    dist (solution ((t j)))
    (euler f x0 t j) +
    dist ((euler' e g x0' t j)) (euler f x0 t j)"
    by (rule dist_triangle2)
  also have "dist (solution ((t j)))
     (euler f x0 t j) \<le>
    sqrt DIM('a) * (B + 1) / 2 * (exp (B' * e1' + 1) - 1) * max_stepsize j"
    using assms convergence[OF assms] t0_float by simp
  also have "dist ((euler' e g x0' t j)) (euler f x0 t j) \<le>
    sqrt DIM('a) * (B + 1) / 2 * (exp (B' * e1' + 1) - 1) * max_stepsize j"
    using assms stability by simp
  finally
  have "dist (solution ((t j))) ((euler' e g x0' t j))
     \<le> sqrt DIM('a) * (B + 1) / 2 * (exp (B' * (e1') + 1) - 1) *
       max_stepsize j +
       sqrt DIM('a) * (B + 1) / 2 * (exp (B' * (e1') + 1) - 1) *
       max_stepsize j" by simp
  thus ?thesis by (simp add: field_simps)
qed

end

end
