Existence of Noise Induced Order, a Computer Aided Proof

TL;DR

This paper provides a rigorous computer-aided proof of Noise Induced Order in a chaotic chemical reaction model, demonstrating how increased noise amplitude can induce a transition from chaos to order.

Contribution

It introduces a novel computer-assisted method to rigorously prove the existence of Noise Induced Order in a specific dynamical system, with potential applicability to similar systems.

Findings

Proves the transition from chaos to order as noise increases.

Establishes Lipschitz continuity of the stationary measure under perturbations.

Shows Lyapunov exponent varies H"older continuously with noise amplitude.

Abstract

We prove the existence of Noise Induced Order in the Matsumoto-Tsuda model, where it was originally discovered in 1983 by numerical simulations. This is a model of the famous Belosouv-Zabotinsky reaction, a chaotic chemical reaction, and consists of a one dimensional random dynamical system with additive noise. The simulations showed that an increase in amplitude of the noise causes the Lyapunov exponent to decrease from positive to negative; we give a mathematical proof of the existence of this transition. The method we use relies on some computer aided estimates providing a certified approximation of the stationary measure in the $L^{1}$ norm. This is realized by explicit functional analytic estimates working together with an efficient algorithm. The method is general enough to be adapted to any piecewise differentiable dynamical system on the unit interval with additive noise. We also prove that the stationary measure of the system varies in a Lipschitz way if the system is perturbed and that the Lyapunov exponent of the system varies in a H\"older way when the noise amplitude increases.

Tables

Table 1. Table 1 : The input and output of our computer aided estimates. In particular the first column shows the size of the noise, the fifth column shows the size of the approximation grid, the last column shows intervals enclosing the exact value for the Lyapunov exponent related to the noise size in the first column. See Section 6 for explanations on all the other values.

$\xi$, noise size	${\delta }_{\text{contr}}$	${\alpha }_{\text{contr}}$	$n_{\text{contr}}$	$\delta$	$\alpha$	${\sum }_{i=0}^{n_{\text{contr}}}{C}_i$	a priori $L^1$ err. on measure	${\delta }_{\text{est}}$	refined $L^1$ err. on measure	rig. estimate on Lypunov exponent
$0.860\times {10}^{-2}$	$2^{-18}$	$0.023$	$56$	$2^{-27}$	$0.044$	$29.54$	$0.445\times {10}^{-4}$	$2^{-14}$	$0.246\times {10}^{-5}$	$-{6.03}_{602}^{536}\times {10}^{-1}$
$0.785\times {10}^{-2}$	$2^{-18}$	$0.017$	$58$	$2^{-27}$	$0.039$	$29.20$	$0.479\times {10}^{-4}$	$2^{-14}$	$0.257\times {10}^{-5}$	$-{5.499}_{901}^{221}\times {10}^{-1}$
$0.721\times {10}^{-2}$	$2^{-18}$	$0.013$	$60$	$2^{-27}$	$0.037$	$28.95$	$0.517\times {10}^{-4}$	$2^{-14}$	$0.268\times {10}^{-5}$	$-{4.985}_{749}^{026}\times {10}^{-1}$
$0.661\times {10}^{-2}$	$2^{-18}$	$0.015$	$60$	$2^{-27}$	$0.04$	$28.75$	$0.562\times {10}^{-4}$	$2^{-14}$	$0.282\times {10}^{-5}$	$-{4.467}_{862}^{104}\times {10}^{-1}$
$0.606\times {10}^{-2}$	$2^{-18}$	$0.022$	$54$	$2^{-27}$	$0.05$	$28.44$	$0.612\times {10}^{-4}$	$2^{-14}$	$0.296\times {10}^{-5}$	$-{3.95}_{636}^{558}\times {10}^{-1}$
$0.556\times {10}^{-2}$	$2^{-18}$	$0.015$	$58$	$2^{-27}$	$0.046$	$28.77$	$0.672\times {10}^{-4}$	$2^{-14}$	$0.311\times {10}^{-5}$	$-{3.462}_{884}^{063}\times {10}^{-1}$
$0.509\times {10}^{-2}$	$2^{-18}$	$0.014$	$60$	$2^{-27}$	$0.048$	$29.19$	$0.747\times {10}^{-4}$	$2^{-14}$	$0.330\times {10}^{-5}$	$-{2.97}_{360}^{272}\times {10}^{-1}$
$0.466\times {10}^{-2}$	$2^{-18}$	$0.019$	$59$	$2^{-27}$	$0.056$	$29.56$	$0.832\times {10}^{-4}$	$2^{-14}$	$0.351\times {10}^{-5}$	$-{2.52}_{568}^{476}\times {10}^{-1}$
$0.426\times {10}^{-2}$	$2^{-18}$	$0.032$	$53$	$2^{-27}$	$0.072$	$29.70$	$0.930\times {10}^{-4}$	$2^{-14}$	$0.374\times {10}^{-5}$	$-{2.10}_{419}^{322}\times {10}^{-1}$
$0.391\times {10}^{-2}$	$2^{-18}$	$0.025$	$57$	$2^{-27}$	$0.07$	$30.68$	$0.104\times {10}^{-3}$	$2^{-14}$	$0.400\times {10}^{-5}$	$-{1.74}_{474}^{370}\times {10}^{-1}$
$0.359\times {10}^{-2}$	$2^{-18}$	$0.023$	$60$	$2^{-27}$	$0.074$	$31.73$	$0.118\times {10}^{-3}$	$2^{-14}$	$0.431\times {10}^{-5}$	$-{1.42}_{844}^{733}\times {10}^{-1}$
$0.329\times {10}^{-2}$	$2^{-18}$	$0.028$	$60$	$2^{-27}$	$0.085$	$32.63$	$0.134\times {10}^{-3}$	$2^{-14}$	$0.467\times {10}^{-5}$	$-{1.1}_{600}^{588}\times {10}^{-1}$
$0.302\times {10}^{-2}$	$2^{-18}$	$0.034$	$60$	$2^{-27}$	$0.097$	$33.55$	$0.152\times {10}^{-3}$	$2^{-14}$	$0.510\times {10}^{-5}$	$-9{.}_{411}^{399}\times {10}^{-2}$
$0.277\times {10}^{-2}$	$2^{-18}$	$0.044$	$58$	$2^{-27}$	$0.11$	$34.26$	$0.173\times {10}^{-3}$	$2^{-14}$	$0.560\times {10}^{-5}$	$-7{.}_{706}^{692}\times {10}^{-2}$
$0.252\times {10}^{-2}$	$2^{-18}$	$0.053$	$56$	$2^{-27}$	$0.13$	$35.04$	$0.198\times {10}^{-3}$	$2^{-14}$	$0.626\times {10}^{-5}$	$-{6.3}_{187}^{035}\times {10}^{-2}$
$0.232\times {10}^{-2}$	$2^{-18}$	$0.053$	$56$	$2^{-27}$	$0.14$	$35.94$	$0.223\times {10}^{-3}$	$2^{-14}$	$0.693\times {10}^{-5}$	$-{5.4}_{472}^{305}\times {10}^{-2}$
$0.212\times {10}^{-2}$	$2^{-18}$	$0.055$	$56$	$2^{-27}$	$0.15$	$36.88$	$0.254\times {10}^{-3}$	$2^{-14}$	$0.784\times {10}^{-5}$	$-{4.7}_{957}^{769}\times {10}^{-2}$
$0.192\times {10}^{-2}$	$2^{-18}$	$0.051$	$58$	$2^{-27}$	$0.16$	$38.21$	$0.293\times {10}^{-3}$	$2^{-14}$	$0.909\times {10}^{-5}$	$-{4.3}_{993}^{776}\times {10}^{-2}$
$0.177\times {10}^{-2}$	$2^{-18}$	$0.053$	$58$	$2^{-27}$	$0.18$	$38.98$	$0.329\times {10}^{-3}$	$2^{-14}$	$0.104\times {10}^{-4}$	$-{4.2}_{897}^{654}\times {10}^{-2}$
$0.162\times {10}^{-2}$	$2^{-18}$	$0.054$	$58$	$2^{-27}$	$0.19$	$39.82$	$0.373\times {10}^{-3}$	$2^{-14}$	$0.121\times {10}^{-4}$	$-{4.3}_{479}^{194}\times {10}^{-2}$
$0.150\times {10}^{-2}$	$2^{-18}$	$0.055$	$58$	$2^{-27}$	$0.2$	$40.53$	$0.419\times {10}^{-3}$	$2^{-14}$	$0.139\times {10}^{-4}$	$-4{.}_{521}^{488}\times {10}^{-2}$
$0.137\times {10}^{-2}$	$2^{-18}$	$0.063$	$57$	$2^{-27}$	$0.23$	$41.00$	$0.476\times {10}^{-3}$	$2^{-14}$	$0.164\times {10}^{-4}$	$-{4.7}_{992}^{613}\times {10}^{-2}$
$0.125\times {10}^{-2}$	$2^{-18}$	$0.066$	$58$	$2^{-27}$	$0.25$	$42.14$	$0.554\times {10}^{-3}$	$2^{-14}$	$0.200\times {10}^{-4}$	$-{5.1}_{750}^{293}\times {10}^{-2}$
$0.115\times {10}^{-2}$	$2^{-18}$	$0.071$	$58$	$2^{-27}$	$0.27$	$43.12$	$0.636\times {10}^{-3}$	$2^{-14}$	$0.239\times {10}^{-4}$	$-5{.}_{545}^{491}\times {10}^{-2}$
$0.105\times {10}^{-2}$	$2^{-18}$	$0.079$	$58$	$2^{-27}$	$0.3$	$44.28$	$0.747\times {10}^{-3}$	$2^{-14}$	$0.294\times {10}^{-4}$	$-{5.9}_{835}^{172}\times {10}^{-2}$
$0.960\times {10}^{-3}$	$2^{-18}$	$0.086$	$58$	$2^{-27}$	$0.34$	$45.41$	$0.876\times {10}^{-3}$	$2^{-14}$	$0.360\times {10}^{-4}$	$-6{.}_{427}^{346}\times {10}^{-2}$
$0.885\times {10}^{-3}$	$2^{-18}$	$0.092$	$58$	$2^{-27}$	$0.37$	$46.43$	$0.102\times {10}^{-2}$	$2^{-14}$	$0.436\times {10}^{-4}$	$-6{.}_{855}^{758}\times {10}^{-2}$
$0.810\times {10}^{-3}$	$2^{-18}$	$0.095$	$58$	$2^{-27}$	$0.4$	$47.44$	$0.119\times {10}^{-2}$	$2^{-14}$	$0.534\times {10}^{-4}$	$-7{.}_{320}^{200}\times {10}^{-2}$
$0.748\times {10}^{-3}$	$2^{-18}$	$0.098$	$58$	$2^{-27}$	$0.43$	$48.25$	$0.138\times {10}^{-2}$	$2^{-14}$	$0.643\times {10}^{-4}$	$-7{.}_{711}^{569}\times {10}^{-2}$
$0.686\times {10}^{-3}$	$2^{-18}$	$0.1$	$57$	$2^{-27}$	$0.46$	$48.65$	$0.162\times {10}^{-2}$	$2^{-14}$	$0.779\times {10}^{-4}$	${-}_{8.061}^{7.888}\times {10}^{-2}$
$0.623\times {10}^{-3}$	$2^{-18}$	$0.099$	$56$	$2^{-27}$	$0.5$	$49.04$	$0.192\times {10}^{-2}$	$2^{-14}$	$0.956\times {10}^{-4}$	$-8{.}_{280}^{068}\times {10}^{-2}$
$0.573\times {10}^{-3}$	$2^{-19}$	$0.07$	$60$	$2^{-27}$	$0.29$	$46.78$	$0.141\times {10}^{-2}$	$2^{-15}$	$0.595\times {10}^{-4}$	$-8{.}_{209}^{076}\times {10}^{-2}$
$0.524\times {10}^{-3}$	$2^{-19}$	$0.085$	$56$	$2^{-27}$	$0.33$	$46.32$	$0.162\times {10}^{-2}$	$2^{-15}$	$0.696\times {10}^{-4}$	$-7{.}_{940}^{784}\times {10}^{-2}$
$0.480\times {10}^{-3}$	$2^{-19}$	$0.07$	$58$	$2^{-27}$	$0.34$	$47.75$	$0.185\times {10}^{-2}$	$2^{-15}$	$0.800\times {10}^{-4}$	$-7{.}_{490}^{312}\times {10}^{-2}$
$0.436\times {10}^{-3}$	$2^{-19}$	$0.066$	$58$	$2^{-27}$	$0.36$	$48.66$	$0.215\times {10}^{-2}$	$2^{-15}$	$0.935\times {10}^{-4}$	$-6{.}_{861}^{653}\times {10}^{-2}$
$0.399\times {10}^{-3}$	$2^{-19}$	$0.06$	$60$	$2^{-27}$	$0.39$	$50.39$	$0.254\times {10}^{-2}$	$2^{-15}$	$0.110\times {10}^{-3}$	${-}_{6.208}^{5.964}\times {10}^{-2}$
$0.368\times {10}^{-3}$	$2^{-19}$	$0.065$	$60$	$2^{-27}$	$0.43$	$51.36$	$0.299\times {10}^{-2}$	$2^{-15}$	$0.128\times {10}^{-3}$	$-5{.}_{623}^{339}\times {10}^{-2}$
$0.343\times {10}^{-3}$	$2^{-20}$	$0.041$	$65$	$2^{-27}$	$0.24$	$49.27$	$0.232\times {10}^{-2}$	$2^{-16}$	$0.976\times {10}^{-4}$	${-}_{5.086}^{4.871}\times {10}^{-2}$
$0.312\times {10}^{-3}$	$2^{-20}$	$0.044$	$65$	$2^{-27}$	$0.26$	$50.33$	$0.270\times {10}^{-2}$	$2^{-16}$	$0.112\times {10}^{-3}$	$-4{.}_{436}^{189}\times {10}^{-2}$
$0.287\times {10}^{-3}$	$2^{-20}$	$0.048$	$65$	$2^{-27}$	$0.29$	$51.29$	$0.309\times {10}^{-2}$	$2^{-16}$	$0.127\times {10}^{-3}$	$-3{.}_{844}^{565}\times {10}^{-2}$
$0.259\times {10}^{-3}$	$2^{-20}$	$0.042$	$67$	$2^{-27}$	$0.31$	$53.14$	$0.368\times {10}^{-2}$	$2^{-16}$	$0.149\times {10}^{-3}$	$-2{.}_{984}^{659}\times {10}^{-2}$
$0.237\times {10}^{-3}$	$2^{-20}$	$0.046$	$67$	$2^{-27}$	$0.35$	$54.32$	$0.431\times {10}^{-2}$	$2^{-17}$	$0.144\times {10}^{-3}$	$-1{.}_{998}^{686}\times {10}^{-2}$
$0.218\times {10}^{-3}$	$2^{-20}$	$0.05$	$67$	$2^{-27}$	$0.38$	$55.49$	$0.504\times {10}^{-2}$	$2^{-17}$	$0.165\times {10}^{-3}$	${-}_{9.274}^{5.731}\times {10}^{-3}$
$0.199\times {10}^{-3}$	$2^{-20}$	$0.051$	$68$	$2^{-27}$	$0.42$	$57.28$	$0.606\times {10}^{-2}$	$2^{-17}$	$0.193\times {10}^{-3}$	${}_{2.988}{}^{7.138}\times {10}^{-3}$
$0.184\times {10}^{-3}$	$2^{-20}$	$0.051$	$69$	$2^{-27}$	$0.46$	$59.07$	$0.725\times {10}^{-2}$	$2^{-17}$	$0.225\times {10}^{-3}$	$1{.}_{362}^{844}\times {10}^{-2}$
$0.168\times {10}^{-3}$	$2^{-20}$	$0.053$	$70$	$2^{-27}$	$0.5$	$61.10$	$0.896\times {10}^{-2}$	$2^{-17}$	$0.272\times {10}^{-3}$	$2{.}_{394}^{977}\times {10}^{-2}$
$0.154\times {10}^{-3}$	$2^{-21}$	$0.042$	$74$	$2^{-27}$	$0.29$	$57.78$	$0.650\times {10}^{-2}$	$2^{-18}$	$0.708\times {10}^{-4}$	$3{.}_{533}^{689}\times {10}^{-2}$
$0.142\times {10}^{-3}$	$2^{-21}$	$0.048$	$74$	$2^{-27}$	$0.32$	$59.11$	$0.757\times {10}^{-2}$	$2^{-18}$	$0.822\times {10}^{-4}$	$4{.}_{282}^{462}\times {10}^{-2}$
$0.129\times {10}^{-3}$	$2^{-21}$	$0.049$	$75$	$2^{-27}$	$0.36$	$60.98$	$0.901\times {10}^{-2}$	$2^{-18}$	$0.988\times {10}^{-4}$	$5{.}_{023}^{239}\times {10}^{-2}$
$0.106\times {10}^{-3}$	$2^{-21}$	$0.058$	$75$	$2^{-27}$	$0.45$	$64.34$	$0.135\times {10}^{-1}$	$2^{-18}$	$0.154\times {10}^{-3}$	$6{.}_{626}^{965}\times {10}^{-2}$
$0.873\times {10}^{-4}$	$2^{-21}$	$0.062$	$75$	$2^{-27}$	$0.55$	$67.55$	$0.209\times {10}^{-1}$	$2^{-18}$	$0.252\times {10}^{-3}$	$8{.}_{365}^{917}\times {10}^{-2}$

Equations

$$T_{a,b,c}(x)=\left\{\begin{array}[]{c}(a+(x-\frac{1}{8})^{\frac{1}{3}})e^{-x}+b,~{}~{}~{}0\leq x\leq 0.3\\ c(10xe^{\frac{-10x}{3}})^{19}+b~{}~{}~{}0.3<x\leq 1\end{array}\right.$$

$$a\in\scalebox{0.77}{$0.506073569036822351319599371053047956980141736828203749380990114218225638827_{6}^{9}$},$$

$$b\in\scalebox{0.77}{$0.02328852830307032054478158044023918735669943648088852646123182739831022528_{158}^{213}$},$$

$$c\in\scalebox{0.77}{$0.121205692738975111744666848150620569782497212127938371936404761693002104361_{5}^{8}$}.$$

$${\lambda }_{\xi }:={\int }_0^1\log |T^{\prime }(x)|d{\mu }_{\xi }.$$

$$x\to T(x)\widehat{+}{\omega }_n$$

$$\pi (x)=\underset{i\in \mathbb{Z}}{\min }|x-2i|.$$

$$a\widehat{+}b:=\pi (a+b)$$

$${\rho }_{\xi }(x)=\frac{1}{\xi }\rho (\frac{1}{\xi }x).$$

$$c=\frac{20}{3^{20}\cdot 7}\cdot\bigg{(}\frac{7}{5}\bigg{)}^{1/3}\cdot e^{187/10}.$$

$$a=\frac{19}{42}\cdot\bigg{(}\frac{7}{5}\bigg{)}^{1/3}.$$

$$L:SM(X)\to SM(X)$$

$$L\mu (A)=\mu (T^{-1}(A)).$$

$$N_{\xi }(f):={\rho }_{\xi }\hat{\ast }f$$

$${\rho }_{\xi }\hat{\ast }f={\pi }^{\ast }({\hat{\rho }}_{\xi }\ast \hat{f})$$

$$L_{\xi }:=N_{\xi }L.$$

$$L_{\xi }{\mu }_{\xi }={\mu }_{\xi }$$

$${\lambda }_{\xi }:={\int }_0^1\log |T^{\prime }(x)|d{\mu }_{\xi }.$$

$${\pi }_{\delta }(f)=\mathbf{E}(f|F_{\delta })$$

$$L_{\delta ,\xi }:={\pi }_{\delta }{N}_{\xi }{\pi }_{\delta }L{\pi }_{\delta }.$$

$$L_{\delta ,\xi }^n=({\pi }_{\delta }{N}_{\xi }{\pi }_{\delta }L{)}^n{\pi }_{\delta },$$

$${\tilde{L}}_{\delta ,\xi }={\pi }_{\delta }{N}_{\xi }L{\pi }_{\delta };$$

$$||L_{\delta ,\xi }^{\overline{n}}{|}_V|{|}_{L^1\to L^1}\le \alpha <1$$

$$\Vert f_{\xi }-f_{\xi ,\delta }{\Vert }_{L^1}\le \frac{1}{1-\alpha }{\Vert (L_{\delta ,\xi }^{\overline{n}}-L_{\xi }^{\overline{n}})f_{\xi }\Vert }_{L^1}.$$

$$\Vert f_{\delta ,\xi }-f_{\xi }{\Vert }_{L^1}$$

$$\Vert f_{\xi }-f_{\xi ,\delta }{\Vert }_{L^1}\le \alpha ||f_{\xi }-f_{\xi ,\delta }{\Vert }_{L^1}+{\Vert (L_{\delta ,\xi }^{\overline{n}}-L_{\xi }^{\overline{n}})f_{\xi }\Vert }_{L^1}$$

$${\Vert (L_{\delta ,\xi }^{\overline{n}}-L_{\xi }^{\overline{n}})f_{\xi }\Vert }_{L^1}\le \Vert ({\pi }_{\delta }-1)f_{\xi }{\Vert }_{L^1}+\left(\sum _{i=0}^{\overline{n}-1}{\Vert L_{\delta ,\xi }{|}_V^i\Vert }_{L^1\to L^1}\right)\times$$

$$\times \left({\Vert N_{\xi }({\pi }_{\delta }-1)L{f}_{\xi }\Vert }_{L^1}+{\Vert N_{\xi }{\pi }_{\delta }L({\pi }_{\delta }-1)f_{\xi }\Vert }_{L^1}\right)$$

$$L_{\delta ,\xi }-L_{\xi }$$

$${\Vert (L_{\delta ,\xi }^{\overline{n}}-L_{\xi }^{\overline{n}})f_{\xi }\Vert }_{L^1}$$

$$\le \sum _{i=0}^{\overline{n}}{\Vert ({\pi }_{\delta }{N}_{\xi }{\pi }_{\delta }L{)}^i({\pi }_{\delta }-1)(N_{\xi }L{)}^{n-i}{f}_{\xi }\Vert }_{L^1}$$

Full text

1 Introduction

The “Noise induced order” phenomenon was discovered in numerical simulations and experiments regarding systems modeled by a deterministic dynamics perturbed by noise. The somewhat surprising phenomenon emerging is that the system appears to be chaotic for very small noise intensity, but when the intensity increases the system begins to have a less and less chaotic behavior, passing from a positive Lyapunov exponent to a negative one. A similar behavior was also found for other indicators of chaos. This paper however will focus only on the Lyapunov exponent.

The phenomenon was first discovered by numerical simulations in [23] in a system related to the famous Belousov-Zhabotinsky reaction (see Figure 1) modeled by a one dimensional map perturbed by additive noise. Real experiments confirmed the existence of the phenomenon appearing in the model of the reaction (see [32], and also [31], [33], [34], [20], [10], [24], for some examples of related works).

Despite the impact that the discovery of such noise induced phenomena had in the nonlinear science and physical literature (more than 390 citations to [23] on Google Scholar at the moment of writing this paper) to the best of our knowledge there is no mathematical literature about Noise Induced Order or rigorous proofs of its actual existence in nontrivial systems.

The mathematical study of this phenomenon is difficult because in the deterministic part of the dynamics (see Figure 2) there is a coexistence of strongly expanding and strongly contracting regions and the prevalence of expanding or contracting behavior for typical orbits is a consequence of the global structure of the dynamics.  We remark that this phenomenon is one dimensional and inherently nonlinear, thus mathematically not much related to the noise induced stabilization studied in [2] and following papers. With the help of some computer aided estimates we prove that the global structure of this random system, allows expansion to prevail when the noise amplitude is very small, but the appearance of quite a large noise allows the contraction to prevail.

Our approach is based on the fact that the presence of the noise simplifies the functional analytic properties of the transfer operator associated to the system, smoothing out fine resolution details, and making it well approximable by a finite resolution and finite dimensional one. This makes a computer aided proof possible, letting the computer manage the complexity of the deterministic part of the system at a finite resolution scale and understanding the global structure of the dynamics. However, the computational power required to perform these computations in a naive way is out of the range of current computers. This is true mostly for proving that the Lyapunov exponent is positive in some case of very small noise range. Indeed small noise corresponds to high resolution needed in the study of the system. Because of this we had to find some mathematical clever way to perform the needed estimates, using different functional spaces. This is the main part of the mathematical work contained in what follows and can be applied to many other dynamical systems perturbed by noise. The algorithm we develop in this work was indeed already used in [9] and [3] for the study of other dynamical systems with additive noise considered as models of certain phenomena in climate science and neuroscience.

It is known that naive computer simulations of chaotic systems may not be reliable in some case (see for example [7],[17],[15],[19],[16] for examples of misleading naive simulations and a general discussion on the problem). Beside the pure mathematical interest of a rigorously proved result and a rigorously certified estimate, the study of inherently reliable methods for the numerical study of chaotic dynamical systems is strongly motivated.

Overview of the results. In this work we consider the model of the Belousov-Zhabotinsky reaction studied in [23] (see also [36] for explanations on how the model can be deduced from the chemical mechanism of the BZ reaction). This is a random dynamical system: a deterministic map with additive noise at each iteration. The deterministic part of the dynamics in the model is driven by a map $T_{a,b,c}:[0,1]\rightarrow[0,1]$ defined by

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

The graph of an example of $T_{a,b,c}$ is shown in Figure 2. Following the Inverval Arithmetics formalism we represent intervals in the real line by subscript and superscript describing the decimal expansion of lower and upper bounds for $x$ so that $x\in 0.klm_{tuv}^{xyz}$ means that $x$ belongs to the interval $[0.klmtuv,0.klmxyz]$. The choice of $a,b,c$ follow the one made in [23], adding some more precision (see Remark 4 for more details).

Remark 1**.**

The interval arithmetics and its certified numerical methods (see [29] for an introduction) allow to obtain rigorous results as output of the computer aided estimates. Our computer aided estimates are implemented in this framework. We used SAGE [26] and the validated numerics packages shipped with it. (The interval package is a binding to MPFI [25].) Part of the numerical linear algebra was done using OpenCL [28] running on Nvidia graphic cards.

At each iteration of the map a uniformly distributed noise perturbation with span of size $\xi$ is applied. Further details on the system are presented in Section 2.

In the paper we prove that when the noise size $\xi$ is contained in the interval $[\xi_{1},1/2]$ where $\xi_{1}=\frac{8.73}{10^{5}}$ this random system has a unique ergodic absolutely continuous stationary measure $\mu_{\xi}$ (see Proposition 45) and consider the associated Lyapunov exponent

[TABLE]

We prove that the behaviour of $\lambda_{\xi}$ in the system is similar to the one found by the numerical investigations of Matsumoto and Tsuda ([23]). In particular there is a transition from positive to negative exponent as the noise amplitude increases. We provide explicit examples of values for the noise amplitude having positive and negative Lyapunov exponent. In particular, the findings of the present work applied to the Belousov-Zhabotinsky model defined above allow to state the following

Theorem 2**.**

Let $\lambda_{\xi}$ be the Lyapunov expontent of the system defined above111For every choice of the coefficients $a,b,c$ of $($1$)$ as in $($2$)$,$($3$)$,$($4$)$. with noise of size $\xi$.  For each $\alpha<1,\ \lambda_{\xi}$ is $\alpha-$Hőlder continuous as a function of $\xi$ when $\xi\in[\frac{8.73}{10^{5}},1/2]$; furthermore for $\xi_{1}=\frac{8.73}{10^{5}}$ and $\xi_{2}=\frac{8.60}{10^{3}}$ it holds

• I1

the Lyapunov exponent $\lambda_{\xi_{1}}\in[\frac{8.365}{10^{2}},\frac{8.917}{10^{2}}]$, hence it is rigorously certified to be positive;

• I2

the Lyapunov exponent $\lambda_{\xi_{2}}\in[\frac{-6.03602}{10},\frac{-6.03536}{10}]$, hence it is rigorously certified to be negative.

Therefore, the system exibits Noise Induced Order.

We refer to Section 6 for the certified results proving the change of sign for the Lyapunov exponent, while we refer to Section 7 (see Corollary 42) for the Hőlder regularity of the Lyapunov exponent. In the paper, we estimate similarly the Lyapunov exponent for many other values of the noise size. The results are reported in Table 6 and described in Figure 4 by a graph.

The method of proof and certification of these results relies on the approximation of the stationary measure $\mu_{\xi}$ of the system up to a certified error in the $L^{1}$ norm222The computation of stationary measures for dynamical systems perturbed by noise was approached from different points of view in [18] where the convergence (without effective bounds on the approximation error) of an approximation scheme based on Fourier analysis was proved for certain classes of maps, and in [6], where the computability of the stationary measure up to a given error was considered in an abstract framework, giving bounds on the computational complexity of the problem. The rigorous computation of the stationary measure for expanding and contracting Iterated Function Systems is considered in [14].. Using this approximation we can compute the Lyapunov exponent of the system up to a small certified error. The methods used however are general enough to be applied to any system formed by piecewise differentiable maps of the interval perturbed by additive noise having a Bounded Variation distribution (see Sections 2 and 3.3 for a more detailed discussion).

Plan of the paper. In Section 2 we describe the systems under study and we introduce the technique we use for the approximation of the transfer operator. The techniques leading to the proof of Theorems 2 are explained in a sequence of settings of decreasing generality where the needed assumptions on the system are listed (see Settings 13, 18). In Section 3 we describe how to find an explicit bound on the approximation error for the computation of the stationary measure. In Subsection 3.3 we describe an efficient procedure which exploits information coming from a coarse knowledge of the stationary measure in a bootstrap argument. This procedure uses several technical lemmas estimating the variation of certain densities and the norms of certain operators; such lemmas are listed in Section 3 and proved in Section 9 to make reading easier. In Section 4 we show an efficient procedure (which is used in the main algorithm) for the estimation of the rate of contraction of a finite rank transfer operator, when applied to zero average measures. This is a quantitative measure of the rate of convergence to equilibrium of the system which is involved in the quantitative estimate of its stability under perturbation. This is important in the estimation of the approximation error. In Section 5 we show the estimates needed to compute the Lyapunov exponent of the system once we know the stationary measure. In Section 6 we apply all these techniques to the system described in Section 2, showing the results of our computer aided estimates, and proving the existence of noise induced order (and in particular, Theorem 2). In Section 7 we consider the stability of the stationary measure to changes in the system’s parameters and of the Lyapunov exponent on changes of the noise amplitude proving the Hőlder continuity. Section 8 contains some definitions and generalities about Random Dynamical Systems we include for completeness and to justify the correctness of the notion of Lyapunov exponent which is estimated in the paper.

Aknowledgements.

The authors thank Yuzuru Sato for the nice seminar where we learned about the existence of Noise Induced Order and further explanations, EU Marie-Curie IRSES “Brazilian-European partnership in Dynamical Systems” (FP7-PEOPLE-2012-IRSES 318999 BREUDS) for support during the research and The Abdus Salam International Centre for Theoretical Physics (ICTP) where the work started. I.N. was partially supported by CNPq, University of Uppsala and KAW grant 2013.0315. I. N. thanks UFRJ, CAPES (through the programs PROEX and the CAPES-STINT project ”Contemporary topics in non uniformly hyperbolic dynamics”). S.G. thanks Leverhulm Trust grant (IN-2014-021, “Statistical properties of non uniformly hyperbolic dynamical systems: computer assisted proofs and rigorous computation”) and GNAMPA/INDAM (“Stabilita e instabilita in sistemi con rumore, un approccio computer assistito.”)

2 The system and its transfer operator

In this section we describe more precisely the system which will be studied in the paper and the associated transfer operators. Basic notions on random dynamical systems, stationary measures and Lyapunov exponents we use in the paper are presented in Section 8.

A random dynamical system with additive noise on $[0,1]$ and reflecting boundary conditions is a random perturbation of a deterministic map, defined by

[TABLE]

where $T:[0,1]\rightarrow[0,1]$ is a Borel measurable map and $\omega_{n}$ is an i.i.d. process distributed according to a probability density $\rho_{\xi}$ and $\hat{+}$ is the ”reflecting boundaries sum” on $[0,1]$ defined as follows.

Definition 3**.**

Let $\pi:{\mathbb{R}}\rightarrow[0,1]$ be the piecewise linear map

[TABLE]

Let $a,b\in{\mathbb{R}}$ then

[TABLE]

where $+$ is the usual sum operator on $\mathbb{R}$. By this $a\hat{+}b\in[0,1].$

In the following we will consider the case where the noise density $\rho_{\xi}$ is the rescaling of some Bounded Variation kernel $\rho\in BV[-\frac{1}{2},\frac{1}{2}]$ with $\int\rho=1$ in the interval $[-\frac{\xi}{2},\frac{\xi}{2}]$, hence

[TABLE]

(see Definition 19 for a recall on the definition of bounded variation.)

As described in the introduction, the model studied in the present paper is the one studied in [23]. We consider a random dynamical system with additive noise, as in (5) where $T=T_{a,b,c}$ is defined in (1). At each iteration of the map a uniformly distributed noise perturbation with span of size $\xi$ with reflecting boundaries is applied.

Remark 4**.**

The parameters $a,b,c$ defined below (1) have been computed using interval arithmetic, in the implementation of our algorithm they are represented by intervals. We explain the motivation for the choice of the parameters in [23]. The parameter $c$ is defined in a way to be nearby to a value for which $T(0.3^{-})=T(0.3^{+})$, making $T$ continuous at $0.3$. The exact value of $c$ giving the continuity can be computed in a closed form as:

[TABLE]

The parameter $a$ is defined similarly in a way so that $T^{\prime}(0.3^{-})=0$, making $T^{\prime}$ continuous at $0.3$. The value of such $a$ can be computed in a closed form as:

[TABLE]

The choice of the parameter $b$ in [23] is motivated by a parallel with the logistic map. Let us denote by $T_{b}$ the map as only the parameter $b$ varies; each $T_{b}$ has a repelling fixed point $p_{b}$. In [23] this explicit value of $b$ is chosen as an approximation of the parameter value for which $T^{4}(0.3)=p_{b}$ following the kneading sequence “RLLL”, i.e., a Misiurewicz condition. We computed a certified interval containing $b$ using a Newton Interval Method [29]. The interval enclosing $b$ is computed giving the result shown at (3).

In our computer assisted estimates we will hence consider the map given at (1) and uniform noise, however the mathematical treatment about approximation of stationary measures in Section 3 and following is more general. We start considering a general random dynamical system where $T$ is measurable and $\rho\in BV$ giving general results and estimates we then improve using more assumptions on the system (i.e. $T$ piecewise smooth) putting ourselves in different general settings which are stated precisely (see Settings 13, 18), to keep the exposition as clear and general as possible.

Remark 5**.**

The reflecting boundary condition at (5) is not influent when the noise amplitude is smaller than the parameter $b$ (as it will be for all the noise amplitudes considered in our computer aided proofs, see Table 6).

**The transfer operator. **We study the statistical properties of the dynamical systems with additive noise, as defined at $($5$)$ through the study of the properties of their associated transfer operators. Let us recall that a measurable map $T:X\rightarrow X$, induces a map

[TABLE]

where $SM(X)$ is the space of Borel signed measures on $X.$ The associated map $L$ is defined in the following way: if $\mu\in SM(X)$ then:

[TABLE]

In the literature, $L$ is also called the pushforward map associated to $T$, sometime denoted by $T^{\ast}$. It is a linear operator on the vector space $SM(X)$ and it is also called the transfer operator associated to $T$. The space of measures with density in $L^{1}([0,1])$ can be seen as a subspace of $SM([0,1])$. If $T$ is nonsingular, $L$ can be considered as an operator $L^{1}([0,1])\rightarrow L^{1}([0,1])$.

The (annealed) transfer operator $L_{\xi}$ associated to the system with noise is given by the composition of the transfer operator $L$ and a reflecting boundary convolution operator $N_{\xi}:SM([0,1])\rightarrow L^{1}([0,1])$ (a suitable modification of the usual convolution), defined by

[TABLE]

where the “reflecting boundaries convolution” $\hat{\ast}$ is defined similarly to the reflecting boundaries sum as

Definition 6**.**

Let $\mu\in SM({\mathbb{R}})$. Let $\pi:{\ \mathbb{R}}\rightarrow[0,1]$ be the piecewise linear map defined at (6) and $\pi^{\ast}$ its associated pushforward map. We consider $\pi^{\ast}\mu\in SM([0,1])$ as the “reflecting boundary” version of $\mu$.

Definition 7**.**

Let $f\in SM([0,1])$, $\rho_{\xi}\in BV[-{\xi},{\xi}]$. Let $\hat{f}\in SM(\mathbb{R})$ defined by $\hat{f}=1_{[0,1]}f$ and $\hat{\rho}_{\xi}\in L^{1}(\mathbb{R})$ by $\hat{\rho}_{\xi}=1_{[-{\xi},{\xi}]}\rho_{\xi}$. We define

[TABLE]

where $\ast$ stands for the usual convolution operator on $\mathbb{R}$.

This boundary reflecting convolution operator is regularizing, in particular $\rho_{\xi}\hat{\ast}f\in BV([0,1])$ if $f\in BV$, and has properties similar to the usual convolution operator. For its basic properties see Subsection 3.3.1.

Definition 8**.**

The annealed transfer operator $L_{\xi}:SM([0,1])\rightarrow L^{1}([0,1])$ associated to a deterministic system with additive noise, as described at (5) is defined as

[TABLE]

Remark 9**.**

The annealed transfer operator is obtained by averaging the transfer operator $L$ over all the possible noise perturbations. We refer to Section 8 for some basic facts on this operator. In the following we will mainly consider $L_{\xi}$ as an operator $L^{1}([0,1])\rightarrow L^{1}([0,1])$.  In the notation we emphasize the dependence of the operator on the amplitude of the noise.

Definition 10**.**

Let $\mu_{\xi}\in L^{1}([0,1])$ be a fixed probability measure for $L_{\xi}$, i.e,

[TABLE]

we will call $\mu_{\xi}$ a stationary measure for the system with additive noise or an invariant measure for $L_{\xi}$.

By the regularizing properties of the convolution by a Bounded Variation kernel, and standard compactness arguments it is easy to see that the transfer operator $L_{\xi}$ corresponding to a map with additive noise has at least one fixed point $f_{\xi}$ in $BV[0,1]$ (see [11], Lemma 23 for more details). Following [23] we will study the Lyapunov exponent of the system as $\xi$ varies. The Lyapunov exponent is defined as follows:

Definition 11**.**

The average Lyapunov exponent associated to noise size $\xi$ is

[TABLE]

When $\mu_{\xi}$ is ergodic the average Lyapunov exponent coincides almost everywhere with the pointwise Lyapunov exponent (see Section 8). As a byproduct of our computer aided estimates, using the results given in Section 4 and 7 we will prove that the systems considered are ergodic (see Proposition 45). This justifies the correctedness of the average Lyapunov exponent as an indicator of the behaviour of the system.

The Ulam Approximation. The main tool for the study of the behavior of (10) in this work is the rigorous approximation of the stationary measure $\mu_{\xi}$. This is done by approximating $L_{\xi}$ by a finite dimensional operator $L_{\delta,\xi}:L^{1}([0,1])\rightarrow L^{1}([0,1])$. The fixed points of $L_{\xi}$ are then approximated by the ones of $L_{\delta,\xi}$ with a certified bound on the approximation error.

Let $\pi_{\delta}:L^{1}([0,1])\rightarrow L^{1}([0,1])$ be a projection on a finite dimensional space defined in the following way: the space $[0,1]$ is discretized by a partition $I_{\delta}$ (with $k$ elements); the projection considered is defined by the conditional expectation

[TABLE]

where $F_{\delta}$ is the $\sigma-$algebra associated to the partition $I_{\delta}.$

The approximated operator is then defined by finite element approach, composing with $\pi$:

[TABLE]

The finite dimensional approximation of an operator based on the conditional expectation, as above, is commonly called Ulam discretization or Ulam method. This method was widely studied in the literature (see, e.g [8],[21],[4],[5],[12]).

Observe that

[TABLE]

taking into account that $\pi_{\delta}^{2}=\pi_{\delta}$. We remark that since $||\pi_{\delta}||_{L^{1}\rightarrow L^{1}}\leq 1$ and $||L||_{L^{1}\rightarrow L^{1}}\leq 1$, then $||L_{\delta,\xi}||_{L^{1}\rightarrow L^{1}}\leq 1$.

Remark 12**.**

Another discretization that could be used is

[TABLE]

*while this definition is reasonable, it is more difficult to implement and would force us to recompute the discretized operator for each size $\xi$ of the noise. Our definition permits us to compute once and for all $\pi_{\delta}L\pi_{\delta}$ which is computationally expensive but leads to a sparse matrix, and then apply the operator $\pi_{\delta}N_{\xi}\pi_{\delta}$ which is independent of the dynamics. *

3 Rigorous approximation of the stationary measure for dynamical systems with additive noise

The certified approximation of the Lyapunov exponent is based on the certified approximation of the stationary measure of the system in the $L^{1}$ norm. This is the main part of our general construction and is described in this section. The algorithm to approximate the stationary measure uses both a priori (analytical) and a posteriori (computer assisted) estimates on the measure and on the transfer operator.* For these estimates to be performed we do not need particular expansion or hyperbolicity properties of the one dimensional map driving the deterministic part of the dynamics.* We now introduce a general and simple algorithm which works for measurable maps perturbed by additive noise with a Bounded Variation kernel, in Section 3.3 we refine this algorithm, assuming that the map is piecewise smooth and getting much better estimates. In this first part of the section we will hence work in the following setting.

Setting 13**.**

Let us suppose $L_{\xi}$ is the transfer operator of a system with additive noise as considered in (5). We suppose the noise is distributed according to a Bounded Variation kernel $\rho_{\xi}$ with support in $[-\frac{\xi}{2},\frac{\xi}{2}]$ and the deterministic part of the system is driven by a measurable map $T$.

Let $L_{\delta,\xi}$ be the Ulam approximation of $L_{\xi}$, defined by projecting on a partition of size $\delta$. Let $f_{\delta,\xi},$ $f_{\xi}\in L^{1}$ respectively be invariant probability measures of $L_{\delta,\xi}$ and $L_{\xi}$. Since the measure $f_{\delta,\xi}$ is a fixed point of the finite dimensional operator $L_{\delta,\xi}$, it can be computed to any precision. We will treat now the issue of estimating $\|f_{\delta,\xi}-f_{\xi}\|_{L^{1}}$ effectively.

Lemma 14**.**

Suppose that for some $\overline{n}\in\mathbb{N}$

[TABLE]

where $V=\{\nu\in L^{1},\nu([0,1])=0\}$ being the space of zero-average measures. Then

[TABLE]

Remark 15**.**

We remark that $L_{\delta,\xi}|_{V}$ is a finite dimensional operator and can be represented by a matrix. Thus is possible for a computer to verify that $||L_{\delta,\xi}^{\overline{n}}|_{V}||_{L^{1}\rightarrow L^{1}}\leq\alpha$ for some $\overline{n}$.

Proof.

(of Lemma 14)  Since both $f_{\xi}$, $f_{\delta,\xi}$ are fixed points we can write

[TABLE]

Then

[TABLE]

implying the statement.

3.1 An informal description of the main algorithm to compute the

stationary measure up to a small given error.

Based on Lemma 14, a strategy to rigorously bound $\|f_{\delta,\xi}-f_{\xi}\|_{L^{1}}$ is the following. The computer will find an $\overline{n}$ such that (12) is satisfied, then (13) will give an estimate for the approximation error. We remark that if $\delta$ is small enough, $\left\|(L_{\delta,\xi}^{\overline{n}}-L_{\xi}^{\overline{n}})f_{\xi}\right\|_{L^{1}}$ has a chance of being small, since it is the difference of two nearby operators, both applied to the same regular (Bounded Variation) measure. This is where the size $\delta$ of the approximation grid has a role in the quality of the approximation. On the other hand, $\overline{n}$ may depend on $\delta$. This is why a priori estimates on the approximation error are not trivial. Lemma 14  provides some a posteriori estimate on the error; the approximation error is known after the computer certifies the $\overline{n}$ and the $\alpha$ for which (12)  is satisfied.

Hence the main algorithm for the approximation of the invariant measure will work as follows:

Given the grid size $\delta$, compute $L_{\delta,\xi}$ and $f_{\xi,\delta}$ up to some prescribed precision. 2. 2.

Find a good $\overline{n}$ and $\alpha$: we estimate $||L_{\delta,\xi}^{n}|_{V}||_{L^{1}\rightarrow L^{1}}$ in an efficient way, finding a good compromise between $\overline{n}$ and $\alpha$, in a way that $\alpha$ is not too close to $1$ and $\overline{n}$ not too big. We remark that the norm of the finite dimensional operator $L_{\delta,\xi}^{n}$ is directly computable in principle, but if $\delta$ is small then the size of the associated matrix is huge and this can be a hard computational task. For this we use a coarse-fine strategy which is explained in Section 4, and which takes into account that the huge matrix representing $L_{\delta,\xi}$ is actually coming from a certain dynamical system with noise. 3. 3.

Find a good estimate for $\left\|(L_{\delta,\xi}^{\overline{n}}-L_{\xi}^{\overline{n}})f_{\xi}\right\|_{L^{1}}$. We remark that in this formula $f_{\xi}$ is not known, but still we can find enough information on it to estimate the difference of operators we are interested in. This will be done by a method using both a priori and a posteriori estimates, using an approximated knowledge of $f_{\xi}$ and its variation. The procedure is explained in Section 3.2 and in the following sections. A simple but not efficient bound (equation (16)) is proved in 3.2.1; we refine the method in Section 3.3 greatly improving the efficency of the estimate with a bootstrap argument. 4. 4.

Estimate the approximation error $\|f_{\xi,\delta}-f_{\xi}\|_{L^{1}}$ by Lemma 14.

3.2 A bound for $\left\|(L_{\delta,\xi}^{\overline{n}}-L_{\xi}^{\overline{n}})f_{\xi}\right\|_{L^{1}}$

Having outlined the main algorithm we now show how to perform the main needed estimates. In this Section we describe how to estimate the quantity appearing on the right hand side of (13). This will be done by splitting this term in different parts which will be treated differently. We start estimating the term as a telescopic sum whose summands will be estimated in the following subsections.

Lemma 16**.**

Let $L_{\xi}$ the transfer operator of the random system and $L_{\delta,\xi}$ its Ulam approximation, as defined in Section 2. Let $f_{\xi}$ be an invariant probability measure for $L_{\xi}$, it holds

[TABLE]

Proof.

The proof is based on a telescopic decomposition. We have

[TABLE]

Performing a similar decomposition to $L_{\delta,\xi}^{\overline{n}}=(\pi_{\delta}N_{\xi}\pi_{\delta}L)^{\overline{n}}\pi_{\delta}$. Pairing in a suitable way the terms we inserted we obtain:

[TABLE]

considering that $f_{\xi}$ is fixed by $N_{\xi}L=L_{\xi}$. Shifting indexes by $1$ in the the first sum, the estimate can be written as

[TABLE]

(notice that $(\pi_{\delta}-1)g$ has always average zero for any $g$, and consequently belongs to $V$).

3.2.1 An initial (a priori) bound for $\|f_{\xi}-f_{\xi,\delta}\|_{L^{1}}$

Now we show a strategy to get a simple effective bound for the approximation error $\|f_{\xi}-f_{\xi,\delta}\|_{L^{1}}$ based Lemma 16, estimating the summands on the right hand side of (14) by quantities which are known from the description of the system or can be computed by its approximated transfer operator. In the next section we will improve the method, using more information on $T$ and $f_{\xi}$ and getting much more efficient estimates.

Lemma 17**.**

Let $f_{\xi}$ a stationary measure for a system defined as in (5) and $f_{\xi,\delta}$ a stationary measure for its Ulam approximation, as defined in Section 2. If there is $\overline{n}$ such that

[TABLE]

then

[TABLE]

where $0\leq C_{i}\leq 1$ are such that $\|L_{\delta,\xi}^{i}|_{V}\|_{L^{1}\rightarrow L^{1}}\leq C_{i}$.

Proof.

The proof of the lemma is based on the following estimates, proved in Corollary 48 and Prop. 49 (Section 9) allowing a first estimate on the right hand side of (13). We have

[TABLE]

[TABLE]

Now we apply (17,18) to the summands of the righ hand side of (14). We see that all the items there have either a $N_{\xi}(1-\pi_{\delta})$ or a $(1-\pi_{\delta})N_{\xi}$ appearing. Indeed since $||f_{\xi}||_{1}\leq 1$

[TABLE]

Similarly the other summands, can be estimated recalling that $\|\pi_{\delta}\|_{L^{1}\rightarrow L^{1}}\leq 1$ and $\|N_{\xi}\|_{L^{1}\rightarrow L^{1}}\leq 1$. Applying (13) we get the statement.

The estimate given at (16) mainly depend on the ratio $\delta\xi^{-1}$ between the partition size and the noise amplitude. This estimate is obtained without any information on the deterministic part of the dynamics, only the information about the contraction rate of the approximated transfer operator $L_{\delta,\xi}$ (to obtain $\overline{n}$ and $\alpha$). This would already allow to obtain a good approximation for the invariant density $f_{\xi}$, in principle, if we had enough computation power to carry on the computation with a very small $\delta$. Unfortunately, in the Matsumoto-Tsuda system, positive Lyapunov exponent appears for very small sizes of the noise making the computation unfeasible even with the help of a supercomputer, due to the growth of the computational complexity as $\xi$ becomes small.

Therefore, we have to apply a more subtle and complicated strategy where the bootstrap argument comes into play, i.e., using some information on $f_{\xi}$ (and in particular about its variation in given intervals) which we can obtain with a preliminary computation.

3.3 A stronger (a posteriori) bound

In this section we analyze better (14) and see how using some more assumptions on $T$, more information on $f_{\xi}$  and the use of the Wasserstein distance, we can drastically improve the estimates given in Lemma 17. We remark that the explicit error bound provided by Lemma 17 is proportional to $\delta\xi^{-1};$ the new error estimate will be a sum where most summands are proportional to $\delta^{2}\xi^{-1}$.

Setting 18**.**

From now on we will suppose we are in the framework of Setting 13 and furthermore we suppose $T$ being piecewise smooth. We suppose that there is a partition $\{P_{i}\}_{1\leq i\leq k}$ such that

• •

each $P_{i}$ is an interval,

• •

on each $P_{i}$ the branch $T_{i}:=T|_{P_{i}}$ is monotonic and $C^{2}$ in the interior of $P_{i}$

• •

The limits of $T_{i}^{\prime}(x)$ as $x$ tend to the frontier of $P_{i}$ exist in $\mathbb{R}\cup\{-\infty,\infty\}$.

We will consider the transfer operator $L$ related to the map and to every branch of the map. For each $L^{1}$ density $g$, we let $L_{i}g$ be the component of $Lg$ coming from the $i$-th monotone branch (the pushforward map related to $T_{i}$), that is

[TABLE]

In this way we have $Lg=\sum_{i}L_{i}g$.

Once the discretized transfer operator $L_{\delta,\xi}$ is computed and has a unique fixed probability measure $f_{\delta,\xi},$ an approximation $\tilde{f}$  for $f_{\delta,\xi}$ can be computed up to any given small error in $L^{1}$ (this is the computation of the fixed point of the big matrix representing $L_{\delta,\xi}$). Let us assume that a numerical approximation $\tilde{f}$ of $f_{\delta,\xi}$ in the $L^{1}$ norm is computed. Using $\tilde{f}$, we will look for a strategy for estimating the total error $\|f_{\xi}-\tilde{f}\|_{L^{1}}$, that includes the appoximation error $\|f_{\xi}-f_{\delta,\xi}\|_{L^{1}}$ (because we approximated on a partition of size $\delta$) and the numerical error $\|f_{\delta,\xi}-\tilde{f}\|_{L^{1}}$. We remark that

[TABLE]

For the estimation of $\|f_{\delta,\xi}-\tilde{f}\|_{L^{1}}$ in our algorithm we apply the same method described in [12]. To find a stronger bound for $\|f_{\xi}-f_{\xi,\delta}\|_{L^{1}}$ let us start again from the estimate given at Lemma 16. We will estimate independently the terms

[TABLE]

appearing at (14), using the information we can extract from the approximation $\tilde{f}$. Except in the first case (that has the smallest weight in the estimate, according to (14)), the estimate will become roughly proportional to $\delta^{2}\xi^{-1}$, greatly improving the quality of the approximation certification. For each of the terms in (21) we prove in Sections 3.3.2, 3.3.3, 3.3.4 bounds of the form

[TABLE]

where $A$ and $B$ depend on $\delta$ and become small when $\delta$ is small. We remark that these bounds depend on the error $||\tilde{f}-f_{\xi}||_{L^{1}}$ itself. This together with precise bounds, based on some approximation $\tilde{f}$, permits us to tighten the bounds on the error $||\tilde{f}-f_{\xi}||_{L^{1}}$, using a so called “bootstrapping” process.

3.3.1 A summary of norms and estimates

Before entering in the details of the estimate of 21, we introduce some of the norms used in the paper, and we summarize some of the bounds that are used in this section and proved in Appendix 9.

Definition 19**.**

Let $X\subset[0,1]$ be a finite union of pairwise disjoint intervals, $X=\bigcup_{j}I_{j}$, $I_{i}\cap I_{j}=\emptyset$ for $i\neq j$. We define the variation on $X$ of the function $f$ (and denote it by $\mathrm{Var}_{X}(f)$) as follows:

• •

when $X$ is an interval (the endpoints may be included or not), the variation is defined as

[TABLE]

(the supremum being over finite increasing sequences of any length $k$ contained in $X$);

• •

when $X$ is a finite union of pairwise disjoint intervals $X=\bigcup_{j}I_{j}$, $I_{i}\cap I_{j}=\emptyset$ for $i\neq j$, the variation is defined as $\mathrm{Var}_{X}(f):=\sum_{j}\mathrm{Var}_{I_{j}}(f)$.

Remark 20**.**

As it is well known, if $f$ has a weak derivative, then $\mathrm{Var}_{X}(f)=\|f^{\prime}\|_{L^{1}(X)}$, see for instance [27, Chap. 4, Prop. 4.2].

We will also consider a norm which is weaker than the $L^{1}$ norm.

Definition 21**.**

Let $f$ be a function in $L^{1}([0,1])$ with zero average, we define the Wasserstein-like norm of $f$, as

[TABLE]

Let $Y\subset[0,1]$ be a finite union of pairwise disjoint intervals $Y=\bigcup_{j}I_{j}$, $I_{i}\cap I_{j}=\emptyset$ for $i\neq j$, we denote by $W(Y)$ the space of  $L^{1}$ functions with support contained in $Y$ having zero average in each $I_{j}$. If $f\in W(Y)$ we define its norm by

[TABLE]

Next proposition contains a summary of the bounds used in next proofs and the location of their proof in the paper; it can be used as an handy guideline throughout the paper, we refer to the cited lemmas and proposition for details.

Proposition 22** (Summary of the bounds).**

Let $\pi_{\delta}$ be the Ulam projection on a homogeneous partition of size $\delta$, as defined in (11), let $I$ be a set which is a finite union of intervals of the partition then:

$||1-\pi_{\delta}||_{\mathrm{Var}\to L^{1}}\leq\delta/2$, Lemma 47, 2. 2.

$||1-\pi_{\delta}||_{L^{1}\to W}\leq\delta/2$, Lemma 51, 3. 3.

$||1-\pi_{\delta}||_{\mathrm{Var}(I)\to W(I)}\leq\delta^{2}/8$, Lemma 52.

Let $N_{\xi}$ be the convolution operator then:

$||N_{\xi}||_{L^{1}\to\mathrm{Var}}\leq\xi^{-1}\mathrm{Var}(\rho)$, Lemma 46, 2. 5.

$||N_{\xi}||_{W\to L^{1}}\leq\xi^{-1}\mathrm{Var}(\rho)$, Lemma 50, 3. 6.

$||(1-\pi_{\delta})N_{\xi}||_{L^{1}\to L^{1}}\leq\frac{1}{2}\delta\xi^{-1}\mathrm{Var}(\rho)$, Corollary 48, 4. 7.

$||N_{\xi}(1-\pi_{\delta})||_{L^{1}\to L^{1}}\leq\frac{1}{2}\delta\xi^{-1}\mathrm{Var}(\rho)$, Proposition 49.

Let $L$ be the transfer operator associated to $T$, and let $L_{i}g$ the component of $Lg$ coming from the $i-th$ branch as defined at (19), then:

$||L||_{W(I)\to W}\leq||T^{\prime}||_{L^{\infty}(I)}$, Lemma 54 2. 9.

in Lemma 55 the local variation inequality is proved:

[TABLE]

3.3.2 Estimate for $\|(\pi_{\delta}-1)f_{\xi}\|_{L^{1}}$

We give now an estimate for the first item of (21).

Lemma 23**.**

Let $\pi_{\delta}$ be the Ulam projection on a homogeneous partition of size $\delta$, then

[TABLE]

for

[TABLE]

Proof.

We estimate:

[TABLE]

where in the last line we used Proposition 22, items 6) and 1).

An upper bound on $\mathrm{Var}(N_{\xi}L\tilde{f})$ will be estimated using the results in Section 9.2.2 and the explicit knowledge of the computed $\tilde{f}$.

3.3.3 Estimate for $||N_{\xi}(1-\pi_{\delta})Lf_{\xi}||_{L^{1}}$

We give now the estimate of the second item of (21). The main idea is to use a coarser partition $\Pi$ made of intervals whose size is an integral multiple of $\delta$. Then, instead of estimating $||N_{\xi}(1-\pi_{\delta})Lf_{\xi}||_{L^{1}}$ globally we bound this quantity on each interval of $\Pi$ exploiting the approximate knowledge of the variation of $f_{\xi}$ in the interval. This drastically improves the quality of the estimate in almost every interval of $\Pi$ and allows the error-checking computation to be performed on a partition which is coarser than the initial partition of size $\delta$.

Lemma 24**.**

Let $\Pi$ be a uniform partition whose parts have size that is an integral multiple of $\delta$, $T$ piecewise monotonic with $L_{i}$ defined as above. We have

[TABLE]

with

[TABLE]

Proof.

We can estimate as

[TABLE]

The first summand can be rewritten as

[TABLE]

splitting on the intervals of the partition, and using Lemma 55 (because each $I\in\Pi$ is a union of intervals of the partition of size $\delta$)

[TABLE]

In the last step we used the property of the $W$ norm stated in Lemma 51 and 52. Adding the second term of (28) we have proved the result.

Therefore, once we have an approximation $\tilde{f}$ we can compute $A_{2}$, $B_{2}$ by a simple algorithm that evaluates the double summation in the last equation of the above Lemma.

Remark 25**.**

When computing $B_{2}$, the minimum could be always the one depending on $\|L_{i}{\tilde{f}}\|_{L^{1}(I)}$. If this happens, the sum adds up to $1$, and the new bound is worse than the a priori estimate, since we are introducing a factor $\|f_{\xi}-\tilde{f}\|_{L^{1}}+1>1$.

In practice, this does not happen. On each $i$-th preimage of an interval $I$ of the partition the new estimate will provide a better bound as soon as

[TABLE]

First of all, the variation of $\tilde{f}$ has an a priori bound by

[TABLE]

Now, we can try to control $\mathrm{Var}_{I}(L_{i}\tilde{f})$ by using the local variation inequality in Proposition 22. If the preimage $T_{i}^{-1}(I)$ does not contain a critical point or a singular point, we expect $(\delta^{2}/8)\cdot\mathrm{Var}_{I}(L_{i}\tilde{f})$ to be small.

For the intervals $I$ where we cannot apply the local variation inequality or where it does not give us good enough bounds we fall back to the a-priori estimate depending on the $L^{1}$ mass rather than on the variation.

3.3.4 An estimate for $||N_{\xi}\pi_{\delta}L(\pi_{\delta}-1)f_{\xi}||_{1}$

We give an estimate of the third item of (21). The general idea is similar to the one explained in the previous section, again, some required estimates are technical lemmas proved in Section 9.

Lemma 26**.**

Let $\Pi$ be a uniform partition whose parts have size that is multiple of $\delta$, we have:

[TABLE]

with

[TABLE]

Proof.

We have

[TABLE]

An algorithm for estimating for $\mathrm{Var}(N_{\xi}L\tilde{f})$ can be found in Lemma 56. The term $\|N_{\xi}L(1-\pi_{\delta})N_{\xi}L\tilde{f}\|_{L^{1}}$ can be split over the intervals $I\in\Pi$ and estimated as

[TABLE]

proving the statement thanks to Lemma 54 (because each $I$ is a union of intervals of the partition of size $\delta$).

Remark 27**.**

To estimate $B_{3}$, we estimate computationally $\mathrm{Var}_{I}(N_{\xi}L\tilde{f})$ for each interval $I\in\Pi$ using the algorithms explained in Sections 9.2.1 and 9.2.2. As in Lemma 24, we obtain a stronger estimate in the interval $I$ as soon as

[TABLE]

Remark that in all our computations ${\delta\xi^{-1}\mathrm{Var}(\rho)}$ needs to be small, since it controls the approximation error (refer to Proposition 22, items 6 and 7). This implies that the inequality above will be true for most of the intervals but those where $T^{\prime}$ becomes very big.

3.3.5 An estimate for the $L^{1}$ error $\|f_{\xi}-\tilde{f}\|_{L^{1}}$

In the previous sections we built the ingredients for estimating $\|f_{\xi}-f_{\xi,\delta}\|_{L^{1}}$, but we want to estimate $\|f_{\xi}-\tilde{f}\|_{L^{1}}$ where $\tilde{f}$ is the output of a computation approximating the fixed point of $L_{\delta,\xi}$. We will do so assuming that we have an estimate of the numerical error $\|f_{\xi,\delta}-\tilde{f}\|_{L^{1}}$ (such an estimate can be found in [12]).

Let $A_{i}$, $B_{i}$ ($i=1,2,3$) be the constants defined as in (25), (24), (30). Plugging (24), (26), (29) (according to Lemmas 24 and 26) into (14), we have that $\|(L_{\delta,\xi}^{\overline{n}}-L_{\xi}^{\overline{n}})f_{\xi}\|_{L^{1}}$ can be bounded as

[TABLE]

where

[TABLE]

Thanks to (13) we have

[TABLE]

for $C=A/(1-\alpha)$ and $D=B/(1-\alpha)$, $\alpha$ is appearing in (13). Therefore, by $($20$)$ we have

[TABLE]

which implies

[TABLE]

4 Contraction speed estimates via coarse-fine methods

In this section we show an efficient way to estimate the rate of contraction of the discretized transfer operator and find the suitable $\overline{n}$ and $\alpha$ described in Section 3. Since $L_{\delta,\xi}$ is represented by a matrix, a first attempt to perform this task would be to iterate and estimate the norm of the iterate. This method is not very effective, since the matrix we should iterate is quite big. For this we implement a strategy in which we get information for the full matrix from the iterates of coarser versions of it. An earlier approach to this problem for deterministic systems can be found in [13].

In Lemmas 50 and 51 and Corollary 48 we have seen that

[TABLE]

We now prove the following lemma which bounds the distance between the powers of $L_{\xi}$ and $L_{\delta,\xi}$, provided that the noise has been applied at least once before the application of $L_{\xi}$ and $L_{\delta,\xi}$.

Lemma 28**.**

Let $\|L_{\delta,\xi}^{i}|_{V}\|_{L^{1}}\leq C_{i}$; let $\sigma$ be a linear operator such that $\sigma^{2}=\sigma$, $||\sigma||_{L^{1}}\leq 1$, and $\sigma\pi_{\delta}=\pi_{\delta}\sigma=\pi_{\delta}$; let $\Lambda=\sigma N_{\xi}\sigma L$.

Then $\forall n\geq 0$

[TABLE]

In particular the Lemma applies if:

$\sigma=Id$ and $\Lambda=L_{\xi}$; 2. 2.

$\sigma=\pi_{\delta^{\prime}}$ and $\Lambda=L_{\delta^{\prime},\xi}$ , for any $\delta^{\prime}$ such that $n\delta^{\prime}=\delta$ with $n\in\mathbb{R}$.

As a consequence, we obtain a way to bound the contraction rate of certain operator $L_{\delta^{\prime},\xi}$ on the zero-average space $V$ using the computed contraction rate for a coarser operator $L_{\delta,\xi}$. We remark that

[TABLE]

We remark that on the left-hand side we have an $n+1$ in (36), which guarantees that the noise has been applied at least once; this permits us to use Lemma 28 to estimate the second summand of the right-hand of (36).

Thus, if one is searching for an $n$ such that $||L_{\delta^{\prime},\xi}^{n+1}|_{V}\|_{L^{1}}<1$ this can be found and certified by using a suitable coarse version $L_{\delta,\xi}$, computing the norm of its iterates and using Lemma 28 in a way that the second hand of (36) is smaller than $1$.

Proof of Lemma 28.

Notice that as a consequence of the hypotheses we have

[TABLE]

because $\sigma g-\pi_{\delta}g=\sigma(1-\pi_{\delta})g$ applying Lemmas 51 and 52.

The proof is along the lines of what has been proved in Section 3.2. Indeed, we have

[TABLE]

and the thesis follows from the fact that $||L||_{L^{1}}\leq 1$, $||N_{\xi}||_{L^{1}}\leq 1$, $||\sigma||_{L^{1}}\leq 1$.

5 Estimating the average of an observable

As a result of the previous sections, we are able to obtain a precise approximation $\tilde{f}$ of $f_{\xi}$ in the $L^{1}$ norm. This is not enough in order to estimate the Lyapunov exponent which we recall can be defined as $\lambda_{\xi}=\int h~{}df_{\xi}$ where $h=\log|T^{\prime}|$. This is because $h$ is not in $L^{\infty}$ in the whole interval. In the Belouzov-Zhabotinsky case, there are two points where $h$ goes to infinity: the critical point and the point where the $|T^{\prime}|$ goes to $+\infty$. Outside of neighborhoods of these two points $h$ is bounded. Moreover, in these two neighborhoods $h$ still has bounded $L^{1}$ norm. This is enough to perform our estimates since we can have $L^{\infty}$ bounds on the stationary measure $f_{\xi}$; this allows to compute the Lyapunov exponent using alternately $L^{1}$ and $L^{\infty}$estimates on $f_{\xi}$ and $h$ in different sets. Therefore, we can join all these observations together to obtain a rigorous approximation of $\lambda_{\xi}$ applying the following strategy:

• •

we select a region $E$ of $[0,1]$ such that $h$ is in $L^{\infty}$ outside $E$;

• •

we estimate, on $E$, the quantities $\|f_{\xi}\|_{L^{\infty}(E)}$ and $\|h\|_{L^{1}(E)}$;

• •

we approximate (keeping rigorously track of the numerical errors) the integral $\int_{[0,1]}h~{}df_{\xi}$ with $\int_{[0,1]\setminus E}h~{}d\tilde{f}$ (discarding the set $E$ from the computation);

• •

we estimate the error in such an approximation in terms of $\|f_{\xi}-\tilde{f}\|_{L^{1}}$, $\|f_{\xi}\|_{L^{\infty}(E)}$, $\|h\|_{L^{1}(E)}$ and $\|h\|_{L^{\infty}([0,1]\setminus E)}$ (see Corollary 29).

In Subsection 5.1 we show how the error can be estimated in terms of the mentioned quantities. The remaining subsections are devoted to estimating $\|f_{\xi}\|_{L^{\infty}(E)}$ and $\|h\|_{L^{1}(E)}$, notice that in the case of uniform noise $\|f_{\xi}\|_{L^{\infty}(E)}\leq 1/\xi$, but we are able to obtain a better estimate via $\tilde{f}$.

5.1 Approximating the average using $L^{1}$ and $L^{\infty}$

estimates

In this section we assume that $\tilde{f}$ is an approximation of $f_{\xi}$ and both are probability measures, therefore $\tilde{f}-f_{\xi}$ has [math] average.

Corollary 29**.**

Let $f$ and $\tilde{f}$ be probability densities on the measure space $(X,m)$, both contained in $L^{1}$ and in $L^{\infty}$. Let $E\subset X$ be a Borel subset, and $H$ be an $L^{1}$ observable that is $L^{\infty}$ in $X\setminus E$. Then

[TABLE]

The proof of the corollary is straightforward, applying the following Lemma on the set $X\setminus E$

Lemma 30**.**

Let $(X,m)$ be a measure space, let $H\in L^{\infty}(X)$, and let $v\in L^{1}(X)$ a function having [math] average. Then we have

[TABLE]

Proof.

Indeed, for a constant $c$ we have

[TABLE]

because $v$ has [math] average, and this is clearly optimized taking $c=(\sup H+\inf H)/2$.

Remark 31**.**

Corollary 29 yields immediately an algorithm for estimating an observable that is $L^{1}$, and is $L^{\infty}$ outside a neighborhood of a finite number of points $s_{i}$ where it goes to $\infty$, as is the case with the observable $\log|T^{\prime}|$ of the system we are studying. In fact, there is a trade-off on the size of $E$, and we attempt different sets $E$ enclosing the $s_{i}$ with intervals of different sizes on order to obtain the tightest possible estimate on the error. Every such choice of the set $E$ yields an approximation of $\int_{X}Hf\ dm$ as $\int_{X\setminus E}H\tilde{f}\ dm$ and a bound for the error.

5.2 $L^{\infty}$ bounds for the stationary measure in an interval.

To estimate the average of the unbonded observable $h$ and apply Corollary 29, in this Subsection we obtain a bound for the $L^{\infty}$ norm of the invariant measure $f_{\xi}$ on intervals of a certain partition.  We derive it as a byproduct of the rigorous estimate of the $L^{1}$ error, and the algorithms explained in Sections 9.2.1 and 9.2.2, that allow to bound $\mathrm{Var}_{I}(N_{\xi}L\tilde{f})$ for each $I\in\Pi$.

Lemma 32**.**

Let $\Pi$ be a uniform partition, for each $I\in\Pi$ we have

[TABLE]

Proof.

Indeed,

[TABLE]

and the estimate follows because $\|N_{\xi}\|_{L^{1}\rightarrow L^{\infty}}\leq\|\rho_{\xi}\|_{L^{\infty}}$.

5.3 $L^{1}$ bounds on $\log|T^{\prime}|$

In this section we compute explicit bounds for the $L^{1}$-norm of $\log|T^{\prime}|$ for the map defined in Section 2 in a neighborhood of the points where it is not bounded, as required to apply Remark 31. As can be deduced from its definition in Section 2, we need to do so in intervals enclosing $x=0.125$ and $x=0.3$. We omit the proofs, as they are all very elementary.

Lemma 33**.**

For $0.125-2^{-6}<u<0.125<v<0.125+2^{-6}$ we have

[TABLE]

Lemma 34**.**

For $0.2<x<0.3$ we have

[TABLE]

for $d_{1}=\frac{1}{3}(\frac{2}{3}\cdot(0.175)^{-5/3}+(0.175)^{-2/3})$ and $d_{2}=[a+(x-0.125)^{1/3}-\frac{1}{3}|x-0.125|^{-2/3}]/(0.3-x)$, where $a=0.50607356\dots$ as in Section 2.

Lemma 35**.**

For $0.3<x<0.303$ we have

[TABLE]

6 Computation details and results

We give here some details about the code performing our computer aided estimates and about the results. The main algorithm is written in Python using the Sage framework  and interval arithmetics ([29]), some critical parts are written in C++ and uses (optionally) the GPU, this in particular is used for the iterartion of large matrices needed to apply the methods of Section 4. Such parts have been optimized to use high performance computing; even so each contraction test has required a time of the of order of a week. 333The contraction time estimates were run on a Asus GeForce GTX 1050Ti, 4GB of Ram GPU installed in a desktop computer with an AMD A4-6300 3.4 Ghz processor and 8 Gb of Ram. The matrices were assembled on a Dell R710 server with 2 sixcore Xeon 5660 2.8 Ghz processors and 24 Gb of Ram. The code we used can be found at

http://im.ufrj.br/~maurizio.monge/wordpress/rigorous_computation_dyn/

Table 6 contains the result of the computer aided estimates we performed and the values of the parameters used in these estimates. In Figure 4 we summarize with a graph the most important information contained in the table. The graph shows intervals enclosing the Lyapunov exponent at the selected noise values. These are the final result of our computer aided estimates; it is worth to remark again that the estimated requiring more computational power are the ones performed to prove that the Lyapunov exponent is positive for small size of the noise.

To ease the understanding of Table 6, we now outline some more details of the implementation of our algorithm and how the parameters take place in the algorithm’s execution; the columns are ordered as they are subsequently used in the algorithm, or deduced from previous quantities and via computations.

One of the main ingredients and goals of the computer aided estimates is to compute the number of iterates of the transfer operator needed to contract the zero average space (see Item 1 of Section 3.1). For this we apply the ”coarse fine” methods explained in Section 4. For each value of the noise amplitude, denoted by $\xi$, we build a coarse discretization $L_{\xi,\delta_{\text{contr}}}$ of the operator $L_{\xi}$, on a partition of coarse size $\delta_{\text{contr}}$. We then compute values of $n_{\text{contr}}$ and $\alpha_{\text{contr}}$ represented in Table 6, which satisfy

[TABLE]

and compute explicit bounds for $||L_{\xi,\delta_{\text{contr}}}^{i}|_{V}||_{L^{1}}$ (whose value does not appear in the table). The algorithm used for these finite dimensional estimates is the same as in [12] and there explained.

We consider then a finer partition size $\delta$ and use the coarse fine estimates of Section 4 to compute $\alpha$ and $\sum_{i=0}^{n_{\text{contr}}}C_{i}$, where $C_{i}=||L_{\xi,\delta}^{i}|_{V}||_{L^{1}}$ and

[TABLE]

The same bounds work for $L_{\xi}$, i.e.

[TABLE]

Notice that the bigger the ${n_{\text{contr}}}$, the worse is going to be the estimate on the $L^{1}$ norm on $(L_{\xi}^{n}-L_{\xi,\delta_{\text{contr}}}^{n})N_{\xi}$, i.e., the error coming from the coarse-fine inequality. While increasing $n_{\text{contr}}$ permits us to find smaller $\alpha_{\text{contr}}$, this may not imply that the corresponding $\alpha$ is smaller. Our algorithms attempts to find the best compromise. Needless to say, in practice this procedure may fail, if unable to detect any contraction in a reasonable time. This might happen for example if the original system is not mixing.

To estimate an upper bound to the $L^{1}$ error in the computation of the stationary measure we use the results shown in Section 3. The column ”A priori…” contains the a priori estimate on the $L^{1}$ error of the approximation of the measure on the partition of size $\delta$ as given using the results in Subsection 3.2.1 while the column ”Refined…” contains the $L^{1}$ error when we use the bootstrapping tecniques of Subsection 3.3.

Once we have a good approximation of the invariant measure, we compute an approximation of the Lyapunov exponent by computing an integral on a partition of size $\delta_{\text{est}}$, using Section 5; the computed intervals enclosing rigorously the Lyapunov exponent are contained in the last column. This allows to prove Theorem 2.