Article Contents
Article Contents

# Rigorous numerics for ODEs using Chebyshev series and domain decomposition

• * Corresponding author: Jan Bouwe van den Berg

The first author was partially supported by NWO-VICI grant 639033109

• In this paper we present a rigorous numerical method for validating analytic solutions of nonlinear ODEs by using Chebyshev-series and domain decomposition. The idea is to define a Newton-like operator, whose fixed points correspond to solutions of the ODE, on the space of geometrically decaying Chebyshev coefficients, and to use the so-called radii-polynomial approach to prove that the operator has an isolated fixed point in a small neighborhood of a numerical approximation. The novelty of the proposed method is the use of Chebyshev series in combination with domain decomposition. In particular, a heuristic procedure based on the theory of Chebyshev approximations for analytic functions is presented to construct efficient grids for validating solutions of boundary value problems. The effectiveness of the proposed method is demonstrated by validating long periodic and connecting orbits in the Lorenz system for which validation without domain decomposition is not feasible.

Mathematics Subject Classification: 37M99, 65G20, 34B15.

 Citation:

• Figure 1.  Approximations of the solution of $u' = u(1-u)$ with $u(0) = \frac{1}{2}$. In blue (partly obscured by red) based on a truncated Chebyshev series for $t\in[0,50]$ with $N = 75$ terms and geometric decay factor $\nu = 1.05$ (see Section 2) with error at most $6.3 \cdot 10^{-10}$ (uniformly). In red based on a truncated Taylor series for $t\in [0,3]$ with $N = 225$ terms and $\nu = 1$ with uniform error control $2.7 \cdot 10^{-4}$. The computation times (including the proof of the error bound) for both cases are comparable. The rigorous error bounds are based on the truncated series only, not on the availability of an explicit expression for the solution

Figure 2.  The connecting orbit from the positive eye to the origin in the Lorenz equation with classical parameters. The time of flight between the local (un)stable manifolds is $L = 30$. The geometric objects colored in red and green are representations of $W_{ \rm{loc}}^{s} \left( q^{+} \right)$ and $W_{ \rm{loc}}^{u} \left( {\bf{0}} \right)$, respectively

Figure 3.  (a) A periodic orbit on the Lorenz attractor of period $L \approx 25.03$ validated with $m = 35$ subdomains. (b) A semi-logarithmic plot of the coefficients in the Chebyshev series on all subdomains for the three components of the solution

Figure 4.  (a) A typical periodic orbit near the homoclinic connection. The period of the orbit is $L \approx 100.25$. Notice the sharp turn of the orbit near the origin. (b) The $x$, $y$ and $z$ components of the orbit. The three components are fairly flat for a relatively long time. These flat parts correspond to the part of the orbit near the equilibrium where the dynamics are slow

Figure 5.  The fifteen components of the validated periodic orbit in the coupled Lorenz system (4). The period of the periodic orbit is $L \approx 1.70$. We have grouped the components, which in the synchronized setting would correspond to the "same component", together. In each case the components with higher indices are colored with a lighter shade of blue

Figure 6.  (a) A validated solution of the initial value problem with initial condition $p_{0} = \begin{bmatrix} 6 & 10 & 8 \end{bmatrix}^{T}$ and integration time $L = 11.2$. (b) The three components of the validated orbit. The dashed lines in grey indicate the six chosen integration times $L \in \{2.3, 4.8, 6.1, 7.6, 9.5, 11.2\}$. At the $k$-th integration time, ordered from small to large, the orbit transitioned $k$ times between the eyes

Figure 7.  The approximate complex singularities of the validated periodic orbit. Note that the time variable has been rescaled to $[0,1]$. The complex singularities were computed with the procedure described in Section 4.2

Figure 8.  (a) A plot of the two grids $\left(i, t_{i} \right)_{i = 0}^{33}$ and $\left(i, \tau_{i} \right)_{i = 0}^{34}$ corresponding to $m = 33$ and $m = 34$, respectively. (b) A plot of the grid-points $t_{32}$, $\tau_{33}$ and the complex singularities (colored in red) which determine the size of the Bernstein ellipses associated to $\left[ t_{32},1 \right]$ and $\left[ \tau_{33},1 \right]$

Figure 9.  The dependence of the period $L$ as a function of $\rho$ obtained via non-rigorous pseudo-arclength continuation

Figure 10.  (a) The complex singularities of the connecting orbit. (b) The grid, determined by the algorithm in Section 4, on which the connecting orbit was validated

Figure 11.  (a) The complex singularities of the validated periodic orbit in the coupled Lorenz system for $K = 5$. (b) A semi-logarithmic plot of the coefficients in the Chebyshev series of the periodic orbit for $K = 5$ for all subdomains and components. The decay rates are approximately the same for each subdomain and component. The black line corresponds to the theoretically predicted decay rate determined by the size of the smallest Bernstein-ellipse. The theoretically predicted decay rate coincides (approximately) with the observed decay rate, which indicates that the approximation of the complex singularities (near the real-axis) was accurate

Figure 12.  The thirty components of the validated periodic orbit in the coupled Lorenz system (4). The period of the periodic orbit is $L \approx 1.68$. We have grouped the components, which in the synchronized setting would correspond to the "same component", together. In each case the components with higher indices are colored with a lighter shade of blue

Figure 13.  (a) Step sizes used by awa for $L = 11.2$. (b) Step sizes used by verifyode for $L = 6.1$. Approximately $45 \%$, $40 \%$ and $10 \%$ of the time steps are located around small neighborhoods of $t = 4$, $t = 5$ and (just after) $t = 6$, respectively. Moreover, these time steps are of minimal size

Table 1.  Numerical results for a validated periodic orbit of period $L \approx 25.0271$ in the classical Lorenz system. In each case the number of modes $N_{i}$ per domain was approximately the same. The number $\bar N$ denotes the (rounded) average number of modes per domain

 $m$ $\bar N$ $\dim \mathcal{X}^{N}_{\nu}$ $\nu$ $\nu_{e}$ $I_{m,\nu}$ $32$ $79$ $7618$ $1.1148$ $1.6413$ $\left[ 4.5581 \cdot 10^{-10}, \ 1.0217 \cdot 10^{-7} \right]$ $33$ $76$ $7498$ $1.1278$ $1.6837$ $\left[ 3.1192 \cdot 10^{-10}, \ 1.5371 \cdot 10^{-7} \right]$ $34$ $63$ $6433$ $1.1432$ $1.8749$ $\left[ 3.8158 \cdot 10^{-10}, \ 1.1950 \cdot 10^{-7} \right]$ $35$ $61$ $6457$ $1.1529$ $1.9056$ $\left[ 3.3529 \cdot 10^{-10}, \ 1.0320 \cdot 10^{-7} \right]$ $36$ $61$ $6610$ $1.1564$ $1.9115$ $\left[ 3.1282 \cdot 10^{-10}, \ 1.1097 \cdot 10^{-7} \right]$

Table 2.  Numerical results for two periodic orbits near the homoclinic connection. The periodic orbit of period $L \approx 100.2554$ was in both cases validated on a grid for which (the same) six subdomains were used to approximate the non-flat part of the orbit

 $L$ $m$ $\dim \mathcal{X}^{N}_{\nu}$ $I_{m,\nu}$ $4.5473$ $1$ $566$ $\left[ 4.4568 \cdot 10^{-11}, \ 8.4403 \cdot 10^{-6} \right]$ $100.2554$ $7$ $5894$ $\left[ 1.5186 \cdot 10^{-11}, \ 7.4915 \cdot 10^{-8} \right]$ $100.2554$ $506$ $8441$ $\left[ 1.5174 \cdot 10^{-11}, \ 4.2914 \cdot 10^{-6} \right]$

Table 3.  Numerical results for the connecting orbit from $q^{+}$ to the origin. The interval $I_{m,\nu}$ is the set of admissible radii on which the radii-polynomials were proven to be strictly negative

 $m$ $\bar N$ $\dim \Pi_{N} \left( \mathcal{X}_{\nu} \right)$ $\nu$ $I_{m,\nu}$ $55$ $42$ $6864$ $1.3710$ $\left[ 6.4412 \cdot 10^{-9}, \ r^{\ast} \right]$

Table 4.  Numerical results for a validated periodic orbit in systems of five and 10 coupled Lorenz vector fields. We took $m = 4$ for both proofs. The interval $I_{m,\nu}$ is the set of admissible radii on which the radii-polynomials were proven to be strictly negative

 $3K$ $L$ $N$ $\dim \Pi_{N} \left( \mathcal{X}_{\nu} \right)$ $\nu$ $I_{m,\nu}$ $15$ $1.7046$ $[57 \; \; \; 54 \; \; \; 59 \; \; \; 55]$ $3376$ $1.1807$ $\left[ 4.39 \cdot 10^{-10}, \ 1.66 \cdot 10^{-6} \right]$ $30$ $1.6839$ $[54 \; \; \; 53 \; \; \; 56 \; \; \; 54 ]$ $6511$ $1.1849$ $\left[ 3.33 \cdot 10^{-10}, \ 2.21 \cdot 10^{-6} \right]$

Table 5.  Numerical results for the validated integration of an IVP in the Lorenz system with initial condition $p_{0} = [10 \; \; 6 \; \; 8]^{T}$ and different integration times $L$. The orbits were validated using the verified integrators awa and verifyode in INTLAB using the default settings. In particular, $10$-th order Taylor approximations were used and the minimal step size was $10^{-4}$. The integer $n_{q^{\pm}}$ denotes the number of transitions between the eyes. The third and fourth column correspond to the number of time steps used for integration. The integrator verifyode failed, due to too large enclosures, for $n_{q^{\pm}} \geq 4$

 Steps Steps $C^{0}$-error $C^{0}$-error $L$ $n_{q^{\pm}}$ awa verifyode awa verifyode $2.3$ $1$ $719$ $150$ $5.3774 \cdot 10^{-11}$ $1.5439 \cdot 10^{-3}$ $4.8$ $2$ $1471$ $2643$ $9.9833 \cdot 10^{-10}$ $9.8674 \cdot 10^{-2}$ $6.1$ $3$ $1873$ $5213$ $1.186 \cdot 10^{-8}$ $4.0664$ $7.6$ $4$ $2342$ $-$ $1.186 \cdot 10^{-8}$ $-$ $9.5$ $5$ $2901$ $-$ $7.1604 \cdot 10^{-8}$ $-$ $11.2$ $6$ $3411$ $-$ $5.1561 \cdot 10^{-7}$ $-$

Table 6.  Numerical results for the validated integration of an IVP in the Lorenz system with initial condition $p_{0} = [10 \; \; 6 \; \; 8]^{T}$ and different integration times $L$ using the proposed method. The integer $n_{q^{\pm}}$ denotes the number of transitions between the eyes. The number $\bar N$ denotes the (rounded) average number of modes per subdomain. The interval $I_{m,\nu}$ is the set of admissible radii on which the radii-polynomials were proven to be strictly negative and $\nu$ is the associated weight with which the validation was performed. Note that the left-endpoints of $I_{m,\nu}$ are rigorous upper bounds for the $C^{0}$-error

 $L$ $n_{q^{\pm}}$ $m$ $\bar{N}$ $\dim \Pi_{N} \left( \mathcal{X}_{\nu} \right)$ $\nu$ $I_{m,\nu}$ $2.3$ $1$ $2$ $118$ $705$ $1.057$ $\left[2.0218 \cdot 10^{-10}, 9.0595 \cdot 10^{-4}\right]$ $4.8$ $2$ $4$ $122$ $1458$ $1.087$ $\left[6.6631 \cdot 10^{-9}, 2.7927 \cdot 10^{-5}\right]$ $6.1$ $3$ $5$ $125$ $1881$ $1.1$ $\left[2.1101 \cdot 10^{-8}, 9.4164 \cdot 10^{-6} \right]$ $7.6$ $4$ $8$ $122$ $2937$ $1.094$ $\left[2.1476 \cdot 10^{-8}, 5.1128 \cdot 10^{-6} \right]$ $9.5$ $5$ $11$ $124$ $4101$ $1.113$ $\left[7.7316 \cdot 10^{-8}, 1.1828 \cdot 10^{-6} \right]$ $11.2$ $6$ $15$ $123$ $5547$ $1.136$ $\left[1.6234 \cdot 10^{-7}, 5.1904 \cdot 10^{-7} \right]$
•  [1] G. Arioli and H. Koch, Computer-assisted methods for the study of stationary solutions in dissipative systems, applied to the Kuramoto-Sivashinski equation, Arch. Ration. Mech. Anal., 197 (2010), 1033-1051.  doi: 10.1007/s00205-010-0309-7. [2] G. Arioli and H. Koch, Existence and stability of traveling pulse solutions of the FitzHugh-Nagumo equation, Nonlinear Anal., 113 (2015), 51-70.  doi: 10.1016/j.na.2014.09.023. [3] G. Arioli and H. Koch, Integration of dissipative partial differential equations: A case study, SIAM J. Appl. Dyn. Syst., 9 (2010), 1119-1133.  doi: 10.1137/10078298X. [4] A. W. Baker, M. Dellnitz and O. Junge, A topological method for rigorously computing periodic orbits using Fourier modes, Discrete Contin. Dyn. Syst., 13 (2005), 901-920.  doi: 10.3934/dcds.2005.13.901. [5] R. Barrio, A. Dena and W. Tucker, A database of rigorous and high-precision periodic orbits of the Lorenz model, Comput. Phys. Comm., 194 (2015), 76-83.  doi: 10.1016/j.cpc.2015.04.007. [6] M. Berz and K. Makino, Rigorous reachability analysis and domain decomposition of Taylor models, in International Workshop on Numerical Software Verification, Lecture Notes in Computer Science, 10381, Springer, Cham, 2017, 90-97. doi: 10.1007/978-3-319-63501-9_7. [7] M. Berz and K. Makino, Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models, Reliab. Comput., 4 (1998), 361-369.  doi: 10.1023/A:1024467732637. [8] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Dover Publications, Inc., Mineola, NY, 2001. [9] M. Breden and J.-P. Lessard, Polynomial interpolation and a priori bootstrap for computer-assisted proofs in nonlinear ODEs, Discrete Contin. Dyn. Syst. Ser. B, 23 (2018), 2825-2858.  doi: 10.3934/dcdsb.2018164. [10] M. Breden, J.-P. Lessard and M. Vanicat, Global bifurcation diagrams of steady states of systems of PDEs via rigorous numerics: A 3-component reaction-diffusion system, Acta Appl. Math., 128 (2013), 113-152.  doi: 10.1007/s10440-013-9823-6. [11] B. Breuer, J. Horák, P. J. McKenna and M. Plum, A computer-assisted existence and multiplicity proof for travelling waves in a nonlinearly supported beam, J. Differential Equations, 224 (2006), 60-97.  doi: 10.1016/j.jde.2005.07.016. [12] F. Bünger, A Taylor model toolbox for solving ODEs implemented in MATLAB/INTLAB, J. Comput. Appl. Math., 368 (2020), 20pp. doi: 10.1016/j.cam.2019.112511. [13] J. Burgos-García, J.-P. Lessard and J. D. Mireles James, Spatial periodic orbits in the equilateral circular restricted four-body problem: Computer-assisted proofs of existence, Celestial Mech. Dynam. Astronom., 131 (2019), 36pp. doi: 10.1007/s10569-018-9879-8. [14] X. Cabré, E. Fontich and R. de la Llave, The parameterization method for invariant manifolds. â…¢. Overview and applications, J. Differential Equations, 218 (2005), 444-515.  doi: 10.1016/j.jde.2004.12.003. [15] CAPD, Computer assisted proofs in dynamics., Available from: http://capd.ii.uj.edu.pl. [16] B. A. Coomes, H. Koçak and K. J. Palmer, Transversal connecting orbits from shadowing, Numer. Math., 106 (2007), 427-469.  doi: 10.1007/s00211-007-0065-2. [17] COSY, Center for Beam Theory and Dynamical Systems., Available from: http://bt.pa.msu.edu/index.htm. [18] S. Day, Y. Hiraoka, K. Mischaikow and T. Ogawa, Rigorous numerics for global dynamics: A study of the Swift-Hohenberg equation, SIAM J. Appl. Dyn. Syst., 4 (2005), 1-31.  doi: 10.1137/040604479. [19] S. Day, J.-P. Lessard and K. Mischaikow, Validated continuation for equilibria of PDEs, SIAM J. Numer. Anal, 45 (2007), 1398-1424.  doi: 10.1137/050645968. [20] T. A. Driscoll, N. Hale and L. N. Trefethen, Chebfun Guide, Pafnuty Publications, Oxford, 2014. Available from: http://www.chebfun.org/. [21] O. Fogelklou, W. Tucker, G. Kreiss and M. Siklosi, A computer-assisted proof of the existence of solutions to a boundary value problem with an integral boundary condition, Commun. Nonlinear Sci. Numer. Simul., 16 (2011), 1227-1243.  doi: 10.1016/j.cnsns.2010.07.008. [22] M. Gameiro, T. Gedeon, W. Kalies, H. Kokubu, K. Mischaikow and H. Oka, Topological horseshoes of traveling waves for a fast-slow predator-prey system, J. Dynam. Differential Equations, 19 (2007), 623-654.  doi: 10.1007/s10884-006-9013-6. [23] M. Gameiro and J.-P. Lessard, Analytic estimates and rigorous continuation for equilibria of higher-dimensional PDEs, J. Differential Equations, 249 (2010), 2237-2268.  doi: 10.1016/j.jde.2010.07.002. [24] P. Gonnet, R. Pachón and L. N. Trefethen, Robust rational interpolation and least-squares, Electron. Trans. Numer. Anal., 38 (2011), 146-167. [25] A. Hungria, J.-P. Lessard and J. D. Mireles James, Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach, Math. Comp., 85 (2016), 1427-1459.  doi: 10.1090/mcom/3046. [26] M. Kashiwagi, KV - A C++ library for verified numerical computation., Available from: http://verifiedby.me/kv/index-e.html. [27] J.-P. Lessard, J. D. Mireles James and J. Ransford, Automatic differentiation for Fourier series and the radii polynomial approach, Phys. D, 334 (2016), 174-186.  doi: 10.1016/j.physd.2016.02.007. [28] J.-P. Lessard, J. D. Mireles James and C. Reinhardt, Computer assisted proof of transverse saddle-to-saddle connecting orbits for first order vector fields, J. Dynam. Differential Equations, 26 (2014), 267-313.  doi: 10.1007/s10884-014-9367-0. [29] J.-P. Lessard and C. Reinhardt, Rigorous numerics for nonlinear differential equations using Chebyshev series, SIAM J. Numer. Anal, 52 (2014), 1-22.  doi: 10.1137/13090883X. [30] R. J. Lohner, Computation of guaranteed enclosures for the solutions of ordinary initial and boundary value problems, in Computational Ordinary Differential Equations, Inst. Math. Appl. Conf. Ser. New Ser., 39, Oxford Univ. Press, New York, 1992,425-435. [31] R. J. Lohner, Enclosing the solutions of ordinary initial and boundary value problems, in Computerarithmetic, Teubner, Stuttgart, 1987,255-286. [32] K. Makino and M. Berz, Verified computations using Taylor models and their applications, in International Workshop on Numerical Software Verification, Lecture Notes in Computer Science, 10381, Springer, Cham, 2017, 3-13. doi: 10.1007/978-3-319-63501-9_1. [33] J. C. Mason and D. C. Hanscomb, Chebyshev Polynomials, Chapman & Hall/CRC, Boca Raton, FL, 2003. [34] R. Pachón, P. Gonnet and J. van Deun, Fast and stable rational interpolation in roots of unity and Chebyshev points, SIAM J. Numer. Anal, 50 (2012), 1713-1734.  doi: 10.1137/100797291. [35] L. N. Rump, INTLAB - INTerval LABoratory, in Developments in Reliable Computing, Kluwer Academic Publishers, Dordrecht, 1999, 77-104. Available from: http://www.ti3.tuhh.de/rump/. [36] L. N. Trefethen, Approximation Theory and Approximation Practice, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2013. [37] J. B. van den Berg, C. M. Groothedde and J. F. Williams, Rigorous computation of a radially symmetric localized solution in a Ginzburg-Landau problem, SIAM J. Appl. Dyn. Syst, 14 (2015), 423-447.  doi: 10.1137/140987973. [38] J. B. van den Berg, J.-P. Lessard and K. Mischaikow, Global smooth solution curves using rigorous branch following, Math. Comp., 79 (2010), 1565-1584.  doi: 10.1090/S0025-5718-10-02325-2. [39] J. B. van den Berg, J. D. Mireles James, J.-P. Lessard and K. Mischaikow, Rigorous numerics for symmetric connecting orbits: Even homoclinics of the Gray-Scott equation, SIAM J. Math. Anal., 43 (2011), 1557-1594.  doi: 10.1137/100812008. [40] J. B. van den Berg and R. Sheombarsing, MATLABcode for rigorous numerics for ODEs using Chebyshev-series and domain decomposition, 2021., Available from: http://www.few.vu.nl/janbouwe/code/domaindecomposition. [41] J. B. van den Berg and R. Sheombarsing, Validated computations for connecting orbits in polynomial vector fields, Indag. Math. (N.S.), 31 (2020), 310-373.  doi: 10.1016/j.indag.2020.01.007. [42] M. Webb, Computing complex singularities of differential equations with Chebfun, SIAM Undergraduate Research Online, 2011. Available from: http://evoq-eval.siam.org/Portals/0/Publications/SIURO/Vol6/COMPUTING_COMPLEX_SINGULARITIES_Differential.pdf?ver=2018-04-06-151849-873. [43] D. Wilczak and P. Zgliczynski, $C^n$-Lohner algorithm, Scheade Informaticae, 20 (2011), 9-46. [44] D. Wilczak and P. Zgliczynski, Heteroclinic connections between periodic orbits in planar restricted circular three-body problem - A computer assisted proof, Comm. Math. Phys., 234 (2003), 37-75.  doi: 10.1007/s00220-002-0709-0. [45] N. Yamamoto, A numerical verification method for solutions of of boundary value problems with local uniqueness by Banach's fixed-point theorem, SIAM J. Numer. Anal, 35 (1998), 2004-2013.  doi: 10.1137/S0036142996304498. [46] P. Zgliczynski, $C^1$-Lohner algorithm, Found. Comput. Math., 2 (2002), 429-465.  doi: 10.1007/s102080010025.

Figures(13)

Tables(6)