Numerical Validations of the Tangent Linear Model for the Lorenz Equations

Abstract: The validity of the tangent linear model (TLM) is studied numerically using the example of the Lorenz equations in this paper. The relationship between the limit of the validity time of the TLM and initial perturbations for the Lorenz equations is investigated using the Monte Carlo sampling method. A new error function between the nonlinear and the linear evolution of the perturbations is proposed. Furthermore, numerical sensitivity analysis is carried to establish the relationship between parameters and the validity of the TLM, such as the initial perturbation, the prediction time, the time step size and so on, by the method of mathematical statistics.


Introduction
The study of chaotic systems has attracted a lot of attentions in the literature. The chaotic system in a population model is studied in depth by the Lorenz et al. [Lorenz (1963[Lorenz ( , 1965[Lorenz ( , 1982; Yassen (2005); Wu, Xie, Fang et al. (2007); Zhao, Xing and Yu (2009)]. The chaotic system is characterized by a "sensitive dependence on initial condition" [Lorenz (1963)], such that the predictability of the future state is often severely limited by the chaotic dynamics of the system. The Lorenz model is one of the most popular models in dynamic systems, and has been studied by Wu et al. [Wu, Xie, Fang et al. (2007); Richter (2001) ;Yang, Chen and Yau (2002); Aniszewska and Rybaczuk (2005); Foias and Jolly (2005); Ding and Li (2012)]. In Shepelev et al. [Shepelev, Strelkova and Anishchenko (2018)], the transition from the regime of spatio-stationary structures and solitary states to complete incoherence in the network of nonlocally coupled Lorenz systems is numerically studied and a typical property of nonhyperbolic chaotic systems is obtained. Faghihnaini et al. [Faghihnaini and Shen (2018)] shows that the feedback loop produces periodic modes for three-and five-dimensional nondissipative Lorenz model. In Bougoffa et al. [Bougoffa, Al-Awfi and Bougouffa (2018)], the passage for Lorenz system to the third order non linear differential equation is established and for some special parameter values the obtained 1 equation can be reduced to the well known equations. The Lorenz model is the simplest model for dynamics of the convective layers and dynamics of closed convection loops. One application of the Lorenz model is the motion of a cell of fluid that is cooled from above and warmed from below. If the temperature difference between bottom and top is small, then the fluid of bottom will slowly rises because of the heating. If the temperature difference is large, then then he fluid from bottom may rise, and the colder fluid from button may drop at the same time [Lorenz (1963)]. The Lorenz model specifies a three-dimensional flow by the following non-dimensional equations: where the parameters σ, r and b represent the Prandtl number, the Rayleigh number, and a geometric factor, respectively. The parameters have the values σ = 10, r = 28 and b = 8/3, for which the well-known "butterfly" attractor exists. In Zhang et al. [Zhang, Mu, Zhou et al. (2017); Zhang, Liao and Zhang (2015)], the dynamical behavior of a generalized Lorenz system is derived based on the stability theory of dynamical systems so that the choice of parameters is reasonable. In addition, the state variables x, y, z represent measures of fluid velocity and the spatial temperature distribution in the fluid layer under gravity. The same sign x and y in the ODEs above indicates that the warm fluid is rising and the cold fluid is descending, which shows that a convection is taking place. The linear approach assumes that the initial perturbation is sufficiently small such that its evolution can be governed approximately by the TLM instead of the nonlinear model. In the Lorenz model, the validity of the TLM depends on the initial perturbation, the time prediction, and the parameters in the model. However, there have been few papers addressing the numerical relationship between the parameters and the validity of the TLM. In the present study, the TLM is effective as long as satisfies a limited condition and we provide a confidence interval so that the validity can be described more precisely. Meanwhile we show the numerical relationship between parameters and the validity of the TLM, by the method of mathematical statistics. Currently, there have been several predictability studies on the Lorenz model. The linear approach assumes that the initial perturbation is sufficiently small such that its evolution can be governed approximately by the TLM of the nonlinear model, and the computation of the linear fastest growing perturbation is reduced to the calculation of the linear singular vector and linear singular value. In order to study nonlinear mechanism of the amplification of the initial perturbations, Mu proposed the concept of the nonlinear singular vectors and nonlinear singular values. In Duan et al. [Duan and Mu (2009)], the concept of the conditional nonlinear optimal perturbation (CNOP) was introduced. The CNOP is the initial perturbation whose nonlinear evolution attains the maximal value of the cost function constructed according to the physical problems of interests at a specified time with physical constraint conditions and is called optimal. In Mu et al. [Mu and Zhang (2006)], the CNOP of a two-dimensional quasigeostrophic model is obtained numerically. In Wang et al. [Wang, Zhang and Zhu (2015)], the CNOP is applied to examine the non-linear and linear evolutions of perturbation in stochastic basic flows. The CNOP can be regarded as the most nonlinearly unstable initial perturbation superimposed on the basic state, which usually characterizes the structure of initial errors possessing the largest impact on the uncertainties at the prediction time. There exists a remarkable difference, if the initial perturbation becomes large, or the period is considerably long, or both. This means the TLM may fail to approximate the original nonlinear model. Hence, we focus on the problem of the validity of the TLM by the Lorenz model. Compared with the TLM and the CNOP, a new objective function is introduced to combine the nonlinear mechanism of amplification of initial perturbations and the TLM of the nonlinear model. The relationship between several parameters and the validity of the TLM of the Lorenz model is studied by using the method of mathematical statistics. The paper is organized as follows. In Section 2, the TLM of the Lorenz equations is introduced. In Sections 3 and 4, we analyze the effect of the parameters on the validity of the TLM for the Lorenz equations. In Section 5, the relationship between the limit of the time prediction of the validity of the TLM and the initial perturbation is established. Finally, some conclusions are given in Section 6.

TLM for the Lorenz equations
In this section, we first derive the TLM for the Lorenz equations. Adding perturbation (δx, δy, δz) to the state variables (x, y, z), we get (2) Using Eq. (1), we obtain the TLM for the Lorenz equations below if we omit the higher order terms δxδz and δxδy, where (δx, δy, δz) is the perturbation of (x, y, z).
Corresponding to the definition of the TLM, the accuracy of program of the TLM is defined by the following formula: where M n represents the nonlinear model Eq.
In Richter [Richter (2001)], the initial state is presented as , which is computed since the Lorenz system is known to exhibit chaotic behavior. The parameters have the values σ = 10, r = 28 and b = 8/3, for which the well-known "butterfly" attractor exists.  Figure 1: Results of the TLM program test Tab. 1 and Fig. 1 all show that R consistently approximates to 1 as α approaches 0. It shows that the program of the TLM is correct and effective. By the way, the reason why α is so small that |R − 1| gradually becomes large is the machine error which occurs for converting decimal digit to binary digit in computer.

The validity of the TLM
When the initial perturbation or the time prediction increases, the validity of the TLM may be destroyed [Farrell (1990)]. In this study, we mainly consider the effect of r, T , ∆t and the initial perturbation on the validity of the TLM, with σ = 10 and b = 8/3. In this paper, we choose the Runge-Kutta discrete scheme [Verwer (1996); Guellal, Grimalt and Cherruault (1997); Zingg and Chisholm (1999)]. Similarly, the last two equations are analyzed as the first equation: where P n represents (x n , y n , z n ), P n denotes ( x n , y n , z n ), and F i represents the ith equation of Eq.
(1). Combining Eq. (6) and Eq. (7), we obtain that Expanding Eq. (8) by the T aylor formula, we get that Similarly, we get that In order to get the perturbation of the nonlinear evolution, making Q n ((x + δx) n , (y + δy) n , (z + δz) n ), then, we obtain that Similarly, we obtain that Then, the nonlinear evolution of perturbation can be described as For the TLM, taking the similar discrete format, we obtain that where LF i denotes the ith equation of Eq.
Then, we can obtain the approximate error between nonlinear and linear evolution of perturbation: Two kinds of initial state are chosen in numerical experiments. We set the parameters σ = 10, b = 8/3, r = 28 and Time step ∆t = 0.01. Experiment 1 : the initial state is chosen as (−3.86, −8.77, 21.2). Fig. 2 shows the evolution of error in the L 2 norm between nonlinear model and linear model. Therefore, the result of numerical experiment is close to theoretical analysis. In fact, not only the norm of approximating error, but also the components of state have the same result. In Fig. 3, it is the evolution of δx, δy, δz, approximation between nonlinear and linear, respectively. Experiment 2 : the initial state is (−6 √ 2, −6 √ 2, 27). The evolution of error in the L 2 norm between nonlinear model and linear model is depicted in Fig. 4 where ( theoretical, respectively. The graphic of (a) is the same to the graphic of (b) so that the numerical scheme is effective. Fig. 5 shows the evolution of δx, δy, δz, approximation between nonlinear and linear and also depicts the result of numerical experiment is close to theoretical analysis. The effect of parameters is studied in this section. Some scholars have studied the validity of the TLM by the evolutions of the CNOP or the linear singular vector, which needs two objective functions. Here, we use a new objective function to show the evolution. Furthermore, the validity of the TLM of the initial perturbation ϕ * is called the conditional nonlinear and linear optimal perturbation with constraint ||φ 0 || ≤ σ and the objective function E(φ 0 ), if and only if We say that the TLM is valid as long as E(φ 0 ) ≤ Φ 0 × η, for a given small number η. E(φ 0 ) is the evolution error between the nonlinear and linear evolution and we denote it by E in the following. At the same time, we give a confidence interval if the parameter makes The TLM is said to be valid if the parameters fall in the confidence interval. When r < 1, there is no chaotic phenomenon in the Lorenz model; When r ≥ 1, there may be chaotic phenomenon in the Lorenz model. When 0 ≤ r ≤ 1, there is only one equilibrium point (0, 0, 0), which is asymptotically stable. When 1 ≤ r, there are two additional equilibria: x = y = ± b(r − 1), z = r − 1, which are also asymptotically stable [Tabor and Weiss (1981); Zhou, Liao, Zheng et al. (2004)]. We are more interested in the situation of r ≥ 1. Firstly, we fix the parameters σ = 10, b = 8/3, t = 200 and ∆t = 0.01. Let the initial perturbation equal 0.01 times the initial states (−3.86, −8.77, 21.2) and (−6 √ 2, −6 √ 2, 27). Consider the parameter r = 1, r = 20 and r = 28. In Fig. 6, we plot the evolution error E for r = 1, r = 20 and r = 28. We obtain that the initial state (−3.86, −8.77, 21.2) is closer to the equilibrium point than the initial state (−6 √ 2, −6 √ 2, 27) for r = 1 and r = 20; on the contrary, the initial state (−6 √ 2, −6 √ 2, 27) is closer to the equilibrium point than the initial state (−3.86, −8.77, 21.2) for r = 28. And the error E is relatively small for different r. As the initial state is (−3.86, −8.77, 21.2), using the method of nonlinear least squares, we obtain that E = 10 −5 3r 2 + 49r − 159 |sin(0.52330r − 0.78495) − 0.31831| + 19 .
As the initial state is (−6 √ 2, −6 √ 2, 27), we get that  Fig. 10, E ≤ Φ 0 × 10 −4 and the TLM is valid as ∆t ≤ 0.025, for the initial state (−3.86, −8.77, 21.2). Meanwhile, when E ≤ Φ 0 × 10 −3 , we can obtain that the TLM is valid whatever ∆t is. [0.025, +∞) is the confidence interval. However, Fig. 11 shows that ∆t increases to 0.12, the TLM is valid for the initial state (−6 √ 2,−6 √ 2, 27). When E ≤ Φ 0 × 10 −4 , ∆t ≤ 0.14, the TLM is valid. Therefore, [0.12, 0.14] is the confidence interval for the initial state (−6 √ 2,−6 √ 2, 27). We note that although the sample size is relatively small, the study is not affected. Thirdly, the TLM is valid, when the initial perturbation is small enough in the pioneering work. There are little studies on the numerical relationship between the initial perturbation and the validity of the TLM. Let σ = 10, b = 8/3, r = 28, T = 120 and ∆t = 0.01. The initial perturbation is initial state times 0.1, 0.01 and 0.001.  subgraphs (a, b), the initial perturbation 0.01 in subgraphs (c, d) and the initial perturbation 0.001 in subgraphs (e, f ) In Fig. 12, they are relationship between E and φ 0 , where φ 0 is the L 2 norm of the initial perturbation. We obtain that E becomes smaller, as the initial perturbation is smaller for the initial state (−3.86, −8.77, 21.2) and (−6 √ 2, −6 √ 2, 27). As the initial state is (−3.86, −8.77, 21.2)  Fig. 13 shows that the TLM is valid as φ 0 ≤ 0.02 and E ≤ Φ 0 × 10 −4 . And the TLM is valid as φ 0 ≤ 0.26 and E ≤ Φ 0 × 10 −3 . Therefore [0.02, 0.26] is the confidence interval for the initial state (−3.86, −8.77, 21.2). Similarly, we obtain that [0.3, 0.86] is the confidence interval for the initial state (−6 √ 2, −6 √ 2, 27). The initial perturbation has obviously effect on the error of approximation between nonlinear and linear model. Finally, let σ = −10, b = 8/3, r = 28, ∆t = 0.01, and initial perturbation is (x, y, z) times 0.01. We mainly analyze the effect of T , i.e., T = 120, T = 500, T = 1000. Fig. 15 shows that E becomes larger when T is larger. It means that the nonlinear evolution of the Lorenz model and the TLM is complicated and is hard to predict ,as the terminal time T is large. And at the same time prediction, E is much smaller, when the initial state is close to the equilibrium point.  Figure 17: The dotted line denotes the scatter diagram of E and T and the full line denotes numerical simulation of the relationship between E and T for the initial stata (−6 √ 2, −6 √ 2, 27). The dash line and pecked line denote E = Φ 0 × 10 −4 and E = Φ 0 × 10 −3 , respectively From Fig. 16 and Fig. 17, E produces the periodic changing with T when T becomes small. It is the result of Lorenz model rather than an accidental phenomenon when the terminal time T is limited. In Fig. 16, the TLM is valid when T ≤ 500. And the confidence interval can be omited since E = Φ 0 × 10 −3 almost equals with E = Φ 0 × 10 −4 and E is much greater as T increases. Similarly, T ≤ 2200, the TLM is valid, in Fig. 17. Whatever T is, the TLM is valid in the condition of E ≤ Φ 0 × 10 −3 .

Relationship between initial perturbation and the prediction time of validity of the TLM
The limit of the time prediction of validity of the TLM is defined as the time T p at which The parameters have the values σ = 10, r = 28, b = 8/3, for which the well-known "butterfly" attractor exists. A long integration using the Runge-Kutta method with a time stepsize ∆t = 0.01 was performed to obtain 800 points within the attractor, with two kinds of initial states. A set of random errors, which have a Gaussian distribution with zero mean and the root-mean-square distance is magnitude φ 0 , were superimposed on the initial states to form an ensemble of perturbed states. An ensemble of 800 errors, which M n (Φ 0 + φ 0 ) − M n (Φ 0 ) − M l (φ 0 ) , at each time step was obtained by integrating solutions of the Lorenz model. The ensemble average of the magnitude of the errors was defined as the average error E. The ensemble average was taken as the geometric mean [Ding and Li (2012)].  Fig. 18 shows T p as a function of lg φ 0 , with the initial state (−3.86, −8.77, 21.2). T p decreases approximately linearly with increasing lg φ 0 for φ 0 ≥ 10 −2.6 ; however, T p decreases more quickly with increasing lg φ 0 for φ 0 ≤ 10 −2.6 . Therefore, we conclude that the relationship between T p and φ 0 in the Lorenz model can result from the logarithmic regression using Eq. (19).
The logarithmic relationship implies that T p is unbounded as φ 0 approaches 0. The time of validity of the TLM can be extended as long as we can reduce the initial perturbation.

Conclusions
In this paper, a new objective function M n (Φ 0 + φ 0 ) − M n (Φ 0 ) − M l (φ 0 ) is defined to combine the nonlinear mechanism of amplification of initial perturbations and the TLM of the nonlinear model for studying the validity of the TLM, so that the validity can be described more precisely. We obtained some estimates of parameters numerically about the validity of the TLM for the Lorenz equations. Finally we obtain the confidence interval to expand the scope of the relevant parameters. In fact, the validity of the TLM is affected strongly by the initial states in the Lorenz model. Our result suggests that the limit of the time prediction of validity of the TLM can be extended as long as the initial error can be reduced.