Direct and inverse relationships between Riccati systems coupled with multiplicative terms

ABSTRACT An analytical and computational framework for the derivation of solitary solutions to biological systems describing the cooperation and competition of species and expressed by the system of Riccati equations coupled with multiplicative terms is presented in this paper. It is demonstrated that relationships between these solitary solutions can be either direct or inverse. Thus, an infinitesimal perturbation of one population would lead to an infinitesimal change in the other population – if only both solitary solutions are coupled with the direct relationship. But, in general, that is not true if solitary solutions are coupled with the inverse relationship – an infinitesimal perturbation of one population may result into a non-infinitesimal change in the other population. Necessary and sufficient conditions for the existence of solitary solutions are derived in the space of the system's parameters and initial conditions.


Introduction
Cooperation and competition of species in biological systems described by nonlinear differential equations has attracted the attention of researchers since the introduction of the paradigmatic Lotka Volterra model at the beginning of the twentieth century [5]. The interest to transient nonlinear effects in competing populations has been also amplified by the development of mathematical and computational techniques dedicated for the analysis of solitary waves in the middle of the twentieth century. Many nonlinear models, including the Korteweg-de-Vries equation, the nonlinear Schrödinger equation, the sine-Gordon equation, the Lax equation, have been studied in the context of the existence and the qualitative description of solitary solutions in these systems [4,26]. The main objective of this paper is to demonstrate that the existence of solitary solutions is also possible in a system of coupled Riccati equations with the multiplicative term.
Coupling of two ordinary differential equations by multiplicative terms is a well-known technique in mathematical models describing dynamics of biological systems. The multiplicative coupling term is used to describe the air-borne disease transmission rate (the mass action term) in the model of a microparasite in a host population with the Allee effect [3,12,16]. Two nonlinear partial differential equations coupled with multiplicative terms are also used to describe the spatio-temporal dynamics of a predator-prey system where the prey per capita growth rate is subject to the Allee effect [14]. The assumption that per capita mortality rate of prey and predator is equal and the rate of biomass production must be consistent with the rate of biomass assimilation helps to simplify the model to a diffusive predator-prey system coupled with multiplicative terms [13].
The travelling-wave reduction transforms this model into a coupled system of ordinary differential equations [8]. The (G /G) expansion method is used to construct analytical solitary solutions to the coupled equations in [8].
A generalization of the (G /G) method called the simplest equation method was presented in [9] and further extended in [1,10]. The simplest equation method is powerful tool that is used to construct exact solutions (including solitary solutions) to nonlinear differential equations that is based on the expression of the solution to the given nonlinear differential equation using the simplest nonlinear differential equation [9]. The exact solitary wave and periodic wave solutions to the Kuramoto-Sivashinsky equation are constructed in [9]. The simplest equation method has been applied to study pattern formation processes on the semiconductor surfaces under ion beam bombardment [11]. The exact solutions to the Painlevé equations are derived using the method in [1].
The modified simplest equation method, presented by Vitanov in [18,19], is a powerful modification of the simplest equation method that can be used to construct exact solutions to many nonlinear differential equations. The modified simplest equation method is applied to obtain exact solutions to a family of nonlinear partial differential equations (PDEs) containing the Kuramoto-Sivashinsky equation, reaction-diffusion equation with density-dependent diffusion, and the reaction-telegraph equation in [24]. The travellingwave solutions to the generalized Rayleigh and Swift-Hohenberg equations are obtained using the modified simplest equation method in [18]. The Bernoulli and Riccati simplest equations are used to construct travelling-wave solutions for a class of nonlinear PDEs with polynomial nonlinearity in [17]. Solitary wave solutions to differential equations containing both odd and even grade monomials are investigated in [21].
Solitary wave solutions are an often observed physical phenomenon. It has recently been shown that the dynamical system describing a two-layer immiscible fluid flow subject to horizontal harmonic vibrations has solitary solutions [2,6]. Solitary solutions of the Zakharov-Kuznetsov equation used to describe the behaviour of weakly nonlinear ionacoustic waves in plasma comprising cold ions and hot isothermal electrons have been recently considered in [27]. They are also observed in the modelling of dispersive shock waves [7].
The solitary solutions to differential equations modelling population dynamics can be constructed using the modified simplest equation method. It is demonstrated in [23] that the exact nonlinear kink and solitary wave solutions to a model system of partial differential equations for description of the spatio-temporal dynamics of interacting populations can be derived. Coupled kink waves in a system of two interacting populations in which the reproduction and intensity of interaction depend on their spatial density are constructed in [22]. Exact travelling-wave solutions to the reaction-diffusion and reaction-telegraph equations that are used in ecology and population dynamics are obtained in [20]. The modified method of the simplest equation is used to solve issues from the dynamics of interacting populations in [25].
As mentioned previously, the main objective of this paper is to derive necessary and sufficient conditions for the existence of solitary solutions in the system of Riccati equations coupled with the multiplicative terms: where A 2 , B 2 = 0. The insight into the structure of solitary solutions to (1) and the relationship between these solutions would be valuable for understanding nonlinear dynamical processes in biological systems coupled with multiplicative terms and extend the existing results on the dynamics of coupled population systems obtained by the modified simplest equation method [20,22,23,25]. This paper is organized as follows. Riccati equation and its generalizations are discussed in Section 2; systems of Riccati equations and solutions to these systems are derived in Section 3; direct and inverse relationships between solitary solutions are analysed in Sections 4 and 5; concluding remarks are given in the final section.

The Riccati equation and its general solution
The Riccati equation with constant coefficients reads [15]: where A 0 , A 1 , A 2 ∈ C; A 2 = 0. Riccati equation (2) can be rearranged into the canonical form: It is well known [15] that the general solution y = y(x, c, s) to (3) reads: where x is the independent variable; s is the initial value of y at x = c: y(c, c, s) = s; σ = A 2 and y 1 = y 2 are two different roots of the polynomial A 2 y 2 + A 1 y + A 0 . Note that y(x, c, y k ) = y k ; k = 1;2.

Remark 1:
Let the solution of Riccati differential equation (3) be given in the form: where a, b,â,b, α, β ∈ C are the parameters of Riccati equation (3). Then the following relationships hold true:

Systems of Riccati differential equations
In this section we discuss the systems of two uncoupled and coupled Riccati differential equations and the properties of their solutions. Let us consider an uncoupled system of two Riccati equations in the canonical form: We will assume that parameters of this system y 1 , Then the solution of the system (5) y = y(x, c, s); z = z(x, c, t) with Cauchy boundary conditions y(c, c, s) = s; z(c, c, t) = t, where c ∈ C, reads: Note that there is one-to-one correspondence between (5) and (6) in the sense that from (5) we construct (6) and vice versa, knowing the expressions (6) we can find all parameters of the system (5).
We give two definitions.

Definition 3.1:
We say that there is direct relationship between functions (6) if there exists such constants α ∈ C \ {0} and β ∈ C such that for all x and c equality holds: Definition 3.2: We say that there is inverse relationship between functions (6) if there exists such constant γ ∈ C \ {0} such that for all x and c equality holds: Let us denote
Proof: Relationships (11) we get directly from (9). The correctness of (7) we obtain by using transformations: Theorem 3.4: Let parameters of the system (5) be and Then the solution (6) of the system (5) satisfies (8), that is, there exists an inverse relationship between functions y and z when ts = γ .

Corollary 3.5:
Let parameters of the system (5) satisfy (12) and (13) and σ = 0. Then the following equality holds true: Moreover, if then (7) holds true for all values of x and s when t = αs + β.
In other words, the inverse relationship (8) and the direct relationship (7) hold true simultaneously if and only if (12), (13) and Cauchy initial conditions hold true.

Remark 2:
Note that in general case the direct relationship doesn't imply the inverse relationship. Furthermore, there exist other relationships between solutions (6) of the system (5). We will not discuss other types of relationships since they don't describe phenomenons of coexistence of several dynamical systems, which can be expressed by solitary solutions (6).
Let us consider the coupled system of two Riccati equations: dy dx = a 0 + a 1 y + a 2 y 2 + λyz, Assume that a 0 , a 1 , a 2 , λ, b 0 , b 1 , b 2 , μ ∈ C and a 2 , b 2 , λ, μ = 0. Suppose that functions y = y(x, c, s) and z = z(x, c, t) with Cauchy boundary conditions y(c, c, s) = s and z(c, c, t) = t is the solution of the system (17).
The following theorem holds true.  (17) are expressed in the form (6). Then those functions at the same time are the solution of an uncoupled system of Riccati equations, which can be written in this general form: where A 0 , A 1 , A 2 , B 0 , B 1 , B 2 ∈ C and A 2 , B 2 = 0. Then from (17) and (18) for all x, s, t the following equalities hold true for functions y = y(x, c, s), z = z(x, c, t): A 0 + A 1 y + A 2 y 2 = a 0 + a 1 y + a 2 y 2 + λyz, or Since both (19) and (20) must be true at the same time then we have two exclusive cases: Case 1: that is, functions y = y(x, c, s), z = z(x, c, t) have the direct relationship (7).
We note that other relationships between functions y = y(x, c, s), z = z(x, c, t), which satisfy identities (19) and (20) do not exist. This proves the necessity part of the theorem.
Sufficiency. Assume that functions y = y(x, c, s), z = z(x, c, t) have the direct relationship (7). Then we can write the system of two Riccati equations in canonical form: from which we get the coupled system of two Riccati equations (17) in general form. Similarly, let functions y = y(x, c, s), z = z(x, c, t) have the inverse relationship (8). Then we have such system of two Riccati equations in canonical form: from which we also get the coupled system of two Riccati equations (17) in general form. This proves the sufficiency part.

Direct relationship
Let (7) hold true. Then the second differential equation of the system (17) yields: Elementary transformations reduce (25) into the form: Then, necessary and sufficient conditions for (26) to coincide with the first differential equation of the system (17) read: Thus, the relationship (7) holds true if and only if: when = 0 then it would be either a degenerate case discussed below or it would be a nonexistent case of solitary solutions. For this reason here and below we assume that those expressions are not zeros and won't discuss this anymore.
The solution to (34) cannot be described by (35) if the relationship between initial conditions (36) does not hold. This fact is illustrated in Figure 2. Thick solid lines represent the solution at c = 0; s = 1.9 and t = 1 (the constraint (36) does not hold true). It can be observed that lim x→∞ y(x,0,1.9) z(x,0,1) = −1 1 , but the transient process is much more complex than described by (35) -note the local minimum of z(x, 0, 1) at around x = 0.2. Dashed lines represent the solitary solution at c = 0; s = 1.9 and t = 2 14 15 in Figure 2 (the constraint (36) does hold then).
The validity of results can be double-checked by the following computational experiment. We will construct approximate numerical solutions to (34) at c = 0; s k = −1 + 0.015k; t l = 1 + 0.01l; k, l = 0, 1, . . . , 200 using constant step forward marching techniques. Let us denote approximate solutionỹ(jh);z(jh); j = 0, 1, 2, . . . where h is the step size. The exact analytical solution (35) is defined on the parameter line (36) in the phase space of initial conditions. But we release the constraint (36) and assume that the solution (35) is valid throughout the plane of initial conditions. Using RK2 method we travel 100 steps (h = 0.1) from starting point x = c with the preselected pair of initial values y = s  and z = t and compute differences between the approximate numerical solution and the exact solution (35). Adding absolute differences for 100 steps produces an error estimate: The distribution of ε(s, t) is shown in Figure 3; numerical values of ε(s, t) higher than 1 are truncated to 1 in order to make the figure more comprehensive. It is clear that errors are almost equal to zero on the line 3t = 2s+5.

Inverse relationship
Suppose that the solution to (17) satisfies condition (8). Then, the second differential equation of the system (17) yields: Necessary and sufficient conditions for (38) to coincide with the first differential equation of the system (17) read: Then, necessary and sufficient conditions for the condition (8) to hold true reads: The parameter γ can be determined from the equality: Then, (17) can be rearranged into the inverse canonical form: Since the equality yz = γ holds for all x, c, s and t, the first differential equation of the system (40) and the equality (8) yield: which can be rearranged into the following form: The second equation of the system (40) is equivalent to (42). Therefore, Finally, the algorithm for the determination of unknown parameters σ 1 , y 1 and y 2 reads: Then, the solution to (40) reads: if and only if t = γ /s. Note that this is the same solution as (6) except it has the constraint t = γ /s.
System's parameters are: a 0 = 2; a 1 = −6; a 2 = 2; λ = 1; b 0 = −2; b 1 = 6; b 2 = −2; μ = −1. Conditions (39) hold true and thus this system can be rearranged into the inverse form. Thus, γ = 2; y 1 = 1; y 2 = 2; z 1 = 2; z 2 = 1, σ 1 = 2, σ 2 = −2. Then the inverse canonical form of (44) reads: The solution to (45) reads if and only if Note that z(x, c, 2/s) = 2/y(x, c, s). The inverse relationship between y and z can be seen in Figure 4. The validity of results can be double-checked by the following computational experiment. We will construct approximate numerical solutions to (44) at c = 0; s k = 1 + 0.005k; t l = 1 + 0.005l; k, l = 0, 1, . . . , 200 using constant step forward marching techniques. Let us denote approximate solutionỹ(jh);z(jh); j = 0, 1, 2, . . . where h is the step size. The exact analytical solution (46) is defined on the parameter line (47). But we release the constraint (47) and assume that the solution (46) is valid throughout the plane of initial conditions. Using RK2 method we travel 100 steps (h = 0.1) from starting point x = c with the preselected pair of initial values y = s and z = t and compute differences between the approximate numerical solution and the exact solution (46). Adding absolute differences for 100 steps produces an error estimate ε(s, t) (37). The distribution of ε(s, t) is shown in Figure 5; numerical values of ε(s, t) higher than 5 are truncated to 5 in order to make the figure more comprehensive. It is clear that errors are almost equal to zero on the curve t = 2/s.

Concluding remarks
The main objective of this paper is to present an analytical and computational framework for the derivation of solitary solutions to the system of Riccati equations coupled with multiplicative terms. It is shown that solitary solutions can exist in the space of system's parameters and initial conditions. Moreover, these solitary solutions can be oriented in a direct or an inverse relationship -necessary and sufficient conditions for this orientation are derived in the explicit form (Equalities (28) and (39)). The existence of direct and inverse relationships between solitary solutions implies other important properties of the coupled system of Riccati equations. As mentioned inSection 1, the coupled system of Riccati equations with multiplicative terms can be used for the description of the dynamics of biological systems where y(x) and z(x) describe two competing populations. Thus, an infinitesimal perturbation of one population would lead to an infinitesimal change in the other population -if only both solitary solutions are coupled with the direct relationship. But, in general, that is not true if solitary solutions are coupled with the inverse relationship -an infinitesimal perturbation of one population may result into a non-infinitesimal change in the other population. And though we do not try to speculate regarding the biological interpretation of such effects (such interpretations would be a definite objective of the future work), such consequences are interesting from purely theoretical point of view.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding information
Financial support from the Lithuanian Science Council under Project No. MIP078/15 is acknowledged.