The Time Delays’ Effects on the Qualitative Behavior of an Economic Growth Model

Abstract

A further generalization of an economic growth model is the main topic of this paper. The paper specifically analyzes the effects on the asymptotic dynamics of the Solow model when two time delays are inserted: the time employed in order that the capital is used for production and the necessary time so that the capital is depreciated. The existence of a unique nontrivial positive steady state of the generalized model is proved and sufficient conditions for the asymptotic stability are established. Moreover, the existence of a Hopf bifurcation is proved and, by using the normal form theory and center manifold argument, the explicit formulas which determine the stability, direction, and period of bifurcating periodic solutions are obtained. Finally, numerical simulations are performed for supporting the analytical results.

1. Introduction

Most of the phenomena occurring in real-world complex systems, especially in the economics systems, have not an immediate effect but appear with some delay. Therefore time delays have been inserted into mathematical models and in particular in models of the applied sciences based on ordinary differential equations; see the recent book [1]. Differential equations specifically with time delays have been proposed in population dynamics [2] for biological systems such as immune system response [3–6] and tumor growth [7–12], in models of social sciences [13], and in economics systems; see, among others, [14–21].

The introduction of a time delay into an ordinary differential equation could change the stability of the equilibrium (stable equilibrium becomes unstable) and could cause fluctuations, and Hopf bifurcation can occur. Indeed global existence of Hopf bifurcations has been proved in many delay mathematical models; see papers [22–24] and references cited therein.

If on one hand, the stability and bifurcation analysis of ordinary differential equations with a single time delay is well outlined in the pertinent literature [25, 26], on the other hand the analysis of the dynamics of ordinary differential equations with multiple time delays is a difficult task [27–29] and the related literature is much limited. In this context, for several classes of ordinary differential equation models with multiple time delays, sufficient and necessary conditions have been established and a complete description of the stability region has been reached; see, among others, [30, 31] and references cited therein.

The present paper is concerned with a further generalization of the Solow model [32]. The generalized model is governed by a delay differential equation with two time delays. Specifically the two time delays refer, respectively, to the time employed in order that the capital is used for production and the necessary time so that the capital is depreciated. The asymptotic analysis performed in this paper shows the existence of a unique nontrivial positive steady state, and sufficient conditions for the asymptotic stability are established. Moreover, the existence of a Hopf bifurcation is proved and, by using the normal form theory and center manifold argument, the explicit formulas which determine the stability, direction, and period of bifurcating periodic solutions are obtained. Finally, some numerical simulations to support the analytical conclusions are carried out.

The rest of the paper is organized into four more sections which follow this introduction. Specifically, Section 2 discusses the derivation of the generalized Solow model with two time delays. Section 3 deals with the analysis of the existence of Hopf bifurcation and stability of the positive equilibrium for the model proposed in Section 1. Section 4 is concerned with the direction and the stability of the Hopf bifurcation. Some numerical simulations are performed in Section 5 with the aim of supporting the analytic results. Finally Section 6 completes the paper with conclusions and some future research perspectives.

2. The Mathematical Model

Recently, Zak has proposed in [15] the following delay ordinary differential equation 𝑘̇(𝑡)=𝑠𝑓(𝑘(𝑡−𝜏))−𝛿𝑘(𝑡−𝜏) (1) to describe the dynamics of the Solow model [32] in which a production technology has a constant finite period 𝜏≥0 linked to the time needed for the installation of capital. Here, 𝑘(𝑡) denotes capital at time 𝑡, 𝑓(𝑘(𝑡)) is a neoclassical production function, namely a function which is continuous, increasing, and strictly concave in capital, and 𝑠∈(0,1) is the constant savings rate. Population is assumed to be constant and normalized to unity. During production a proportion of the capital stock, 𝛿∈(0,1), depreciates at the same gestation period 𝜏.

The assumption that the growth of the amount of capital at time 𝑡 is a function of the total output of capital at time 𝑡−𝜏 has been showed in [15] to be the source of cyclic behaviour in the economic system (1). Generalizing this idea, we proposed the following delay ordinary differential equation: 𝑘̇(𝑡)=𝑠𝑓𝑘𝑡−𝜏₁−𝛿𝑘𝑡−𝜏₂, (2) that is, the generalized Solow model with two delays.

It is worth stressing that for 𝜏₁=𝜏₂=𝜏=0, we recover the Solow model equation [32]. In this model, the positive equilibrium is asymptotically stable in the absence of delay. For 𝜏₁=𝜏₂=𝜏>0, (2) reduces to the delayed Solow model proposed by Zak in [15].

3. Local Stability and Hopf Bifurcation

The mathematical model (2) has exactly the same equilibrium points of the corresponding system with zero delays. Hence, there exists a unique positive equilibrium 𝑘_∗, where 𝑠𝑓(𝑘_∗)=𝛿𝑘_∗. To determine the stability of this equilibrium and Hopf bifurcation, we linearize (2) around 𝑘_∗. The result is a linear delay differential equation of the form 𝑘̇(𝑡)=−𝐴₁𝑘𝑡−𝜏₁−𝑘_∗−𝐴₂𝑘𝑡−𝜏₂−𝑘_∗, (3) where 𝐴₁=𝑠𝑓^′𝑘_∗−𝛿<0,𝐴₂=𝛿>0. (4) It is well known that the stability of the equilibrium is determined by the spectrum of the eigenvalues of the linearization, which can be found as the roots of the characteristic equation 𝜆=−𝐴₁𝑒^−𝜆𝜏1−𝐴₂𝑒^−𝜆𝜏2. (5) We recall that an equilibrium point of an equation is stable if all eigenvalues of its linearization have negative real part. It changes its stability type when eigenvalues cross the imaginary axis of the complex plane. We first note that 𝜆=0 is not a root of (5) because this would imply 𝑠𝑓^′(𝑘_∗)=0, contradicting the fact that 𝑓^′(𝑘_∗)>0. Next, the distribution of the roots of (5) should be investigated. However, the analysis of the sign of the real parts of eigenvalues is very complicated because of the presence of two different delays, 𝜆₁ and 𝜆₂, in (5). Therefore, we will use a method consisting of determining the stability of the equilibrium when one delay is equal to zero, and, using similar analytic arguments as in the work by Ruan and Wei [33], we will deduce conditions for the stability of the equilibrium when both time delays are nonzero.

3.1. Case 1: 𝜏₁≥0 and 𝜏₂=0

The characteristic equation (5) reduces to 𝜆=−𝐴₁𝑒^−𝜆𝜏1−𝐴₂. (6) If 𝜏₁=0, (6) has the unique root 𝜆=−𝐴₁−𝐴₂=−𝑠𝑓^′(𝑘_∗)<0. Thus, the equilibrium 𝑘_∗ is locally asymptotically stable. Consequently, when 𝜏₁ increases, the stability of the steady state can only be lost if pure imaginary roots appear. Hence we look for purely imaginary roots ±𝑖𝜔, 𝜔>0, of (6). Let 𝑖𝜔 be a purely imaginary root of (6). Then, separating real and imaginary parts, 𝜔 satisfies 𝜔=𝐴₁sin𝜔𝜏₁,𝐴₂=−𝐴₁cos𝜔𝜏₁. (7) It follows that 𝜔²=𝐴²₁−𝐴²₂=𝐴₁+𝐴₂𝐴₁−𝐴₂=𝑠𝑓^′𝑘_∗𝑠𝑓^′𝑘_∗−𝛿−𝛿<0. (8) Hence, (6) has no positive root. Then, we can conclude the following result about the asymptotic stability of the equilibrium of (2).

Proposition 1. Let 𝜏₂=0. Then the positive equilibrium 𝑘_∗ of (2) is locally asymptotically stable.

3.2. Case 2: 𝜏₁=0 and 𝜏₂≥0

The characteristic equation (5) becomes 𝜆=−𝐴₁−𝐴₂𝑒^−𝜆𝜏2. (9)

Setting 𝜏₂=0, we know that the equilibrium 𝑘_∗ is locally asymptotically stable. Let 𝜆=𝑖𝜔, 𝜔>0, be a root of (9). Separating real and imaginary parts, we have the following two equations: 𝜔=𝐴₂sin𝜔𝜏₂,𝐴₁=𝐴₂cos𝜔𝜏₂. (10) Adding the squares of both hand sides of (10), it follows that 𝜔 must be a root of the following equation: 𝜔²=𝐴²₂−𝐴²₁=−𝐴₁−𝐴₂𝐴₁−𝐴₂=−𝑠𝑓^′𝑘_∗𝑠𝑓^′𝑘_∗−𝛿−𝛿>0. (11) Hence, (11) admits only the positive root 𝜔₀=𝐴²₂−𝐴²₁. (12) Define from (10) 𝜏⁰₂=1𝜔₀𝐴arccos¹𝐴₂, (13) and 𝜏^𝑛₂=𝜏⁰₂+2𝑛𝜋𝜔₀,𝑛=0,1,2,…. (14) Then, it is immediate to check that 𝑖𝜔₀ is a simple root of (9) when 𝜏₂=𝜏^𝑛₂.

Lemma 2. Let 𝜆(𝜏₂)=𝑎(𝜏₂)+𝑖𝜔(𝜏₂) be the complex root of (9) near 𝜏₂=𝜏^𝑛₂ satisfying 𝑎(𝜏^𝑛₂)=0 and 𝜔(𝜏^𝑛₂)=𝜔₀. Then the transversal condition 𝜆𝜏sign𝑑Re^𝑛₂𝑑𝜏₂>0 (15) is satisfied.

Proof. Differentiating the characteristic equation (9) with respect to 𝜏₂, we obtain 𝑑𝜆𝑑𝜏₂=𝐴₂𝑒^−𝜏2𝜆𝜏₂𝑑𝜆𝑑𝜏₂+𝜆. (16) This gives 𝑑𝜆𝑑𝜏₂⁻¹1=−𝜆𝜆+𝐴₁−𝜏₂𝜆. (17) Hence, we have 𝑑()signRe𝜆𝑑𝜏_{2𝜆=𝑖𝜔0,𝜏2=𝜏𝑛2}=signRe𝑑𝜆𝑑𝜏₂⁻¹_{𝜆=𝑖𝜔0,𝜏2=𝜏𝑛2}1=sign𝜔²₀+𝐴²₁. (18) This completes the proof.

Bearing the above analysis in mind and the Hopf bifurcation theorem for functional differential equations due to Hale and Verduyn Lunel (see p. 246, Theorem 1.1 of the book [34]), we have the following result.

Theorem 3. Let 𝜏₁=0. Then the positive equilibrium 𝑘_∗ of (2) is locally asymptotically stable when 𝜏₂∈[0,𝜏⁰₂) and unstable when 𝜏₂>𝜏⁰₂. Furthermore, (2) undergoes a Hopf bifurcation at the positive equilibrium 𝑘_∗ when 𝜏₂=𝜏^𝑛₂, 𝑛=0,1,2,…, where 𝜏^𝑛₂ is defined as in (14).

3.3. Case 3: 𝜏₁>0 and 𝜏₂>0

Equation (5) has purely imaginary roots ±𝑖𝜔, where 𝜔>0, if the following equations are satisfied. Consider 𝜔−𝐴₁sin𝜔𝜏₁=𝐴₂sin𝜔𝜏₂,𝐴₁cos𝜔𝜏₁=−𝐴₂cos𝜔𝜏₂. (19) Squaring and adding up both equations in (19) yield sin𝜔𝜏₁=𝑔(𝜔), (20) where 𝑔(𝜔)=𝜔²+𝐴²₁−𝐴²₂2𝜔𝐴₁. (21)

Lemma 4. For every arbitrary 𝜏₁>0, (20) has a finite number of positive solutions for 𝜔∈[𝐴₁+𝐴₂,−𝐴₁+𝐴₂].

Proof. The inequality |sin(𝜔𝜏₁)|≤1 and 𝐴>1 imply 𝜔∈[𝐴₁+𝐴₂,−𝐴₁+𝐴₂]. The function 𝑔(𝜔) has the properties 𝑔(𝐴₁+𝐴₂)=1, 𝑔(−𝐴₁+𝐴₂)=−1, 𝑔(𝐴²₂−𝐴²₁)=0, and 𝑑𝑔(𝜔)=𝜔𝑑𝜔²+𝐴²₂−𝐴²₁2𝜔²𝐴₁𝑑<0,²𝑔(𝜔)𝑑𝜔²=𝐴²₁−𝐴²₂𝜔³𝐴²₁<0. (22) A graphical inspection on the intersections of the functions sin(𝜔𝜏₁) and 𝑔(𝜔) gives the statement.

Remark 5. If 𝜏₁∈(0,𝜋/(−𝐴₁+𝐴₂)), (20) has only one positive solution.

For any 𝜏₁>0, (20) has a finite number of positive zeros 𝜔_𝑗, 𝑗=1,2,…,𝑚. It is clear that for every arbitrary chosen 𝜏₁>0 and for each 𝜔_𝑗 we have an infinite number of 𝜏₂ such that 𝐴₁cos(𝜔_𝑗𝜏₁)=−𝐴₂cos(𝜔_𝑗𝜏₂). For all 𝑗=1,2,…,𝑚, we define 𝜏^𝑗₂=min𝜏₂>0∶𝐴₁cos𝜔_𝑗𝜏₁=−𝐴₂cos𝜔_𝑗𝜏₂. (23) In addition, we set 𝜏⁰₂=min𝜏^𝑗₂∶𝑗=1,2,…,𝑚 (24) and denote 𝜔₀=𝜔_𝑗 for 𝑗 such that 𝜏^𝑗₂=𝜏⁰₂. Next, we check the condition which guarantees that the purely imaginary roots pass through the imaginary axis at ±𝑖𝜔₀. Let 𝜆(𝜏₂)=𝑎(𝜏₂)+𝑖𝜔(𝜏₂) be the root of (5) near 𝜏₂=𝜏⁰₂ such that 𝑎(𝜏⁰₂)=0 and 𝜔(𝜏⁰₂)=𝜔₀. By direct computation we have 𝑑𝜆𝑑𝜏₂=𝐴₁𝑒^−𝜆𝜏1𝜏₁𝑑𝜆𝑑𝜏₂+𝐴₂𝑒^−𝜆𝜏2𝜏₂𝑑𝜆𝑑𝜏₂+𝜆 (25) and obtain 𝑑𝜆𝑑𝜏₂⁻¹=−1+𝜆𝜏₁𝜆𝐴₂𝑒^−𝜆𝜏2+𝜏₁−𝜏₂𝜆. (26) Hence, it follows that 𝑑()signRe𝜆𝑑𝜏_{2𝜆=𝑖𝜔0,𝜏2=𝜏02}=signRe𝑑𝜆𝑑𝜏₂⁻¹_{𝜆=𝑖𝜔0,𝜏2=𝜏02}𝜔=signsin₀𝜏⁰₂+𝜏₁𝜔₀cos𝜔₀𝜏⁰₂𝐴₂𝜔₀=signsin𝜔₀𝜏⁰₂+𝜏₁𝜔₀cos𝜔₀𝜏⁰₂. (27) If for given 𝜏₁ both sin(𝜔₀𝜏⁰₂) and cos(𝜔₀𝜏⁰₂) are positive (resp., negative), then signsin𝜔₀𝜏⁰₂+𝜏₁𝜔₀cos𝜔₀𝜏⁰₂>0resp.,<0 (28) and the purely imaginary roots ±𝑖𝜔₀ of (5) move to the right (resp., left) half plane when the bifurcation parameter 𝜏₂ increases. So, we have the following transversality condition.

Proposition 6. Let Ω₁=[0,𝜋/2] and Ω₂=[𝜋,3𝜋/2].(1)If 𝜔₀𝜏⁰₂∈Ω₁, then sign[𝑑(Re𝜆)/𝑑𝜏₂]_{𝜆=𝑖𝜔0,𝜏2=𝜏02}>0.(2)If 𝜔₀𝜏⁰₂∈Ω₂, then sign[𝑑(Re𝜆)/𝑑𝜏₂]_{𝜆=𝑖𝜔0,𝜏2=𝜏02}<0.

Lemma 7. 𝜆=±𝑖𝜔₀ are simple roots of (5) when 𝜔₀𝜏⁰₂∈Ω₁∪Ω₂.

Proof. If 𝜆=𝑖𝜔₀ is a repeated root for (5), then 1=𝜏₁𝐴₁𝑒^{−𝑖𝜔0𝜏1}+𝜏⁰₂𝐴₂𝑒^{−𝑖𝜔0𝜏02} holds true. Using (5), this leads to −𝜔₀𝜏₁=(𝜏⁰₂−𝜏₁)𝐴₂sin(𝜔₀𝜏⁰₂) and 1=(𝜏⁰₂−𝜏₁)𝐴₂cos(𝜔₀𝜏⁰₂). Thus, we must have −𝜔₀𝜏₁=tan(𝜔₀𝜏⁰₂). If 𝜔₀𝜏⁰₂∈Ω₁, this means 𝜔₀=0, while if 𝜔₀𝜏⁰₂∈Ω₂, this identity does not hold. The conclusion is immediate.

From the discussion above, and recalling that for any 𝜏₁ and 𝜏₂=0 all roots of (5) have strictly negative real parts, the following theorem about stability and Hopf bifurcation of (2) is immediately obtained.

Theorem 8. Let 𝜏₁>0 and Ω₁, Ω₂ be defined as in Proposition 6. Then(1)equation (2) undergoes a Hopf bifurcation at 𝑘_∗ when 𝜏₂=𝜏⁰₂,(2)if 𝜔₀𝜏⁰₂∈Ω₁, then the nontrivial equilibrium 𝑘_∗ to (2) is locally asymptotically stable for 𝜏₂∈[0,𝜏⁰₂) and unstable for 𝜏₂>𝜏⁰₂,(3)if 𝜔₀𝜏⁰₂∈Ω₂, then the nontrivial equilibrium 𝑘_∗ to (2) is locally asymptotically stable for 𝜏₂∈[0,+∞).

4. Direction and Stability of the Hopf Bifurcation

In this section, we study the direction of bifurcations and the stability of bifurcating periodic solutions of (2) at 𝜏⁰₂ by using the method based on the normal form theory and center manifold theory introduced by Hassard et al. [35].

For notational convenience, let 𝜏₂=𝜏⁰₂+𝜇, 𝜇∈ℝ. Then 𝜇=0 is the Hopf bifurcation point for (2). First, we transform (2) into a functional differential equation in 𝐶([−𝜏₁,0],ℝ), which is the Banach space of continuous real-valued functions that map [−𝜏₁,0] into ℝ, and endowed with the norm 𝜑=sup_{−𝜏1≤𝜗≤0}𝜑(𝜗). (29) Set 𝑥(𝑡)=𝑘(𝑡)−𝑘_∗. Then, rewriting (2) in terms of 𝑥(𝑡) and considering its Taylor expansion at the trivial equilibrium up to the third order, we get 𝑥̇(𝑡)=−𝑠𝑓^′𝑘_∗−𝛿𝑥𝑡−𝜏₁−𝛿𝑥𝑡−𝜏₂+𝑠𝑓^′′𝑘_∗2𝑥𝑡−𝜏₁²+𝑠𝑓^′′′𝑘_∗𝑥3!𝑡−𝜏₁³+⋯. (30) For 𝜑∈𝐶([−𝜏₁,0],ℝ), define the linear operator 𝐿_𝜇𝜑=−𝑠𝑓^′𝑘_∗−𝛿𝜑−𝜏₁−𝛿𝜑−𝜏₂ (31) and the nonlinear operator 𝐹𝜇,𝜑=𝑠𝑓^′′𝑘_∗2𝜑−𝜏₁²+𝑠𝑓^′′′𝑘_∗𝜑3!−𝜏₁³+⋯. (32) By the Riesz representation theorem, there exists a bounded variation function 𝜂(𝜗,𝜇) with 𝜗∈[−𝜏₁,0] such that 𝐿_𝜇𝜑=⁰_−𝜏1𝑑𝜂𝜗,𝜇𝜑(𝜗), (33) with 𝜂=(𝜗)𝜗,𝜇−𝛿Γ,𝜗∈−𝜏₂,0,𝑠𝑓^′𝑘_∗−𝛿Γ𝜗+𝜏₁,𝜗∈−𝜏₁,−𝜏₂, (34) where Γ is the Dirac delta function. For 𝜑∈𝐶([−𝜏₁,0],ℝ), define 𝐴𝜇(𝜗)𝜑=𝑑𝜑𝑑𝜗,𝜗∈−𝜏₁,0_0−𝜏1,𝜑(𝑢)𝑅𝜇𝑑𝜂𝑢,𝜇,𝜗=0,𝜑=0,𝜗∈−𝜏₁,𝐹,0𝜇,𝜑,𝜗=0. (35) Then (30) is equivalent to 𝑥̇_𝑡=𝐴𝜇𝑥_𝑡+𝑅𝜇𝑥_𝑡, (36) where 𝑥_𝑡(𝜗)=𝑥(𝑡+𝜗) for 𝜗∈[−𝜏₁,0]. For 𝜓∈𝐶([0,𝜏₁],ℝ), define the operator 𝐴^∗ as 𝐴^∗(0)𝜓(𝑢)=−𝑑𝜓(𝑢)𝑑𝑢,𝑢∈0,𝜏₁,_0−𝜏1𝑑𝜂(𝑟,0)𝜓(−𝑟),𝑢=0, (37) and the bilinear inner product 𝜓,𝜑=𝜓(0)𝜑(0)−⁰_{𝜗=−𝜏1}^𝜗_𝜉=0𝜓(𝜉−𝜗)𝑑𝜂(𝜗,0)𝜑(𝜉)𝑑𝜉, (38) where the over bar denotes complex conjugation. Then 𝐴(0) and 𝐴^∗(0) are adjoint operators. By the discussion in the previous section, we know that ±𝑖𝜔₀ are eigenvalues of 𝐴(0). Thus, they are also eigenvalues of 𝐴^∗(0). We need to compute the eigenvector of 𝐴(0) and 𝐴^∗(0) corresponding to 𝑖𝜔₀ and −𝑖𝜔₀, respectively. A direct computation shows that their eigenvectors are 𝑞(𝜗)=𝑒^{𝑖𝜔0𝜗},𝜗∈−𝜏₁,𝑞,0_∗(𝑢)=𝐵𝑒^{𝑖𝜔0𝑢},𝑢∈0,𝜏₁, (39) respectively. We have ⟨𝑞^∗,𝑞⟩=0. In order to ensure ⟨𝑞^∗,𝑞⟩=1, we choose 𝐵 as 1𝐵=1−𝜏₁𝑠𝑓^′𝑘_∗−𝛿𝑒^{𝑖𝜔0𝜏1}−𝜏⁰₂𝛿𝑒^{𝑖𝜔0𝜏02}. (40) Next, we compute the coordinates to describe the center manifold 𝒞 at 𝜇=0. Let 𝑥_𝑡 be the solution of (36) when 𝜇=0. Define 𝑧(𝑡)=𝑞^∗,𝑥_𝑡,𝑊(𝑡,𝜗)=𝑥_𝑡(𝜗)−2Re𝑧(𝑡)𝑞(𝜗). (41) On the center manifold 𝒞, we have 𝑊(𝑡,𝜗)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜗), with 𝑊(𝑧,𝑧,𝜗)=𝑊₂₀(𝜗)𝑧²2+𝑊₁₁(𝜗)𝑧𝑧+𝑊₀₂(𝜗)𝑧²2+⋯, (42) where 𝑧 and 𝑧 are local coordinates for 𝒞 in the direction of 𝑞^∗ and 𝑞^∗. For any 𝑥_𝑡∈𝒞, since 𝜇=0, we find 𝑧(̇𝑡)=𝑖𝜔₀𝑧(𝑡)+𝑞^∗(0)𝐹₀(𝑧,𝑧), (43) with 𝐹₀(𝑧,𝑧)=𝐹0,𝑊(𝑧,𝑧,0)+2Re𝑧(𝑡)𝑞(0). (44) We rewrite (43) as 𝑧(̇𝑡)=𝑖𝜔₀𝑧(𝑡)+𝑔(𝑧,𝑧), (45) where 𝑔(𝑧,𝑧)=𝑞^∗(0)𝐹0,𝑊(𝑧,𝑧,0)+2Re𝑧(𝑡)𝑞(0)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (46) It follows from (36) and (43) that 𝑊̇=𝑥̇_𝑡−𝑧𝑞̇−𝑧̇𝑞=𝐴(0)𝑊−2Re𝑞^∗(0)𝐹₀𝑞(𝜗),𝜗∈−𝜏₁,𝐴(0),0𝑊−2Re𝑞^∗(0)𝐹₀𝑞(0)+𝐹₀(,𝜗=0,=𝐴𝑊+𝐻𝑧,𝑧,𝜗), (47) where 𝐻(𝑧,𝑧,𝜗)=Re𝑔(𝑧,𝑧)𝑞(0()+𝐹0,𝑊𝑧,𝑧,𝜗)+2Re𝑧(𝑡)𝑞(𝜗)=𝐻₂₀(𝜗)𝑧²2+𝐻₁₁(𝜗)𝑧𝑧+𝐻₀₂(𝜗)𝑧²2+⋯. (48) Expanding the above series and comparing the coefficients, we get 𝐴(0)−2𝑖𝜔₀𝑊₂₀(𝜗)=−𝐻₂₀(𝜗),𝐴(0)𝑊₁₁(𝜗)=−𝐻₁₁(𝜗),𝐴(0)+2𝑖𝜔₀𝑊₀₂(𝜗)=−𝐻₀₂(𝜗). (49) Now 𝑞^∗(0)=𝐵 implies 𝑔(𝑧,𝑧)=𝐵𝑠𝑓^′′𝑘_∗2𝑥²_𝑡−𝜏1+𝐵𝑠𝑓^′′′𝑘_∗𝑥3!^{3𝑡−𝜏1}+⋯. (50) Noticing that 𝑥_𝑡(𝜗)=𝑊(𝑧,𝑧,𝜗)+𝑧𝑞(𝜗)+𝑧𝑞(𝜗), we derive 𝑥_𝑡−𝜏1=𝑊₂₀−𝜏₁𝑧²2+𝑊₁₁−𝜏₁𝑧𝑧+𝑊₀₂−𝜏₁𝑧²2+⋯+𝑒^{−𝑖𝜔0𝜏1}𝑧(𝑡)+𝑒^{𝑖𝜔0𝜏1}𝑧(𝑡). (51) Then substituting this into (50) and comparing the coefficients with (46), the following hold: 𝑔₂₀=𝐵𝑠𝑓^′′𝑘_∗𝑒^{−2𝑖𝜔0𝜏1},𝑔₁₁=𝐵𝑠𝑓^′′𝑘_∗,𝑔₀₂=𝐵𝑠𝑓^′′𝑘_∗𝑒^{2𝑖𝜔0𝜏1},𝑔₂₁=𝐵𝑠𝑓^′′𝑘_∗2𝑊₁₁−𝜏₁𝑒^{−𝑖𝜔0𝜏1}+𝑊₂₀−𝜏₁𝑒^{𝑖𝜔0𝜏1}+𝑠𝑓^′′′𝑘_∗𝑒^{−𝑖𝜔0𝜏1}. (52) In order to compute 𝑔₂₁, we need to know 𝑊₂₀(𝜗) and 𝑊₁₁(𝜗). For 𝜗∈[−𝜏₁,0], we have 𝐻(𝑧,𝑧,𝜃)=−2Re𝑞^∗(0)𝐹₀𝑞(𝜃)=−𝑔𝑞(𝜃)−𝑔𝑞(𝜃)=−𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2𝑞(𝜃)−+⋯𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+⋯𝑞(𝜃). (53) Comparing the coefficients with those in (48) yields 𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐻₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃). (54) From (49) 𝑊̇₂₀(𝜃)=2𝑖𝜔₀𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃)𝑊̇₁₁(𝜃)=𝑔₁₁𝑞(𝜃)+𝑔₁₁𝑞(𝜃). (55) Solving 𝑊₂₀(𝜃), we have 𝑊₂₀(𝜃)𝑔=−²⁰𝑖𝜔₀𝑒^{𝑖𝜔0𝜃}−𝑔₀₂3𝑖𝜔₀𝑒^{−𝑖𝜔0𝜃}+𝐸₁𝑒^{2𝑖𝜔0𝜃}, (56) and, similarly for 𝑊₁₁(𝜃), 𝑊₁₁(𝜃)=𝑔₁₁𝑖𝜔₀𝑒^{𝑖𝜔0𝜃}−𝑔₁₁𝑖𝜔₀𝑒^{−𝑖𝜔0𝜃}+𝐸₂. (57) Here 𝐸₁ and 𝐸₂ are constants to be determined by setting 𝜃=0 in 𝐻(𝑧,𝑧,𝜃). A direct computation shows 𝐸₁=𝑠𝑓^′′𝑘_∗𝑒^{−2𝑖𝜔0𝜏1}2𝑖𝜔₀+𝑠𝑓^′𝑘_∗−𝛿𝑒^{−2𝑖𝜔0𝜏1}+𝛿𝑒^{−2𝑖𝜔0𝜏02},𝐸₂=𝐵𝑠𝑓^′𝑘_∗. (58) Then all 𝑔_𝑖𝑗 have been obtained, and thus we can compute the quantities 𝑐₁(0)=𝑖2𝜔₀𝑔₁₁𝑔₂₀−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏⁰₂,𝛽₂=2Re𝑐₁(0,𝑇₂)𝑐=−Im₁(0)+𝜇₂Im𝜆^′𝜏⁰₂𝜔₀, (59) which determine the properties of bifurcating periodic solutions at the critical value 𝜏⁰₂. From the discussion above, we have the following result.

Theorem 9. Let 𝑘_∗ be the unique positive equilibrium of the model (2). Then one has the following. (1)𝜇₂ determines the direction of the Hopf bifurcation when 𝜏=𝜏⁰₂: if 𝜇₂>0 (resp., 𝜇₂<0), then the Hopf bifurcation is supercritical (resp., subcritical) and the bifurcating periodic solution exists for 𝜏>𝜏⁰₂ (resp., 𝜏<𝜏⁰₂) in a sufficiently small 𝜏⁰₂-neighbourhood.(2)𝛽₂ determines the stability of the bifurcating periodic solution: if 𝛽₂<0 (resp., 𝛽₂>0) the bifurcating periodic solution is locally asymptotically stable (resp., unstable).(3)𝑇₂ determines the period of the bifurcating periodic solution: if 𝑇₂>0 (resp., 𝑇₂<0) the period increases (resp., decreases).

5. Numerical Simulations

This section is concerned with some numerical simulations of the mathematical model (2) with the aim of exploring the analytical results. The model is characterized by four nonnegative parameters (𝑠, 𝛿, 𝜏₁, 𝜏₂) and the function 𝑓. In what follows we restrict our attention to the Cobb-Douglas function; this function reads 𝑓(𝑘(𝑡))=𝐴[𝑘(𝑡)]^𝛼, (60) with 𝛼∈(0,1) and 𝐴>0.

The first set of simulations refers to the following case: 𝑠=0.11,𝛿=0.8,𝛼=0.1,𝐴=1,𝜏₁=3,𝜏₂=2, (61) where 𝜏₁>𝜏₂. Figure 1(a) shows that the function 𝑘(𝑡) reaches the stationary state (the equilibrium 𝑘_∗). This equilibrium is stable; see Figure 1(b), where the evolution of 𝑘^′(𝑡) versus 𝑘(𝑡) is depicted for 𝑡∈[0,100].

(a)

(b)

Figure 1 

The time evolution of the function 𝑘(𝑡) for 𝑠=0.11, 𝛿=0.8, 𝛼=0.1, 𝐴=1, 𝜏1=3, and 𝜏2=2 (a). The 𝑘′(𝑡) versus 𝑘(𝑡), for 𝑡∈[0,100] (b).

The dynamics depicted by Figure 1 does not change from the qualitative viewpoint when 𝛼∈[0.02,1). Indeed when the parameter 𝛼 varies in the interval [00.2,1) the time necessary for reaching the stationary state increases but the behavior is that of Figure 1. When 𝛼<0.2 the time evolution of 𝑘(𝑡) is shown in Figure 2. In this simulation the time length has been increased in order to better visualize the evolution. As Figure 2(b) shows, the equilibrium is now instable.

(a)

(b)

Figure 2 

The time evolution of the function 𝑘(𝑡) for 𝑠=0.11, 𝛿=0.8, 𝛼=0.01, 𝐴=1, 𝜏1=3, and 𝜏2=2 (a). The 𝑘′(𝑡) versus 𝑘(𝑡), for 𝑡∈[0,300] (b).

The second set of simulations refers to the following case: 𝑠=0.11,𝛿=0.8,𝛼=0.1,𝐴=1,𝜏₁=1,𝜏₂=2. (62) These simulations take into account the case 𝜏₁<𝜏₂. As Figure 3 shows, oscillations occur for a long time (about 𝑡∈[0,1500]), with respect to the previous case, before reaching the equilibrium. In this case the equilibrium is stable but a very long time is necessary to reach it; see Figure 3(b). It is sufficient to increase 𝛼 for obtaining the stationary state rapidly. These simulations suggest that when the value of 𝛼 increases then the time necessary to reach the stationary state decreases. Moreover, if the difference 𝜏₂−𝜏₁>0 increases, then the stability of the equilibrium is lost for all 𝛼∈(0.1).

(a)

(b)

Figure 3 

The time evolution of the function 𝑘(𝑡) for 𝑠=0.11, 𝛿=0.8, 𝛼=0.1, 𝐴=1, 𝜏1=1, and 𝜏2=2 (a). The 𝑘′(𝑡) versus 𝑘(𝑡), for 𝑡∈[0,1300] (b).

Finally, we would show some numerical simulations related to the evolution of 𝑘^′(𝑡) versus 𝑘(𝑡). Figure 4 shows the instability of the equilibrium for 𝑠=0.41, 𝛿=0.35, 𝛼=0.8, 𝐴=1, 𝜏₁=10, 𝜏₂=2, see Figure 4(a); it is worth stressing that the equilibrium is reached for 𝑡≥3500. In this case, when we increase the magnitude of 𝑠 we are able to reach the equilibrium more rapidly, see Figure 4(b), which is obtained for 𝑠=0.91 and 𝑡∈[0,3500]. The equilibrium is rapidly reached if we also decrease the magnitude of 𝛼.

(a)

(b)

Figure 4 

The time evolution of 𝑘′(𝑡) versus 𝑘(𝑡) for 𝑠=0.41, 𝛿=0.35, 𝛼=0.8, 𝐴=1, 𝜏1=10, 𝜏2=2, and 𝑡∈[0,3500] (a). The time evolution of 𝑘′(𝑡) versus 𝑘(𝑡) for 𝑠=0.91, 𝛿=0.35, 𝛼=0.8, 𝐴=1, 𝜏1=10, 𝜏2=2, and 𝑡∈[0,3500] (b).

6. Conclusions and Research Perspectives

In the present paper, a generalization of the Solow model by inserting two time delays has been considered. The delays, respectively, represent the time employed in order that the capital is used for production and the necessary time so that the capital is depreciated. Specifically, an asymptotic analysis has been performed referring to the stability analysis of the steady state and the conditions under which a Hopf bifurcation appears.

According to the analysis developed in this paper, the stability of the positive equilibrium changes as the time delays vary. Indeed if 𝜏₂=0, the positive equilibrium is always locally asymptotically stable; if 𝜏₁=0 the positive equilibrium can be locally asymptotically stable or unstable and a Hopf bifurcation occurs; the dynamics is more complicated when the two time delays are both different from zero (see Theorem 8 where as the reader can see the investigation of stability switches becomes quite complicated). This shows that the time delays play an important role in the dynamics of the model. Then, based on the analysis of the existence of the Hopf bifurcation, by using the center manifold theory and the normal form method, an explicit algorithm for determining the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions has been derived. This means that one can obtain the important quantities which determine the properties of bifurcating periodic solutions at the critical value; see Theorem 9. According to our results, we can say that the model with two independent time delays has much more complicated dynamics than the model with only one time delay. That is why it seems to be more realistic.

The introduction of time delays can be also performed in the mathematical model developed in [36] for the mammary carcinoma. Indeed the stability analysis developed in the present paper can help to reach more results in the cancer-immune system competition. From a biological point of view, the Hopf bifurcation means that for small values of parameters the nontrivial stationary solution to the model in [36] is stable, and we do not observe radical changes in the competition. Otherwise, nontrivial stationary solution can oscillate and the amplitude of the oscillations about the stationary solution remains constant. This case simply corresponds to the situation when the competition oscillates in time.

The Hopf bifurcation analysis developed in this paper must be revised if the mathematical models are not based on ordinary differential equations. Recently an increasing number of partial differential equation models for tumor growth or therapy have been developed; see the references section of paper [12] and the references cited in the recent review paper [37].

Moreover thermostated integrodifferential equations have been proposed in papers [38–43] for the modeling of biological systems, vehicular traffic, crowd and swarm dynamics, and economic systems subjected to external force fields. The introduction of the Gaussian isokinetic thermostat ensures the reaching of stationary states whose existence has been proved in [44]. The introduction of multiple time delays in thermostated equations, their stability, and bifurcation analysis is a future research perspective.

It is worth stressing that also the Boltzmann equation with the one-dimensional Bhatnagar-Gross-Krook relaxation type operator [45] and the Kac equation have been coupled with a Gaussian isokinetic thermostat; the existence of stationary solutions is ensured also within these frameworks; see papers [46–48].

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The first author acknowledges the support by the FIRB project RBID08PP3J-Metodi matematici e relativi strumenti per la modellizzazione e la simulazione della formazione di tumori, competizione con il sistema immunitario, e conseguenti suggerimenti terapeutici. The authors acknowledge the financial support of MEDAlics-Research Center for Mediterranean Relations.