Interval of effective time-step size for the numerical computation of nonlinear ordinary differential equations

Abstract The computational uncertainty principle states that the numerical computation of nonlinear ordinary differential equations (ODEs) should use appropriately sized time steps to obtain reliable solutions. However, the interval of effective step size (IES) has not been thoroughly explored theoretically. In this paper, by using a general estimation for the total error of the numerical solutions of ODEs, a method is proposed for determining an approximate IES by translating the functions for truncation and rounding errors. It also illustrates this process with an example. Moreover, the relationship between the IES and its approximation is found, and the relative error of the approximation with respect to the IES is given. In addition, variation in the IES with increasing integration time is studied, which can provide an explanation for the observed numerical results. The findings contribute to computational step-size choice for reliable numerical solutions.


Introduction
Many works have shown the time-step sensibility of nonlinear dynamical systems. Chou (2000, 2001) and Li (2000) proposed the computational uncertainty principle (CUP) for nonlinear systems of ordinary differential equations (ODEs) under a finite machine precision. The CUP states that using different time-step sizes usually results in different effective computation times (ECTs) and that the maximal ECT (MECT), achieved using the optimal step size (OS), gives the best result. Wang and Huang (2006) focused on Lorenz systems, and reported that the maximum prediction time sensitively relies on the timestep size under certain conditions. Teixeira, Reynolds, and Judd (2007) found the time-step size to affect not only Lorenz systems but also a quasi-geostrophic model. Liu et al. (2015) studied the Global/Regional Assimilation and Prediction System mesoscale numerical forecast, and gave a preliminary explanation of the applicability of OS theory to complicated partial differential equations (PDEs).
The CUP presented by Chou (2000, 2001) theoretically explained the time-step sensibility of nonlinear ODEs, which has been cited by many other researches (Hu and Chou 2004;Li and Wang 2008;Liu et al. 2015;Wang, Li, and Li 2012;Wang, Liu, and Li 2014). In particular, based on the CUP, Wang, Li, and Li (2012) deduced a general ECT function of step size, which explained the experimental formulae proposed by Teixeira, Reynolds, and Judd (2007).
Through a large number of numerical experiments, Li, Zeng, and Chou (2000) introduced the concept of the interval of effective step size (IES) of ODEs. Presenting the IES profiles obtained from numerical results (Figure 1), Li, Zeng, and Chou (2000) suggested that numerical solutions are reliable when step sizes belong to the IES. In such cases, if we know the theoretical formulae of lower and upper bounds of the IES corresponding to a certain error tolerance, it will guide the choices of effective step sizes in computations. However, there has been little relevant prior research in this regard. This paper explores the IES for nonlinear ODEs based on the studies of Chou (2000, 2001).
2 ) denote the IES at integral time t under a 2. Method for determining U * t , an approximation of U t First, defining (h cross , E cross ) as the intersection of the functions E 1 (h) and E 2 (h), one gets where C 12 = C 1 /C 2 . Besides, Chou (2000, 2001) stated that Ẽ (h) reaches its minimum E min when the step size h takes the value of OS, and when the OS denoted by H, there are Then, we simultaneously translate the functions E 1 (h) and E 2 (h) so as to move the coordinates of their intersection from (h cross , E cross ) to the lowest point (H, E min and E * 2 (h) equal ̃t respectively to obtain two new equations whose solutions are Taking the situation of p = C 1 = C 2 = 1 as an example, the above process is shown in Figure 2. 3. Relationship between U t and U * t From the above definitions we find: as step size h decreases, Ẽ (h) initially monotonically decreases to its lowest point (H, E min ) before monotonically increasing; E * 1 (h) is a monotonically decreasing function, whereas E * 2 (h) is a monotonically increasing function of h, and their intersection is (H, E min ); it is easy to prove that when h < H, E * 1 (h)>Ẽ(h) is always true, and when h > H, E * 2 (h)>Ẽ(h) is true. Given these, we have: From the statements above we know that U * t ⊂ U t when ̃t > E min , and U * t = U t = {H} when ̃t = E min ; however, when ̃t < E min , both U t and U * t are empty sets. These results indicate that U * t ⊆ U t is always true, which suggests that U * t is suitable for serving as an approximate interval U t . In , and E min = C 2 (2p + 1) ; wheñt < E min , h t,1 and h t,2 do not exist, and h * t,1 > h * t,2 , which does not conform to the definition of U * t .
given error tolerance δ. To obtain U t , it is necessary to give a general formula of the numerical error E(t, h) for the solutions of nonlinear ODEs. In numerical calculation, E(t, h) is usually composed of three parts: truncation error, which is caused by differential equation discretization (Gear 1971;Stoer and Bulirsch 1993); round-off error, which is due to limitations of computer precision Chou 2000, 2001); and initial error (Ding and Li 2008a, 2008b. From Li, Zeng, and Chou (2001, Equations (60) and (83)), it can be shown that where E 1 (h) = C 1 h −0.5 is relevant to the round-off error; E 2 (h) = C 2 h p is relevant to the truncation error, and p is the order of the numerical method; e (0) is relevant to the initial error, and The way to estimate C 1 and C 2 , and details of other parameters, are given in Li, Zeng, and Chou (2001, Equations (60) and (83)). Letting ̃t = /C(t)-Ne (0) and Ẽ (h) = E 1 (h)+E 2 (h), Equation (1) indicates that h t,1 and h t,2 should be the solutions of the equation For a fixed value of t, Equation (3) is a nonlinear equation associated with h, which can be solved numerically to obtain approximate values of h t,1 and h t,2 by methods such as fixed-point iteration and Newtonian iteration (Suli and Mayers 2003); however, it is usually hard to provide function expressions for h t,1 and h t,2 with these methods. This article aims to derive explicit formulae for h t,1 and h t,2 , so as to give a general approximate explicit expression for U t . (1) (3) E(h) =̃t.  Figure 1. ies profiles obtained using the optimal searching method, when computing the solutions of the x-component of the Lorenz equation using the fourth-order runge-Kutta method for the initial value (5, 5, 10) and r = 28 and for 121 different step sizes in the range 10 −7 -10 −1 . source: Li, Zeng, and chou (2000, plate i-2(c)).
notes: here, the step size h is plotted as a logarithm (to base 10) and time is non-dimensional. the grey solid line is for machine single precision and the black dotted line is for double precision.
addition, to obtain a non-empty set U * t , we suppose that ̃t ≥ E min in the following discussion.
Next, we estimate the error of the approximation U * t with respect to U t . For this purpose, let Δ t,1 = |h * t,1 −h t,1 | and Δ t,2 = |h * t,2 −h t,2 |. Assuming that ̃t ≥ E min , the relative errors of h    Figure 3. schematic representation of the variations in the ies U t (solid line) and its approximate interval U * t (dotted line) with increasing integration time t.