# Theory Runge_Kutta

section‹Runge-Kutta methods›
theory Runge_Kutta
imports
"HOL-Analysis.Analysis"
One_Step_Method
"HOL-Library.Float"
Affine_Arithmetic.Executable_Euclidean_Space
Ordinary_Differential_Equations.Multivariate_Taylor
begin

subsection ‹aux›

lemma scale_back: "(r, r *⇩R x) = r *⇩R (1, x)" "(0, r *⇩R x) = r *⇩R (0, x)"
by simp_all

lemma integral_normalize_bounds:
fixes s t::real
assumes "t ≤ s"
assumes "f integrable_on {t .. s}"
shows [symmetric]: "(s - t) *⇩R integral {0 .. 1} (λx. f ((s - t) *⇩R x + t)) = integral {t..s} f"
proof cases
assume "s > t"
hence "s - t ≠ 0" "0 ≤ s - t" by simp_all
from assms have "(f has_integral integral {t .. s} f) (cbox t s)"
by (auto simp: integrable_integral)
from has_integral_affinity[OF this ‹s - t ≠ 0›, of t]
have "((λx. f ((s - t) * x + t)) has_integral (1 / ¦s - t¦) *⇩R integral {t..s} f)
((λx. (x - t) / (s - t))  {t..s})"
using ‹s > t›
also
have "t < s ⟹ 0 ≤ x ⟹ x ≤ 1 ⟹ x * (s - t) + t ≤ s" for x
by (auto simp add: algebra_simps dest: mult_left_le_one_le[OF ‹0 ≤ s - t›])
then have "((λx. (x - t) / (s - t))  {t..s}) = {0 .. 1}"
using ‹s > t›
by (auto intro!: image_eqI[where x="x * (s - t) + t" for x]
simp: divide_simps)
finally
have "integral {0..1} (λx. f ((s - t) * x + t)) = (1 / ¦s - t¦) *⇩R integral {t..s} f"
by (rule integral_unique)
then show ?thesis
using ‹s > t› by simp
qed (insert assms, simp)

lemma
has_integral_integral_eqI:
"f integrable_on s ⟹ integral s f = k ⟹ (f has_integral k) s"

lemma convex_scaleR_sum2:
assumes "x ∈ G" "y ∈ G" "convex G"
assumes "a ≥ 0" "b ≥ 0" "a + b ≠ 0"
shows "(a *⇩R x + b *⇩R y) /⇩R (a + b) ∈ G"
proof -
have "(a / (a + b)) *⇩R x + (b / (a + b)) *⇩R y ∈ G"
using assms
by (intro convexD) (auto simp: divide_simps)
then show ?thesis
by (auto simp: algebra_simps divide_simps)
qed

lemma sum_by_parts_ivt:
assumes "finite X"
assumes "convex G"
assumes "⋀i. i ∈ X ⟹ g i ∈ G"
assumes "⋀i. i ∈ X ⟹ 0 ≤ c i"
obtains y where "y ∈ G" "(∑x∈X. c x *⇩R g x) = sum c X *⇩R y" | "G = {}"
proof (atomize_elim, cases "sum c X = 0", goal_cases)
case pos: 2
let ?y = "(∑x∈X. (c x / sum c X) *⇩R g x)"
have "?y ∈ G" using pos
by (intro convex_sum)
(auto simp: sum_divide_distrib[symmetric]
intro!: divide_nonneg_nonneg assms sum_nonneg)
thus ?case
by (auto intro!: exI[where x = ?y] simp: scaleR_right.sum pos)
qed (insert assms, auto simp: sum_nonneg_eq_0_iff)

lemma
integral_by_parts_near_bounded_convex_set:
assumes f: "(f has_integral I) (cbox a b)"
assumes s: "((λx. f x *⇩R g x) has_integral P) (cbox a b)"
assumes G: "⋀x. x ∈ cbox a b ⟹ g x ∈ G"
assumes nonneg: "⋀x. x ∈ cbox a b ⟹ f x ≥ 0"
assumes convex: "convex G"
assumes bounded: "bounded G"
shows "infdist P ((*⇩R) I  G) = 0"
proof (rule dense_eq0_I, cases)
fix e'::real assume e0: "0 < e'"
assume "G ≠ {}"
from bounded obtain bnd where bnd: "⋀y. y ∈ G ⟹ norm y < bnd" "bnd > 0"
by (meson bounded_pos gt_ex le_less_trans norm_ge_zero)
define e where "e = min (e' / 2) (e' / 2 / bnd)"
have e: "e > 0" using e0
by (auto simp add: e_def intro!: divide_pos_pos ‹bnd > 0›)
from
has_integral[of f I a b,                THEN iffD1, OF f, rule_format, OF e]
has_integral[of "λx. f x *⇩R g x" P a b, THEN iffD1, OF s, rule_format, OF e]
obtain d1 d3
where d1: "gauge d1"
"⋀p. p tagged_division_of cbox a b ⟹ d1 fine p ⟹
norm ((∑(x, k)∈p. content k *⇩R f x) - I) < e"
and d3: "gauge d3"
"⋀p. p tagged_division_of cbox a b ⟹ d3 fine p ⟹
norm ((∑(x, k)∈p. content k *⇩R f x *⇩R g x) - P) < e"
by auto
define d where "d x = d1 x ∩ d3 x" for x
from d1(1) d3(1)
have "gauge d" by (auto simp add: d_def [abs_def])
from fine_division_exists[OF this, of a b]
obtain p where p: "p tagged_division_of cbox a b" "d fine p"
by metis
from tagged_division_of_finite[OF p(1)]
have "finite p" .

from ‹d fine p› have "d1 fine p" "d3 fine p"
by (auto simp: d_def [abs_def] fine_Int)
have f_less: "norm ((∑(x, k)∈p. content k *⇩R f x) - I) < e"
(is "norm (?f - I) < _")
by (rule d1(2)[OF p(1)]) fact
have "norm ((∑(x, k)∈p. content k *⇩R f x *⇩R g x) - P) < e"
(is "norm (?s - P) < _")
by (rule d3(2)[OF p(1)]) fact

hence "dist (∑(x, k)∈p. content k *⇩R f x *⇩R g x) P < e"
also
let ?h = "(λx k y. (content k * f x) *⇩R y)"
let ?s' = "λy. sum (λ(x, k). ?h x k y) p"
let ?g = "λ(x, k). g x"
let ?c = "λ(x, k). content k * f x"
have Pi: "⋀x. x ∈ p ⟹ ?g x ∈ G" "⋀x. x ∈ p ⟹ ?c x ≥ 0"
using nonneg G p
using tag_in_interval[OF p(1)]
by (auto simp: intro!: mult_nonneg_nonneg)
obtain y where y: "y ∈ G" "?s = ?s' y"
by (rule sum_by_parts_ivt[OF ‹finite p› ‹convex G› Pi])
(auto simp: split_beta' scaleR_sum_left ‹G ≠ {}›)
note this(2)
also have "(∑(x, k)∈p. (content k * f x) *⇩R y) = ?f *⇩R y"
by (auto simp: scaleR_left.sum intro!: sum.cong)
finally have "dist P ((∑(x, k)∈p. content k *⇩R f x) *⇩R y) ≤ e"
moreover have "dist (I *⇩R y) ((∑(x, k)∈p. content k *⇩R f x) *⇩R y) ≤ norm y * e"
using f_less
by (auto simp add: dist_real_def mult.commute [of _ "norm y"]
intro!: mult_left_mono)
ultimately
have "dist P (I *⇩R y) ≤ e + norm y * e"
with _ have "infdist P ((*⇩R) I  G) ≤ e + norm y * e"
using y(1)
by (intro infdist_le2) auto
also have "norm y * e < bnd * e"
by (rule mult_strict_right_mono)
(auto simp: ‹e > 0› less_imp_le intro!: bnd ‹y ∈ G›)
also have "bnd * e ≤ e' / 2"
using ‹e' > 0› ‹bnd > 0›
by (auto simp: e_def min_def divide_simps)
also have "e ≤ e' / 2" by (simp add: e_def)
also have "e' / 2 + e' / 2 = e'" by simp
finally show "¦infdist P ((*⇩R) I  G)¦ ≤ e'"
by (auto simp: infdist_nonneg)

lemma
integral_by_parts_in_bounded_closed_convex_set:
assumes f: "(f has_integral I) (cbox a b)"
assumes s: "((λx. f x *⇩R g x) has_integral P) (cbox a b)"
assumes G: "⋀x. x ∈ cbox a b ⟹ g x ∈ G"
assumes nonneg: "⋀x. x ∈ cbox a b ⟹ f x ≥ 0"
assumes bounded: "bounded G"
assumes closed: "closed G"
assumes convex: "convex G"
assumes nonempty: "cbox a b ≠ {}"
shows "P ∈ (*⇩R) I  G"
proof -
let ?IG = "(*⇩R) I  G"
from bounded closed have "bounded ?IG" "closed ?IG"
have "G ≠ {}" using nonempty G by auto
then show ?thesis
using ‹closed ?IG›
by (subst in_closed_iff_infdist_zero)
(auto intro!: assms compact_imp_bounded integral_by_parts_near_bounded_convex_set)
qed

lemma
integral_by_parts_in_bounded_set:
assumes f: "(f has_integral I) (cbox a b)"
assumes s: "((λx. f x *⇩R g x) has_integral P) (cbox a b)"
assumes nonneg: "⋀x. x ∈ cbox a b ⟹ f x ≥ 0"
assumes bounded: "bounded (g  cbox a b)"
assumes nonempty: "cbox a b ≠ {}"
shows "P ∈ (*⇩R) I  closure (convex hull (g  cbox a b))"
proof -
have "x ∈ cbox a b ⟹ g x ∈ closure (convex hull g  cbox a b)" for x
by (meson closure_subset hull_subset imageI subsetCE)
then show ?thesis
by (intro integral_by_parts_in_bounded_closed_convex_set[OF f s _ nonneg _ _ _ nonempty])
(auto intro!: bounded_closure bounded_convex_hull bounded convex_closure
simp: convex_convex_hull)
qed

lemma snd_imageI: "(a, b) ∈ R ⟹ b ∈ snd  R"
by force

lemma in_minus_Collect: "a ∈ A ⟹ b ∈ B ⟹ a - b ∈ {x - y|x y. x ∈ A ∧ y ∈ B}"
by blast

lemma closure_minus_Collect:
fixes A B::"'a::real_normed_vector set"
shows
"{x - y|x y. x ∈ closure A ∧ y ∈ closure B} ⊆ closure {x - y|x y. x ∈ A ∧ y ∈ B}"
proof -
have image: "(λ(x, y). x - y)  (A × B) = {x - y|x y. x ∈ A ∧ y ∈ B}" for A B::"'a set"
by auto
have "{x - y|x y. x ∈ closure A ∧ y ∈ closure B} = (λ(x, y). x - y)  closure (A × B)"
unfolding closure_Times
by (rule image[symmetric])
also have "… ⊆ closure ((λ(x, y). x - y)  (A × B))"
by (rule image_closure_subset)
(auto simp: split_beta' intro!: subsetD[OF closure_subset]
continuous_at_imp_continuous_on)
also note image
finally show ?thesis .
qed

lemma convex_hull_minus_Collect:
fixes A B::"'a::real_normed_vector set"
shows
"{x - y|x y. x ∈ convex hull A ∧ y ∈ convex hull B} = convex hull {x - y|x y. x ∈ A ∧ y ∈ B}"
proof -
have image: "(λ(x, y). x - y)  (A × B) = {x - y|x y. x ∈ A ∧ y ∈ B}" for A B::"'a set"
by auto
have "{x - y|x y. x ∈ convex hull A ∧ y ∈ convex hull B} = (λ(x, y). x - y)  (convex hull (A × B))"
unfolding convex_hull_Times
by (rule image[symmetric])
also have "… = convex hull ((λ(x, y). x - y)  (A × B))"
apply (rule convex_hull_linear_image)
by unfold_locales (auto simp: algebra_simps)
also note image
finally show ?thesis .
qed

lemma set_minus_subset:
"A ⊆ C ⟹ B ⊆ D ⟹ {a - b |a b. a ∈ A ∧ b ∈ B} ⊆ {a - b |a b. a ∈ C ∧ b ∈ D}"
by auto

lemma (in bounded_bilinear) bounded_image:
assumes "bounded (f  s)"
assumes "bounded (g  s)"
shows "bounded ((λx. prod (f x) (g x))  s)"
proof -
from nonneg_bounded obtain K
where K: "⋀a b. norm (prod a b) ≤ norm a * norm b * K" and "0 ≤ K"
by auto
from assms obtain F G
where F: "⋀x. x ∈ s ⟹ norm (f x) ≤ F"
and G: "⋀x. x ∈ s ⟹ norm (g x) ≤ G"
and nonneg: "0 ≤ F" "0 ≤ G"
by (auto simp: bounded_pos intro: less_imp_le)
have "norm (prod (f x) (g x)) ≤ F * G * K" if x: "x ∈ s" for x
using F[OF x] G[OF x] nonneg ‹0 ≤ K›
by (auto intro!: mult_mono mult_nonneg_nonneg order_trans[OF K])
thus ?thesis
by (auto simp: bounded_iff)
qed

lemmas bounded_scaleR_image = bounded_bilinear.bounded_image[OF bounded_bilinear_scaleR]
and bounded_blinfun_apply_image = bounded_bilinear.bounded_image[OF bounded_bilinear_blinfun_apply]

lemma bounded_plus_image:
fixes f::"'a ⇒ 'b::real_normed_vector"
assumes "bounded (f  s)"
assumes "bounded (g  s)"
shows "bounded ((λx. f x + g x)  s)"
proof -
from assms obtain F G
where F: "⋀x. x ∈ s ⟹ norm (f x) ≤ F"
and G: "⋀x. x ∈ s ⟹ norm (g x) ≤ G"
by (auto simp: bounded_iff)
have "norm (f x + g x) ≤ F + G" if x: "x ∈ s" for x
using F[OF x] G[OF x]
by norm
thus ?thesis
by (auto simp: bounded_iff)
qed

lemma bounded_uminus_image[simp]:
fixes f::"'a ⇒ 'b::real_normed_vector"
shows "bounded ((λx. - f x)  s) ⟷ bounded (f  s)"
apply (subst (2) bounded_uminus[symmetric])
unfolding image_image
by simp

lemma bounded_minus_image:
fixes f::"'a ⇒ 'b::real_normed_vector"
assumes "bounded (f  s)"
assumes "bounded (g  s)"
shows "bounded ((λx. f x - g x)  s)"
using bounded_plus_image[of f s "λx. - g x"] assms
by simp

lemma bounded_Pair_image:
fixes f::"'a ⇒ 'b::real_normed_vector"
fixes g::"'a ⇒ 'c::real_normed_vector"
assumes "bounded (f  s)"
assumes "bounded (g  s)"
shows "bounded ((λx. (f x, g x))  s)"
proof -
from assms obtain F G
where F: "⋀x. x ∈ s ⟹ norm (f x) ≤ F"
and G: "⋀x. x ∈ s ⟹ norm (g x) ≤ G"
by (auto simp: bounded_iff)
have "norm (f x, g x) ≤ F + G" if x: "x ∈ s" for x
using F[OF x] G[OF x]
by (intro order_trans[OF norm_Pair_le]) norm
thus ?thesis
by (auto simp: bounded_iff)
qed

lemma closed_scaleR_image_iff:
fixes f::"'a ⇒ 'b::real_normed_vector"
shows "closed ((λx. r *⇩R (f x))  R) ⟷ (r = 0 ∨ closed (f  R))" (is "?l ⟷ _ ∨ ?r")
proof safe
assume ?r
from closed_scaling[OF this] show ?l
by (auto simp: image_image)
next
assume l: ?l
{
assume "r ≠ 0"
with closed_scaling[OF l, of "inverse r"]
have "?r"
by (auto simp: image_image)
} then show "¬ ?r ⟹ r = 0" by blast

lemma closed_translation_iff[simp]:
fixes y::"'a::real_normed_vector"
shows "closed ((λx. f x + y)  S) ⟷ closed (f  S)"
"closed ((λx. y + f x)  S) ⟷ closed (f  S)"
using closed_translation[of "f  S" y] closed_translation[of "(+) y  f  S" "- y"]
by (auto simp: add.commute image_image cong: image_cong)

lemma closed_minus_translation_iff[simp]:
fixes y::"'a::real_normed_vector"
shows "closed ((λx. f x - y)  S) ⟷ closed (f  S)"
using closed_translation_iff(1)[of f "- y" S]
by simp

lemma convex_scaleR_image_iff:
fixes f::"'a ⇒ 'b::real_normed_vector"
shows "convex ((λx. r *⇩R (f x))  R) ⟷ (r = 0 ∨ convex (f  R))" (is "?l ⟷ _ ∨ ?r")
proof safe
assume ?r
from convex_scaling[OF this] show ?l
by (auto simp: image_image)
next
assume l: ?l
{
assume "r ≠ 0"
with convex_scaling[OF l, of "inverse r"]
have "?r"
by (auto simp: image_image)
} then show "¬ ?r ⟹ r = 0" by blast

lemma convex_translation_iff[simp]:
fixes y::"'a::real_normed_vector"
shows "convex ((λx. f x + y)  S) ⟷ convex (f  S)"
"convex ((λx. y + f x)  S) ⟷ convex (f  S)"
using convex_translation[of "f  S" y] convex_translation[of "(+) y  f  S" "- y"]
by (auto simp: add.commute image_image cong: image_cong)

lemma convex_minus_translation_iff[simp]:
fixes y::"'a::real_normed_vector"
shows "convex ((λx. f x - y)  S) ⟷ convex (f  S)"
using convex_translation_iff(1)[of f "- y" S]
by simp

text‹\label{sec:rk}›

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

declare sum.cong[fundef_cong]
fun rk_eval :: "(nat⇒nat⇒real) ⇒ (nat⇒real) ⇒ (real ⇒ 'a::real_vector ⇒ 'a) ⇒ real ⇒ real ⇒ 'a ⇒ nat ⇒ 'a" where
"rk_eval A c f t h x j =
f (t + h * c j) (x + h *⇩R (∑l=1 ..< j. A j l *⇩R rk_eval A c f t h x l))"

primrec rk_eval_dynamic :: "(nat⇒nat⇒real) ⇒ (nat⇒real) ⇒ (real×'a::{comm_monoid_add, scaleR} ⇒ 'a) ⇒ real ⇒ real ⇒ 'a ⇒ nat ⇒ (nat ⇒ 'a)" where
"rk_eval_dynamic A c f t h x 0 = (λ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 *⇩R (∑l=1..j. A (Suc j) l *⇩R K l))))"

definition rk_increment where
"rk_increment f s A b c h t x = (∑j=1..s. b j *⇩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 (∑j=1..s. b j *⇩R rk_eval A c f t h x j)"

definition euler_increment where
"euler_increment f = rk_increment f 1 (λi j. 0) (λi. 1) (λ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 (λi j. 0) (λi. 1) (λ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 (λi j. if i = 2 ∧ j = 1 then x else 0)
(λi. if i = 1 then 1 - 1 / (2 * x) else 1 / (2 * x)) (λ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 \label{sec:rk-euler-cons}›

lemma euler_increment:
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:
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:
assumes t: "t ∈ {t0..T}"
assumes lipschitz: "∀t∈{t0..T}. L-lipschitz_on D' (λx. f t x)"
shows "L-lipschitz_on D' (euler_increment f h t)"
using t lipschitz
by (simp add: lipschitz_on_def euler_increment del: One_nat_def)

lemma rk2_increment:
shows "rk2_increment p f h t x =
(1 - 1 / (p * 2)) *⇩R f t x +
(1 / (p * 2)) *⇩R f (t + h * p) (x + (h * p) *⇩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 \label{sec:setconsistent}›

context derivative_on_prod
begin

lemma euler_consistent_traj_set:
fixes t
assumes ht: "0 ≤ h" "t + h ≤ u"
assumes T: "{t..u} ⊆ T"
assumes x': "⋀s. s ∈ {t..u} ⟹ (x has_vector_derivative f s (x s)) (at s within {t..u})"
assumes x: "⋀s. s ∈ {t..u} ⟹ x s ∈ X"
assumes R: "⋀s. s ∈ {t..u} ⟹
discrete_evolution (euler_increment f) (t + h) t (x t) + (h⇧2 / 2) *⇩R (f' (s, x s)) (1, f s (x s)) ∈ R"
assumes bcc: "bounded R" "closed R" "convex R"
shows "x (t + h) ∈ R"
proof cases
assume "h = 0"
with assms show ?thesis
by (auto simp: discrete_evolution_def)
next
assume "h ≠ 0"
from this ht have "t < u" by simp
from ht have line_subset: "(λta. t + ta * h)  {0..1} ⊆ {t..u}"
by (auto intro!: order_trans[OF add_left_mono[OF mult_left_le_one_le]])
hence line_in: "⋀s. 0 ≤ s ⟹ s ≤ 1 ⟹ t + s * h ∈ {t..u}"
by (rule subsetD) auto
from ht have subset: "{t .. t + h} ⊆ {t .. u}" by simp
let ?T = "{t..u}"
from ht have subset: "{t .. t + h} ⊆ {t .. u}" by simp
from ‹t < u› have "t ∈ ?T" by auto
from ‹t < u› have tx: "t ∈ T" "x t ∈ X" using assms by auto
from tx assms have "0 ≤ norm (f t (x t))" by simp
have x_diff: "⋀s. s ∈ ?T ⟹ x differentiable at s within ?T"
by (rule differentiableI, rule x'[simplified has_vector_derivative_def])

note [continuous_intros] =
continuous_on_compose2[OF has_derivative_continuous_on[OF f'] continuous_on_Pair, simplified]
continuous_on_compose2[OF has_derivative_continuous_on[OF x'[unfolded has_vector_derivative_def]], simplified]

let ?p = "(λt. f' (t, x t) (1, f t (x t)))"
define diff where "diff ≡ λn::nat. if n = 0 then x else if n = 1 then λt. f t (x t) else ?p"
have diff_0[simp]: "diff 0 = x" by (simp add: diff_def)
have diff: "(diff m has_vector_derivative diff (Suc m) ta) (at ta within {t..t + h})"
if mta: "m < 2" "t ≤ ta" "ta ≤ t + h" for m::nat and ta::real
proof -
have image_subset: "(λxa. (xa, x xa))  {t..u} ⊆ {t..u} × X"
using assms by auto
note has_derivative_subset[OF _ image_subset, derivative_intros]
note f'[derivative_intros]
note x'[simplified has_vector_derivative_def, derivative_intros]
have [simp]: "⋀c x'. c *⇩R f' (ta, x ta) x' = f' (ta, x ta) (c *⇩R x')"
using mta ht assms
by (force intro!: f' linear_cmul[symmetric] has_derivative_linear)
have f_comp': "((λ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)
then show ?thesis
using mta ht
by (auto simp: diff_def intro!: has_vector_derivative_within_subset[OF _ subset] x')
qed

from Taylor_has_integral[of 2 diff x t "t + h", OF _ _ diff] ‹0 ≤ h›
have Taylor: "((λxa. (t + h - xa) *⇩R f' (xa, x xa) (1, f xa (x xa))) has_integral x (t + h) - (x t + h *⇩R f t (x t))) {t..t + h}"

have *: "h⇧2 / 2 = content {t..t + h} *⇩R (t + h) - (if t ≤ t + h then (t + h)⇧2 / 2 - t⇧2 / 2 else 0)"
using ‹0 ≤ h›
by (simp add: algebra_simps power2_eq_square divide_simps)
have integral: "((-) (t + h) has_integral h⇧2 / 2) (cbox t (t + h))"
unfolding * cbox_interval
using ‹0 ≤ h›
by (auto intro!: has_integral_diff ident_has_integral[THEN has_integral_eq_rhs]
has_integral_const_real[THEN has_integral_eq_rhs])
from Taylor_has_integral[of 2 diff x t "t + h", OF _ _ diff] ‹0 ≤ h›
have Taylor: "((λxa. (t + h - xa) *⇩R f' (xa, x xa) (1, f xa (x xa))) has_integral x (t + h) - (x t + h *⇩R f t (x t))) {t..t + h}"

define F' where "F' ≡ (λy. (2 / h⇧2) *⇩R (y - discrete_evolution (euler_increment f) (t + h) t (x t)))  R"
have "x (t + h) - (x t + h *⇩R f t (x t)) ∈ (*⇩R) (h⇧2 / 2)  F'"
apply (rule integral_by_parts_in_bounded_closed_convex_set[OF integral Taylor[unfolded interval_cbox]])
subgoal using R ‹h ≠ 0› ‹0 ≤ h› subset by (force simp: F'_def)
by (auto intro!: bounded_scaleR_image bounded_minus_image closed_injective_image_subspace bcc ‹0 ≤ h›
simp: F'_def image_constant_conv closed_scaleR_image_iff convex_scaleR_image_iff ‹h ≠ 0›)
then show ?thesis
using ‹h ≠ 0›
by (auto simp add: discrete_evolution_def euler_increment F'_def)
qed

end

lemma numeral_6_eq_6: "6 = Suc (Suc (Suc (Suc (Suc (Suc 0)))))"
by linarith

context begin

interpretation blinfun_syntax .

definition "heun_remainder1 x f f' f'' t h s
= (h ^ 3 / 6) *⇩R (f'' (h * s + t, x (h * s + t)) $(1::real, f (h * s + t, x (h * s + t)))$ (1::real, f (h * s + t, x (h * s + t))) +
f' (h * s + t, x (h * s + t)) $(0::real, f' (h * s + t, x (h * s + t))$ (1, f (h * s + t, x (h * s + t)))))"

definition "heun_remainder2 p x f f'' t h s =
(h ^ 3 * p / 4) *⇩R f'' (t + s * h * p, x t + (s * h * p) *⇩R f (t, (x t))) $(1::real, f (t, (x t)))$ (1::real, f (t, (x t)))"

lemma rk2_consistent_traj_set:
fixes x ::"real ⇒ 'a::banach" and t
and f' :: "real × 'a ⇒ (real × 'a) ⇒⇩L 'a"
and g' :: "real × 'a ⇒ (real × 'a) ⇒ 'a"
and f'' :: "real × 'a ⇒ (real × 'a) ⇒⇩L (real × 'a) ⇒⇩L 'a"
and g'' :: "real × 'a ⇒ (real × 'a) ⇒ (real × 'a) ⇒⇩L 'a"
assumes ht: "0 ≤ h" "t + h ≤ u"
assumes T: "{t..u} ⊆ T" and X_nonempty: "X ≠ {}" and convex_X: "convex X"
assumes x': "⋀s. s ∈ {t..u} ⟹ (x has_vector_derivative f (s, x s)) (at s within {t..u})"
assumes f': "⋀tx. tx ∈ T × X ⟹ (f has_derivative g' tx) (at tx)"
assumes f'': "⋀tx. tx ∈ T × X ⟹  (f' has_derivative g'' tx) (at tx)"
assumes g': "⋀tx. tx ∈ T × X ⟹ g' tx = f' tx"
assumes g'': "⋀tx. tx ∈ T × X ⟹  g'' tx = f'' tx"
assumes f''_bounded: "bounded (f''  (T × X))"
assumes x: "⋀s. s ∈ {t..u} ⟹ x s ∈ X"
assumes p: "0 < p" "p ≤ 1"
assumes step_in: "x t + (h * p) *⇩R f (t, (x t)) ∈ X"
assumes ccR: "convex R" "closed R"
assumes R: "⋀s1 s2. 0 ≤ s1 ⟹ s1 ≤ 1 ⟹ 0 ≤ s2 ⟹ s2 ≤ 1 ⟹
discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) +
heun_remainder1   x f f' f'' t h s1 -
heun_remainder2 p x f    f'' t h s2 ∈ R"
shows "x (t + h) ∈ R"
proof cases
assume "h = 0"
with R[of 0 0]
show ?thesis
by (auto simp: discrete_evolution_def heun_remainder1_def heun_remainder2_def)
next
have f': "⋀tx. tx ∈ T × X ⟹ (f has_derivative blinfun_apply (f' tx)) (at tx)"
using f' g'
by simp
have f'': "⋀tx. tx ∈ T × X ⟹  (f' has_derivative blinfun_apply (f'' tx)) (at tx)"
using f'' g''
by simp
assume "h ≠ 0"
from this ht have "t < u" by simp
have [simp]: "p ≠ 0" using p by simp
from ‹h ≥ 0› ‹h ≠ 0› have "h > 0" by simp

let ?r = "λa. f'' (t + a, x t + a *⇩R f (t, x t)) (1, f (t, x t))
(1, f (t, x t))"
let ?q = "λ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 = "λtq tr. (h^3) *⇩R ((1/6)*⇩R ?q tq - (p / 4) *⇩R ?r tr)"

from ht have line_subset: "(λta. t + ta * h)  {0..1} ⊆ {t..u}"
by (auto intro!: order_trans[OF add_left_mono[OF mult_left_le_one_le]])
hence line_in: "⋀s. 0 ≤ s ⟹ s ≤ 1 ⟹ t + s * h ∈ {t..u}"
by (rule subsetD) auto
from ht have subset: "{t .. t + h} ⊆ {t .. u}" by simp
let ?T = "{t..u}"
from ht have subset: "{t .. t + h} ⊆ {t .. u}" by simp
from ‹t < u› have "t ∈ ?T" by auto
from ‹t < u› have tx: "t ∈ T" "x t ∈ X" using T ht x by auto

from tx assms have "0 ≤ norm (f (t, x t))" by simp
have x_diff: "⋀s. s ∈ ?T ⟹ x differentiable at s within ?T"
by (rule differentiableI, rule x'[simplified has_vector_derivative_def])
let ?p = "(λt. f' (t, x t) (1, f (t, x t)))"
note f'[derivative_intros]
note f''[derivative_intros]
note x'[simplified has_vector_derivative_def, derivative_intros]

have x_cont: "continuous_on {t..u} x"
by (rule continuous_on_vector_derivative) (rule x')
have f_cont: "continuous_on (T × X) f"
apply (rule has_derivative_continuous_on)
apply (rule has_derivative_at_withinI)
by (rule assms)
have f'_cont: "continuous_on (T × X) f'"
apply (rule has_derivative_continuous_on)
apply (rule has_derivative_at_withinI)
by (rule assms)
note [continuous_intros] =
continuous_on_compose2[OF x_cont]
continuous_on_compose2[OF f_cont]
continuous_on_compose2[OF f'_cont]

from f' f''
have f'_within: "tx ∈ T × X ⟹ (f has_derivative f' tx) (at tx within T × X)"
and f''_within: "tx ∈ T × X ⟹ (f' has_derivative f'' tx) (at tx within T × X)" for tx
by (auto intro: has_derivative_at_withinI)

from f'' have f''_within: "tx ∈ T × X ⟹ (f' has_derivative ($) (f'' tx)) (at tx within T × X)" for tx by (auto intro: has_derivative_at_withinI) note [derivative_intros] = has_derivative_in_compose2[OF f'_within] has_derivative_in_compose2[OF f''_within] have p': "⋀s. s ∈ {t .. u} ⟹ (?p has_vector_derivative ?q s) (at s within ?T)" unfolding has_vector_derivative_def using T x by (auto intro!: derivative_eq_intros simp: scale_back blinfun.bilinear_simps algebra_simps simp del: scaleR_Pair) define diff where "diff n = (if n = 0 then x else if n = 1 then λt. f (t, x t) else if n = 2 then ?p else ?q)" for n :: nat have diff_0[simp]: "diff 0 = x" by (simp add: diff_def) have diff: "(diff m has_vector_derivative diff (Suc m) ta) (at ta within {t..t + h})" if mta: "m < 3" "t ≤ ta" "ta ≤ t + h" for m::nat and ta::real proof - have image_subset: "(λxa. (xa, x xa))  {t..u} ⊆ {t..u} × X" using assms by auto note has_derivative_in_compose[where f="(λ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]: "⋀c x'. c *⇩R f' (ta, x ta) x' = f' (ta, x ta) (c *⇩R x')" using mta ht assms T x by (force intro!: f' linear_cmul[symmetric] has_derivative_linear) have "((λ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 T x by (force intro!: derivative_eq_intros has_derivative_subset[OF f']) then show ?thesis using mta ht by (auto simp: diff_def intro!: has_vector_derivative_within_subset[OF _ subset] x' p') qed from Taylor_has_integral[of 3 diff x t "t + h", OF _ _ diff] have "((λx. ((t + h - x) ^ 2 / 2) *⇩R diff 3 x) has_integral x (t + h) - x t - h *⇩R (f (t, x t)) - (h ^ 2 / (2::nat)) *⇩R (?p t)) (cbox t (t + h))" using ht ‹h≠0› by (auto simp: field_simps of_nat_Suc Pi_iff numeral_2_eq_2 numeral_3_eq_3 numeral_6_eq_6 power2_eq_square diff_def scaleR_sum_right) from has_integral_affinity[OF this ‹h ≠ 0›, of t, simplified] have "((λx. ((h - h * x)⇧2 / 2) *⇩R diff 3 (h * x + t)) has_integral (1 / ¦h¦) *⇩R (x (t + h) - x t - h *⇩R f (t, x t) - (h⇧2 / 2) *⇩R f' (t, x t)$ (1, f (t, x t))))
((λx. x / h - t / h)  {t..t + h})"
by simp
also have "((λx. x / h - t / h)  {t..t + h}) = {0 .. 1}"
using ‹h ≠ 0› ‹h ≥ 0›
by (auto simp: divide_simps intro!: image_eqI[where x="x * h + t" for x])
finally have "((λx. ((h - h * x)⇧2 / 2) *⇩R diff 3 (h * x + t)) has_integral
(1 / ¦h¦) *⇩R (x (t + h) - x t - h *⇩R f (t, x t) - (h⇧2 / 2) *⇩R f' (t, x t) $(1, f (t, x t)))) {0..1}" . from has_integral_cmul[OF this, of h] have Taylor: "((λx. (1 - x)⇧2 *⇩R ((h^3 / 2) *⇩R ?q (h * x + t))) has_integral (x (t + h) - x t - h *⇩R f (t, x t) - (h⇧2 / 2) *⇩R f' (t, x t)$ (1, f (t, x t))))
{0..1}" (is "(?i_Taylor has_integral _) _")
using ‹h ≥ 0› ‹h ≠ 0›
by (simp add: diff_def divide_simps algebra_simps power2_eq_square power3_eq_cube)
have line_in': "h * y + t ∈ T"
"x (h * y + t) ∈ X"
"t ≤ h * y + t" "h * y + t ≤ u"
if "y ∈ cbox 0 1" for y
using line_in[of y] that T
by (auto simp: algebra_simps x)
let ?integral = "λx. x^3/3 - x^2 + x"
have intsquare: "((λx. (1 - x)⇧2) has_integral ?integral 1 - ?integral 0) (cbox 0 (1::real))"
unfolding cbox_interval
by (rule fundamental_theorem_of_calculus)
(auto intro!: derivative_eq_intros
simp: has_vector_derivative_def power2_eq_square algebra_simps)
have f_Taylor: "((λs. (1 - s) *⇩R f'' (x + s *⇩R h) h h) has_integral f (x + h) - f x - f' x $h) {0..1}" if line_in: "(λs. x + s *⇩R h)  {0..1} ⊆ T × X" for x h::"real*'a" proof - from that have *: "y ∈ T × X" if "y ∈ closed_segment x (x + h)" for y using that by (force simp: closed_segment_def algebra_simps intro: image_eqI[where x = "1 - x" for x]) define Df where "Df x i h1 h2 = (if i = 0 then f x else if i = 1 then f' x h2 else f'' x h2 h1)" for x h1 h2 and i::nat have "((λy. ((1 - y) ^ (2 - 1) / fact (2 - 1)) *⇩R Df (x + y *⇩R h) 2 h h) has_integral f (x + h) - (∑i<2. (1 / fact i) *⇩R Df x i h h)) {0..1}" apply (rule multivariate_Taylor_has_integral[of 2 Df h f x "T × X"]) subgoal by simp subgoal by (simp add: Df_def) subgoal premises prems for a i d proof - consider "i = 0" | "i = 1" | "i = 2" using prems by arith then show ?thesis by cases (use prems in ‹auto simp: Df_def[abs_def] blinfun.bilinear_simps intro!: derivative_eq_intros intro: *›) qed subgoal using "*" by blast done then show ?thesis by (simp add: Df_def eval_nat_numeral algebra_simps) qed let ?k = "λt. f ((t, x t) + (h * p) *⇩R (1, f (t, x t)))" have line_in: "(λs. (t, x t) + s *⇩R ((h * p) *⇩R (1, f (t, x t))))  {0..1} ⊆ T × X" proof (clarsimp, safe) fix s::real assume s: "0 ≤ s" "s ≤ 1" have "t + s * (h * p) = t + s * p * h" by (simp add: ac_simps) also have "… ∈ {t .. u}" using ‹0 < p› ‹p ≤ 1› s by (intro line_in) (auto intro!: mult_nonneg_nonneg mult_left_le_one_le mult_le_one) also note ‹… ⊆ T› finally show "t + s * (h * p) ∈ T" . show "x t + (s * (h * p)) *⇩R f (t, x t) ∈ X" using convexD_alt[OF ‹convex X› tx(2) step_in s] by (simp add: algebra_simps) qed from f_Taylor[OF line_in, simplified] have k: "((λs. (1 - s) *⇩R ((h⇧2 * p⇧2) *⇩R f'' (t + s * (h * p), x t + (s * (h * p)) *⇩R f (t, x t))$
(1, f (t, x t)) $(1, f (t, x t)))) has_integral ?k t - f (t, x t) - f' (t, x t)$ (h * p, (h * p) *⇩R f (t, x t))) {0..1}"
(is "(?i has_integral _) _")
unfolding scale_back blinfun.bilinear_simps
have rk2: "discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) =
x t + h *⇩R f (t, x t) -
(h / (2 * p)) *⇩R f (t, x t) +
(h / (p * 2)) *⇩R ?k t"
(is "_ = ?rk2 t")
unfolding rk2_increment_def discrete_evolution_def rk_increment_def
apply (subst rk_eval.simps)
supply rk_eval.simps[simp del]
apply (subst rk_eval.simps)
done
also have "… =
x t + h *⇩R f (t, x t) + (h / (2 * p)) *⇩R (f' (t, x t) ((h * p), (h * p) *⇩R f (t, x t)))
+ (h / (p * 2)) *⇩R integral {0 .. 1} ?i"
unfolding integral_unique[OF k]
also have "(h / (2 * p)) *⇩R f' (t, x t) (h * p, (h * p) *⇩R f (t, x t)) = (h⇧2 / 2) *⇩R ?p t"
by (simp add: scale_back blinfun.bilinear_simps power2_eq_square
del: scaleR_Pair)
finally
have "integral {0 .. 1} ?i =
(discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) -
x t - h *⇩R f (t, x t) -
(h⇧2 / 2) *⇩R ?p t) /⇩R (h / (p * 2))"
with _ have "(?i has_integral
(discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) -
x t - h *⇩R f (t, x t) -
(h⇧2 / 2) *⇩R ?p t) /⇩R
(h / (p * 2))) {0 .. 1}"
using k
by (intro has_integral_integral_eqI) (rule has_integral_integrable)
from has_integral_cmul[OF this, of "h / (p * 2)"]
have discrete_Taylor:
"((λs. (1 - s) *⇩R ((h^3 * p / 2) *⇩R
f'' (t + s * (h * p), x t + (s * (h * p)) *⇩R f (t, x t)) $(1, f (t, x t))$
(1, f (t, x t)))) has_integral
(discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) -
x t - h *⇩R f (t, x t) -
(h⇧2 / 2) *⇩R f' (t, x t) (1, f (t, x t)))) {0 .. 1}"
(is "(?i_dTaylor has_integral _) _")
using ‹h > 0›
by (simp add: algebra_simps diff_divide_distrib power2_eq_square power3_eq_cube)
have integral_minus: "((-) 1 has_integral 1/2) (cbox 0 (1::real))"
by (auto intro!: has_integral_eq_rhs[OF has_integral_diff] ident_has_integral)

have bounded_f: "bounded ((λxa. f (h * xa + t, x (h * xa + t)))  {0..1})"
using ‹0 ≤ h›
by (auto intro!: compact_imp_bounded compact_continuous_image continuous_intros mult_nonneg_nonneg
simp: line_in')
have bounded_f': "bounded ((λxa. f' (h * xa + t, x (h * xa + t)))  {0..1})"
using ‹0 ≤ h›
by (auto intro!: compact_imp_bounded compact_continuous_image continuous_intros
simp: line_in')
have bounded_f'': "bounded ((λxa. f'' (h * xa + t, x (h * xa + t)))  {0..1})"
apply (subst o_def[of f'', symmetric])
apply (subst image_comp[symmetric])
apply (rule bounded_subset[OF f''_bounded])
by (auto intro!: image_eqI line_in')
have bounded_f''_2:
"bounded ((λxa. f'' (t + xa * (h * p), x t + (xa * (h * p)) *⇩R f (t, x t)))  {0..1})"
apply (subst o_def[of f'', symmetric])
apply (subst image_comp[symmetric])
apply (rule bounded_subset[OF f''_bounded])
using line_in
by auto
have 1: "x (t + h) - x t - h *⇩R f (t, x t) - (h⇧2 / 2) *⇩R f' (t, x t) $(1, f (t, x t)) ∈ (*⇩R) (1 / 3)  closure (convex hull (λxa. (h ^ 3 / 2) *⇩R (f'' (h * xa + t, x (h * xa + t))$
(1,
f (h * xa + t, x (h * xa + t))) $(1, f (h * xa + t, x (h * xa + t))) + f' (h * xa + t, x (h * xa + t))$
(0,
f' (h * xa + t, x (h * xa + t)) $(1, f (h * xa + t, x (h * xa + t))))))  cbox 0 1)" by (rule rev_subsetD[OF integral_by_parts_in_bounded_set[OF intsquare Taylor[unfolded interval_cbox]]]) (auto intro!: bounded_scaleR_image bounded_plus_image bounded_blinfun_apply_image bounded_Pair_image bounded_f'' bounded_f' bounded_f simp: image_constant[of 0]) have 2: "discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) - x t - h *⇩R f (t, x t) - (h⇧2 / 2) *⇩R f' (t, x t)$ (1, f (t, x t)) ∈
(*⇩R) (1 / 2)  closure (convex hull
(λs. (h ^ 3 * p / 2) *⇩R
f''
(t + s * (h * p),
x t +
(s * (h * p)) *⇩R f (t, x t)) $(1, f (t, x t))$
(1, f (t, x t)))
cbox 0 1)"
by (rule integral_by_parts_in_bounded_set[OF integral_minus discrete_Taylor[unfolded interval_cbox]])
(auto intro!: bounded_scaleR_image bounded_blinfun_apply_image
bounded_f''_2 simp: image_constant[of 0])
have "x (t + h) - discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t) ∈
{a - b|a b.
a ∈
closure
(convex hull (*⇩R) (1/3) 
(λxa. (h ^ 3 / 2) *⇩R
(f'' (h * xa + t, x (h * xa + t)) $(1, f (h * xa + t, x (h * xa + t)))$
(1,
f (h * xa + t, x (h * xa + t))) +
f' (h * xa + t, x (h * xa + t)) $(0, f' (h * xa + t, x (h * xa + t))$
(1,
f (h * xa + t,
x (h * xa + t))))))
cbox 0 1) ∧
b ∈ closure (convex hull (*⇩R) (1 / 2) 
(λs. (h ^ 3 * p / 2) *⇩R
f''
(t + s * (h * p),
x t +
(s * (h * p)) *⇩R f (t, x t)) $(1, f (t, x t))$
(1, f (t, x t)))
cbox 0 1)}"
using in_minus_Collect[OF 1 2]
unfolding closure_scaleR convex_hull_scaling
by auto
also note closure_minus_Collect
also note convex_hull_minus_Collect
also
define F' where "F' ≡ (λy. y - discrete_evolution (rk2_increment p (λt x. f (t, x))) (t + h) t (x t))  R"
have "closure
(convex hull
{xa - y |xa y.
xa ∈ (*⇩R) (1 / 3)
(λxa. (h ^ 3 / 2) *⇩R
(f'' (h * xa + t, x (h * xa + t)) $(1, f (h * xa + t, x (h * xa + t)))$
(1, f (h * xa + t, x (h * xa + t))) +
f' (h * xa + t, x (h * xa + t)) $(0, f' (h * xa + t, x (h * xa + t))$
(1, f (h * xa + t, x (h * xa + t)))))) 
cbox 0 1 ∧
y ∈ (*⇩R) (1 / 2)
(λs. (h ^ 3 * p / 2) *⇩R
f'' (t + s * (h * p), x t + (s * (h * p)) *⇩R f (t, x t)) $(1, f (t, x t))$
(1, f (t, x t))) 
cbox 0 1}) ⊆ F'"
apply (rule closure_minimal)
apply (rule hull_minimal)
apply (safe)
subgoal for _ _ _ _ _ s1 s2
using R[of s1 s2]
by (force simp: F'_def heun_remainder1_def heun_remainder2_def algebra_simps)
by (auto intro!: ccR simp: F'_def)
finally
show "x (t + h) ∈ R"
by (auto simp: F'_def)
qed

end

locale derivative_norm_bounded = derivative_on_prod T X f f' for T and X::"'a::euclidean_space set" and f f' +
fixes B B'
assumes nonempty: "T ≠ {}" "X ≠ {}"
assumes X_bounded: "bounded X"
assumes convex: "convex T" "convex X"
assumes f_bounded: "⋀t x. t∈T ⟹ x∈X ⟹ norm (f t x) ≤ B"
assumes f'_bounded: "⋀t x. t∈T ⟹ x∈X ⟹ onorm (f' (t,x)) ≤ B'"
begin

lemma f_bound_nonneg: "0 ≤ B"
proof -
from nonempty obtain t x where "t ∈ T" "x ∈ X" by auto
have "0 ≤ norm (f t x)" by simp
also have "… ≤ B" by (rule f_bounded) fact+
finally show ?thesis .
qed

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

sublocale g?: global_lipschitz _ _ _ B'
proof
fix t assume "t ∈ T"
show "B'-lipschitz_on X (f t)"
proof (rule lipschitz_onI)
show "0 ≤ B'" using f'_bound_nonneg .
fix x y
let ?I = "T × X"
have "convex ?I" by (intro convex convex_Times)
moreover have "⋀x. x ∈ ?I ⟹ ((λ(t, x). f t x) has_derivative f' x) (at x within ?I)"
"⋀x. x ∈ ?I ⟹ onorm (f' x) ≤ B'"
using f' f'_bounded  by (auto simp add: intro!: f'_bounded has_derivative_linear)
moreover assume "x ∈ X" "y ∈ X"
with ‹t ∈ T› have "(t, x) ∈ ?I" "(t, y) ∈ ?I" by simp_all
ultimately have "norm ((λ(t, x). f t x) (t, x) - (λ(t, x). f t x) (t, y)) ≤ B' * norm ((t, x) - (t, y))"
by (rule differentiable_bound)
then show "dist (f t x) (f t y) ≤ B' * dist x y"
qed
qed

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

lemma euler_C_nonneg: "euler_C TYPE('a) ≥ 0"
using f_bounded f_bound_nonneg f'_bound_nonneg

lemma euler_consistent_traj:
fixes t
assumes T: "{t..u} ⊆ T"
assumes x': "(x has_vderiv_on (λs. f s (x s))) {t..u}"
assumes x: "⋀s. s ∈ {t..u} ⟹ x s ∈ X"
shows "consistent x t u (euler_C (TYPE('a))) 1 (euler_increment f)"
proof
from x' have x': "⋀s. s ∈ {t..u} ⟹ (x has_vector_derivative f s (x s)) (at s within {t..u})"
fix h::real
assume ht: "0 < h" "t + h ≤ u" hence "t < u" "0 < h⇧2 / 2" by simp_all
let ?d = "discrete_evolution (euler_increment f) (t + h) t (x t)"
have "x (t + h) ∈ (λb. ?d + (h⇧2 / 2) *⇩R b)  cbox (- (B' * (B + 1)) *⇩R One) ((B' * (B + 1)) *⇩R One)"
proof (rule euler_consistent_traj_set[OF _ ‹t + h ≤ u› T x' x])
fix s
assume "s ∈ {t .. u}"
then have "?d + (h⇧2 / 2) *⇩R (f' (s, x s)) (1, f s (x s)) ∈
(λb. ?d + (h⇧2 / 2) *⇩R b)  ((λs. (f' (s, x s)) (1, f s (x s))) {t .. u})"
by auto
also have "… ⊆ (λb. ?d + (h⇧2 / 2) *⇩R b)  cbox (- (B' * (B + 1)) *⇩R One) ((B' * (B + 1)) *⇩R One)"
proof (rule image_mono, safe)
fix s assume "s ∈ {t .. u}"
with T have "norm (f' (s, x s) (1, f s (x s))) ≤ onorm (f' (s, x s)) * norm (1::real, f s (x s))"
by (force intro!: onorm has_derivative_bounded_linear f' x)
also have "… ≤ B' * (B + 1)"
using T ‹s ∈ _›
by (force intro!: mult_mono f'_bounded f_bounded f'_bound_nonneg x order_trans[OF norm_Pair_le])
finally have "f' (s, x s) (1, f s (x s)) ∈ cball 0 (B' * (B + 1))"
by (auto simp: dist_norm mem_cball)
also note cball_in_cbox
finally show "f' (s, x s) (1, f s (x s)) ∈ cbox (- (B' * (B + 1)) *⇩R One) ((B' * (B + 1)) *⇩R One)"
by simp
qed
finally
show "?d + (h⇧2 / 2) *⇩R (f' (s, x s)) (1, f s (x s))
∈ (λb. ?d + (h⇧2 / 2) *⇩R b)  cbox (- (B' * (B + 1)) *⇩R One) ((B' * (B + 1)) *⇩R One)" .
qed (auto intro!: less_imp_le[OF ‹0 < h›] bounded_plus_image bounded_scaleR_image bounded_cbox
closed_scaling convex_scaling convex_box simp: image_constant_conv closed_translation_iff)
then have "x (t + h) - discrete_evolution (euler_increment f) (t + h) t (x t) ∈
(*⇩R) (h⇧2 / 2)  cbox (- (B' * (B + 1)) *⇩R One) ((B' * (B + 1)) *⇩R One)"
by auto
also have "… = cbox (- ((h⇧2 / 2) * (B' * (B + 1))) *⇩R One) (((h⇧2 / 2) * (B' * (B + 1))) *⇩R One)"
using f_bound_nonneg f'_bound_nonneg
by (auto simp add: image_smult_cbox box_eq_empty mult_less_0_iff)
also
note centered_cbox_in_cball
finally show "dist (x (t + h)) (discrete_evolution (euler_increment f) (t + h) t (x t))
≤ euler_C(TYPE('a)) * h ^ (1 + 1)"
by (auto simp: euler_C_def dist_norm algebra_simps norm_minus_commute power2_eq_square mem_cball)
qed

lemma derivative_norm_bounded_subset:
assumes X'_ne: "X' ≠ {}" and X'_subset: "X' ⊆ X" and  "convex X'"
shows "derivative_norm_bounded T X' f f' B B'"
proof -
interpret derivative_on_prod T X' f f'
using X'_subset
by (rule derivative_on_prod_subset)
show ?thesis
using X'_subset
by unfold_locales
(auto intro!: nonempty X'_ne bounded_subset[OF X_bounded X'_subset] convex f_bounded f'_bounded
‹convex X'›)
qed

end

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

locale euler_consistent =
derivative_norm_bounded T X' f f' B B'
for T f X X' B f' B' +
fixes solution t0 x0 r e
assumes sol: "(solution solves_ode f) T X" and sol_iv: "solution t0 = x0" and iv_defined: "t0 ∈ T"
assumes domain_subset: "X ⊆ X'"
assumes interval: "T = {t0 - e .. t0 + e}"
assumes lipschitz_area: "⋀t. t ∈ T ⟹ cball (solution t) ¦r¦ ⊆ X'"
begin

lemma euler_consistent_solution:
fixes t'
assumes t': "t' ∈ {t0 .. t0 + e}"
shows "consistent solution t' (t0 + e) (euler_C(TYPE('a))) 1 (euler_increment f)"
proof (rule euler_consistent_traj)
show "{t'..t0 + e} ⊆ T" using t' interval by simp
with solves_odeD(1)[OF sol] show "(solution has_vderiv_on (λs. f s (solution s))) {