A New Method for Parameter Sensitivity Analysis of Lorenz Equations

Abstract

A new method for parameter sensitivity analysis of Lorenz equations is presented. The sensitivity equations are derived based on the staggered methods. Experimental results indicate that it is possible to determine effects of parameters on model variables so that we can eliminate the less effective ones. Robustness can also be verified in some confidence intervals by simply looking at the corresponding phase portraits. This enables us to control the system. Although the stability properties of the Lorenz equations are studied extensively, to the best knowledge of the authors, the PSA of Lorenz equations has not been considered which is the main goal of this paper.

1. Introduction

Parameter sensitivity analysis (PSA) of large-scale differential algebraic systems is important in many engineering and scientific applications, including biology, chemistry, and economics. Problems such as population dynamics, network modeling, and chemical reactors coming from different branches of science have many parameters whose values may not be known accurately. Infinitesimal changes in most of these model input parameters change the future behavior of the systems partially or completely. Consequently, one can observe an uncontrolled and chaotic behavior of the system. In the present day, one has an opportunity to adjust these parameter values accordingly and make some list of parameters with respect to their effect on model. For instance, if a parameter is less effective than the other parameters, the designer of the model can eliminate that parameter. The analysis of this effectiveness is called parameter sensitivity analysis. Consequently, algorithms which perform PSA in an efficient and rapid manner are invaluable to researchers in many fields.

In this paper, a new method for parameter sensitivity analysis of Lorenz equations is presented. The sensitivity equations are derived based on the staggered methods. Experimental results indicate that it is possible to determine effects of parameters on model variables so that we can eliminate the less effective ones. Robustness can also be verified in some confidence intervals by simply looking at the corresponding phase portraits. This enables us to control the system. Although the stability properties of the Lorenz equations are studied extensively, to the best knowledge of the authors, the PSA of Lorenz equations has not been considered which is the main goal of this paper.

The structure of this paper is as follows. In Section 2, we overview the concept of parameter sensitivity analysis. In Sections 3 and 4, we study the chaotic and sensitivity analysis of Lorenz equations. We complete the paper by some simulation results.

2. Parameter Sensitivity Analysis

It is difficult to construct a model without any parameter. In fact, in reality, the problems coming from different branches of science such as engineering, biology, ecology, and meteorology have many parameters. With the help of faster computers of today, one has a chance to adjust them and make some list of parameters with respect to their effect on model. If a parameter is less effective than the others we, the designer of the model, can eliminate it. The analysis of this effectiveness is called “sensitivity.” When qualitative estimates of sensitivity are desired, a mathematical model of phenomena is desired or at least a relationship is required. Infinitesimal changes in all (or some) of the model input parameters change the future of mathematical design partially or completely (in some cases). The important thing here is the sensitivity of a single component compared to other input variables changed a little bit simultaneously. By saying single component we mean the parameters in the model whose values may not be accurately known. However, such a model brings the questions concerning stability, optimality, sensitivity, and so forth. In this work we concentrated on only the sensitivity analysis of a concrete example, namely, Lorenz equations.

PSA generates essential information for parameter estimation, optimization, control, model simplification, and experimental design. In the literature, staggered direct method, simultaneous corrector method, adjoint method, and staggered corrector method are some of the well-known methods for parameter sensitivity analysis. We can give [1] as a general reference for most of these methods. Some popular software packages for the same task can be listed as ASAP, DASPK, and DASKADJOINT.

In the theory of PSA, another important concept is the index structure which could be defined as the number of differentiations needed for transforming a DAE to an ODE. Intuitively, it is clear that all ODEs have index 0. What defines the index is up to the constraints given in systems. For example, let us consider a simple predator-prey model 𝑦^′𝑥=3𝑥−𝑦,₂+𝑦=𝑔(𝑥), (1) where 𝑔 is differentiable. Taking derivative of constraint equality, we get 𝑦^′=3𝑥−𝑦, and 2𝑥𝑥^′+𝑦^′=𝑔^′(𝑥)⇒𝑥^′=(𝑔^′(𝑥)−3𝑥+𝑦)/2𝑥. Hence, new ODEs are given as 𝑦^′𝑥=3𝑥−𝑦,_′=𝑔^′(𝑥)−3𝑥+𝑦.2𝑥 (2) To obtain this, it takes one differentiation. Thus, the model has index 1.

In order to capture the main idea of the PSA, let us consider the general form of the parameter-dependent DAEs given by 𝐹̇𝑥(0)𝑡,𝑥,𝑥,𝑝=0,=𝑥₀𝑝, (3) where 𝑥∈𝑅^𝑛 and 𝑝∈𝑅^𝑝.

It is not always the case but assume that we have index 0 or 1 DAEs and convert (3) to explicit form of ODEs: ̇,𝑥(0)𝑥=𝑓𝑡,𝑥,𝑝=𝑥₀𝑝, (4) where (𝑥,𝑡,𝑝)∈𝑅^𝑛×𝑅×𝑅^𝑝.

Sensitivity analysis requires the calculation of the term, namely, 𝑠_𝑖, defined as the derivative of 𝑥 with respect to 𝑝_𝑖; that is, 𝑠_𝑖∶=𝜕𝑥/𝜕𝑝_𝑖. Since we are interested in partial derivatives, we can treat one parameter after another, while keeping the remaining ones fixed. Therefore, the derivative of (4) with respect to parameter 𝑝_𝑖 is 𝜕𝜕𝑝_𝑖𝜕𝑥=𝜕𝜕𝑡𝜕𝑝_𝑖𝑓⟺𝜕𝑡,𝑥,𝑝𝜕𝑡𝜕𝑥𝜕𝑝_𝑖=𝜕𝑓𝜕𝑝_𝑖+𝜕𝑓𝜕𝑥𝜕𝑥𝜕𝑝_𝑖. (5) Replacing 𝑠_𝑖=𝜕𝑥/𝜕𝑝_𝑖 into the right-hand side of (5), the 𝑖th sensitivity equation becomes 𝑠_𝑖̇=𝑓_𝑝𝑖+𝑓_𝑥𝑠_𝑖. (6) Since we have 𝑠_𝑖=𝜕𝑥/𝜕𝑝_𝑖=[𝜕𝑥₁/𝜕𝑝_𝑖,𝜕𝑥₂/𝜕𝑝_𝑖,…,𝜕𝑥_𝑛/𝜕𝑝_𝑖]^𝑇=[𝑠_𝑖1,𝑠_𝑖2,…,𝑠_𝑖𝑛]^𝑇, the previous process introduces 𝑛⋅𝑝 additional differential equations. Finally, defining a new variable 𝑥𝑠𝑧∶=, (7) where 𝑠=𝑠₁,𝑠₂,…,𝑠_𝑝=𝑠₁₁,𝑠₁₂,…,𝑠_1𝑛,𝑠₂₁,𝑠₂₂,…,𝑠_2𝑛,…,𝑠_𝑝1,𝑠_𝑝2,…,𝑠_𝑝𝑛, (8) we write ̇𝑥̇𝑠̇=𝑓𝑧=𝑡,𝑥,𝑝𝜕𝑓𝑑𝑝_𝑖+𝜕𝑓𝜕𝑝_𝑖𝑠=𝑓𝑡,𝑥,𝑝𝜕𝑓𝜕𝑝₁+𝜕𝑓𝜕𝑝₁𝑠₁𝜕𝑓𝜕𝑝₂+𝜕𝑓𝜕𝑝₂𝑠₂⋮𝜕𝑓𝜕𝑝_𝑝+𝜕𝑓𝑠𝜕𝑝^𝑝. (9) The initial condition takes the form ∀𝑖,𝑠_𝑖(0)=𝜕𝑥₀𝜕𝑝_𝑖(0)=𝑥⟹𝑧𝜕𝑥⁰₀𝜕𝑝_𝑖=𝑥₀𝑠₁₁(0)𝑠₂₁(0)𝑠₃₁(0)⋮𝑠_𝑛1(0)⋮𝑠_1𝑝(0)𝑠_2𝑝(0)⋮𝑠_𝑛𝑝(0). (10)

3. Lorenz Equations

The Lorenz equations invented by E. N. Lorenz, a meteorologist and a pioneer of chaos theory, are typical examples of equations for system of differential algebraic equations that can be written as 𝑋^′𝑌^′𝑍^′=𝜎(𝑌−𝑋)𝑟𝑋−𝑌−𝑋𝑍−𝑏𝑍+𝑋𝑌, (11) where 𝑟>0 is Rayleigh number, 𝜎 is Prandtl number corresponding to temperature difference between two horizontal plates in convection problem, and 𝑏 is a positive number. These equations arise in studies of convection and instability in planetary atmospheres, models of lasers and dynamos, and so forth. Although the stability and bifurcation properties of the Lorenz equations are studied in the literature [2, 3], to the best knowledge of the authors, the parameter sensitivity analysis of Lorenz equations has not been considered so far which is the main goal of this paper.

The Lorenz equations are nonlinear due to the terms 𝑥𝑦 and 𝑥𝑧. They are also symmetric equations, because the equations are invariant under (𝑥,𝑦)→(−𝑥,−𝑦). Thus, if (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) is a solution of Lorenz equations, so is (−𝑥(𝑡),−𝑦(𝑡),𝑧(𝑡)). System of the Lorenz equations is dissipative. In other words, volumes in phase-space contract under the flow and 𝜎 and 𝑏 are usually known as dissipation parameters. Next, we compute the fixed points of Lorenz equations. Letting each term of (11) be equal to 0, that is, letting 𝑥̇=𝑦̇=𝑧̇=0, we get the following identities: 𝑥=𝑦,𝑟𝑥−𝑦−𝑥𝑧=0,𝑥𝑦=𝑏𝑧, (12) which implies that (0,0,0) is a fixed point of the Lorenz equations. If 𝑥=𝑦≠0, then from the appropriate ones of the previous equations we obtain 𝑧=𝑥²/𝑏,  𝑟𝑥−𝑥−𝑥(𝑥²/𝑏)=0 which gives us 𝑥²=𝑏(𝑟−1) by bearing in mind that 𝑥≠0. If 𝑟>1, then 𝑥^∗=𝑦^∗=±𝑏(𝑟−1) and 𝑧^∗=𝑟−1 that are the fixed-points of Lorenz equations. While 𝑟→1⁺, these equations generate a so-called pitchfork bifurcation. Now, we are making behavior analysis of these fixed points. Depending on the parameter 𝑟, we have three critical points; namely, 𝑎𝑥=𝑦=𝑧=0,₁∶𝑥=𝑦=+[𝑏(𝑟−1)]^1/2,𝑎₁∶𝑥=𝑦=−[𝑏(𝑟−1)]^1/2. (13)

In this study, behaviors dependent on initial conditions are not studied, and they are fixed as 𝑥₀=0, 𝑦₀=1, and 𝑧₀=0. In the following figures different trajectories are given with respect to different 𝑟 values.

The behavior in Figure 1 continues up to a value of 𝑟=24.08. After that, it becomes more complicated and chaotic; for example, for 𝑟=27 some periodic and aperiodic motions are observed as seen in Figure 2.

(a)

(b)

Figure 1 

For 𝑟=3, projection of three-dimensional phase space onto 𝑧-𝑦 space (a) and corresponding solutions (b).

(a)

(b)

(c)

(d)

Figure 2 

For 𝑟=27. The 𝑥 in (a) and 𝑦 in (b) are periodic, whereas 𝑧 in (c) is aperiodic. (d) 𝑦𝑧 space for critical values 𝑎1 and 𝑎2.

Further explanations of these and stability features of Lorenz equations might be seen, for instance, at [4]. In the next section, we study parameter sensitivity analysis of Lorenz equations.

4. Parameter Sensitivity Analysis of Lorenz Equations

Let us write the Lorenz equations having some initial conditions in the following way: 𝑥̇₁𝑥̇₂𝑥̇₃=𝑝₁𝑥₂−𝑥₁𝑝₂𝑥₁−𝑥₂−𝑥₁𝑥₃−𝑝₃𝑥₃+𝑥₁𝑥₂𝑥₁(0)𝑥=0,₂(0)𝑥=1,₃(0)=0. (14)

Our new variable 𝑧 defined in (7) is given as follows: 𝑥𝑧=¹𝑥₂𝑥₃𝑠₁𝑠₂𝑠₃, (15) where 𝑠₁=𝑑𝑥𝑑𝑝₁=𝑑𝑥₁𝑑𝑝₁𝑑𝑥₂𝑑𝑝₁𝑑𝑥₃𝑑𝑝₁^𝑇=𝑠₁₁𝑠₁₂𝑠₁₃^𝑇𝑠₂=𝑑𝑥𝑑𝑝₂=𝑑𝑥₁𝑑𝑝₂𝑑𝑥₂𝑑𝑝₂𝑑𝑥₃𝑑𝑝₂^𝑇=𝑠₂₁𝑠₂₂𝑠₂₃^𝑇𝑠₁=𝑑𝑥𝑑𝑝₃=𝑑𝑥₁𝑑𝑝₃𝑑𝑥₂𝑑𝑝₃𝑑𝑥₃𝑑𝑝₃^𝑇=𝑠₃₁𝑠₃₂𝑠₃₃^𝑇. (16)

Finally, the sensitivity equations take the form [𝑧]̇=𝑥̇₁𝑥̇₂𝑥̇₃𝑠̇₁₁𝑠̇₁₂𝑠̇₁₃𝑠̇₂₁𝑠̇₂₂𝑠̇₂₃𝑠̇₃₁𝑠̇₃₂𝑠̇₃₃=𝑝₁𝑥₂−𝑥₁𝑝₂𝑥₁−𝑥₂−𝑥₁𝑥₃−𝑝₃𝑥₃+𝑥₁𝑥₂𝑥₂−𝑥₁+𝑝₁𝑠₁₂−𝑠₁₁𝑝₂𝑠₁₁−𝑠₁₂−𝑠₁₁𝑥₃−𝑥₁𝑠₁₃−𝑝₃𝑠₁₃+𝑠₁₁𝑥₂+𝑥₁𝑠₁₂𝑝₁𝑠₂₂−𝑠₂₁𝑥₁+𝑝₂𝑠₂₁−𝑠₂₂−𝑠₂₁𝑥₃−𝑥₁𝑠₂₃−𝑝₃𝑠₂₃+𝑠₂₁𝑥₂+𝑥₁𝑠₂₂𝑝₁𝑠₃₂−𝑠₃₁𝑝₂𝑠₃₁−𝑠₃₂−𝑠₃₁𝑥₃−𝑥₁𝑠₃₃−𝑥₃−𝑝₃𝑠₃₃+𝑠₃₁𝑥₂+𝑥₁𝑠₃₂. (17)

Note that the initial conditions do not depend on parameters, so the new initial conditions 𝑠_𝑖𝑗=0  for all 𝑖,𝑗=1,2,3. In other words, the initial conditions for the new variable are given as 𝑧₀=[0,1,0,0,0,0,0,0,0,0,0,0]^𝑇. (18)

In the next section we are making some simulations in order to visualize the results.

5. Computational Results

In order to visualize the results, we made many different simulations as phase portraits, sensitivity analysis, and relation between the components. In this section, we present only some of the simulation results.

After integrating the sensitivity equations, we get 𝑠_𝑖’s as a function of time so that one can analyze the change in the solution with respect to perturbations in the parameters. In these experiments, for time interval 𝑡=[0100], we solved the system by Matlab ODE solver, namely, ode45, based on an explicit Runge-Kutta Method. After integrating the sensitivity equations, we get 𝑠_𝑖’s as a function of time so that one can analyze the sensitivity in phase portraits.

For 𝑟=30, the qualitative behaviors of the system might be seen in Figure 3.

Figure 3 

For 𝑟=30, the behavior of the systems at 𝑥-𝑧 plane.

Before entering the chaotic region which starts from the value 𝑟=24.08, all nine sensitivity components demonstrate the same behavior. For example, taking 𝑟=20, we have well stable solutions, but the relation between the sensitivity variables is highly nonlinear, that is, a significantly important result for this well-known system. This is illustrated in Figure 4.

(a)

(b)

Figure 4 

(a) 𝑠13 and 𝑠31 components. (b) Phase portrait of 𝑠21 and 𝑠31.

Remember that 𝑠₂₁ and 𝑠₃₁ represent sensitivity components of 𝑥₁ with respect to 𝑝₂ and 𝑝₃. Figure 4(b) tells us that altering 𝑝₂ and 𝑝₃ in a simultaneous manner can affect the controllability of the system completely. In Figure 5, phase portrait for 𝑠₃₁ and 𝑥₁ is given.

Figure 5 

Phase portrait for 𝑠31 and 𝑥1.

In the chaotic region, what happened to the sensitivity equations is that, first, they seemed to be very complicated and irregular when they are considered as a function of time as seen in Figure 6.

(a)

(b)

Figure 6 

Solutions of sensitivity equations: (a) 𝑠22 and (b) 𝑠13.

However, in the phase portraits of 𝑠₂₂ and 𝑠₁₃, we obtained a completely linear relation between them as seen in Figure 7.

(a)

(b)

Figure 7 

6. Conclusion

When qualitative estimates of sensitivity are desired, a mathematical model of phenomena is desired or at least a relationship is required. However, such a model brings the questions concerning stability, optimality, sensitivity, and so forth. In this work, we concentrated only on the PSA of Lorenz equations. As we saw in the small application, it is possible to determine effects of parameters on model variables so that we can eliminate the less effective ones. Robustness can also be verified in some confidence intervals by just looking at the phase portraits. This enables us to control the system. This method is efficient if the number of variables is much more than that of the parameters. In a future work, we plan to study PSA for Van der Pool equations.