Constructing Soliton and Kink Solutions of PDE Models in Transport and Biology

We present a review of our recent works directed towards discovery of a periodic, kink-like and soliton-like travelling wave solutions within the models of transport phenomena and the mathematical biology. Analytical description of these wave patterns is carried out by means of our modification of the direct algebraic balance method. In the case when the analytical description fails, we propose to approximate invariant travelling wave solutions by means of an infinite series of exponential functions. The effectiveness of the method of approximation is demonstrated on a hyperbolic modification of Burgers equation.


Introduction
It is well known that methods of obtaining general solutions for the majority of nonlinear evolutionary PDEs actually do not exist. The already known methods of integration, such as the inverse scattering technique, are applied merely to so called completely integrable equations [1]. These equations form a very important but relatively small family, which does not include e.g. nonlinear dissipative models, describing the spatio-temporal patterns formation and evolution [2]. One of the few alternatives to numerical methods of the investigations of those models that are not completely integrable are the symmetry methods [3,4] which are effective when the system under consideration admits a non-trivial group of transformations. Yet, these methods cannot be treated as universal tools for obtaining analytical solutions. In fact, the self-similar reduction, when applied to a nonlinear PDE, gives another nonlinear differential equation with a fewer number of independent variables and it is difficult to expect that the procedure will lead to quadratures in every particular case. Obviously the chances to obtain this way the analytical description of solutions with the required properties and asymptotic behavior (such as e.g. periodic, kink-like or soliton-like travelling wave solutions) are, generally speaking, very small. In the case when the self-similar reduction leads to a system of autonomous ODEs the answer to the question on whether the solutions with the given properties do exist can be obtained by means of qualitative theory methods [5,6], which allow for studying of whole families of solutions with the given symmetry. The qualitative theory methods are very effective if one wants to know whether the family of self-similar solutions possesses solutions with the given properties. Yet in many cases mere information about their existence is not sufficient. Very often one needs to have an analytical expression for the given solution and this is the case when it is necessary to investigate the influence of small addends upon the stability and the evolution properties of the invariant solutions.
Compared to classical symmetry methods, more powerful tools for finding out solutions with the given properties, can be provided by so-called direct algebraic balance method, based upon representing solutions in the form of rational algebraic combinations of properly chosen functions, closed with respect to the algebraic operations and differentiation [7,8,9,10]. However, within this method one must solve, as a rule, very large systems of nonlinear algebraic equations. Regardless of this purely technical reason, the ansatz-based method cannot be treated as the remedy, because there are known evolution PDEs which possess analytical solutions that cannot be obtained, at least, within a particular set of elementary functions (hyperbolic, trigonometric, polynomial, etc.). And it is rather doubtful that a fully universal ansatz for finding out the analytical solutions to nonlinear PDEs can ever be found. Nevertheless, it is possible to answer the question of whether the solution of this or that sort does exist within the given set of invariant solutions, whenever the procedure of a self-similar reduction leads to the system of autonomous ODEs. And in the situation when the direct algebraic balance method, based on the representing solution we are looking for in the form of finite combination of the given functions fails, we propose to express it by means of the infinite series.
The aim of this study is to review our recent work directed towards the description of the already mentioned wave patterns, i.e. the kink-like and soliton-like travelling wave solutions, in various non-local models, simulating the transport phenomena. In our earlier work we investigated merely the set of travelling wave (TW) solutions, being invariant with respect to the group of the spatio-temporal translations, nevertheless, systems obtained via the symmetry reduction prove to be sufficiently difficult to analyze and solve exactly. Therefore we use a complex approach including a qualitative analysis, our own interpretation of the ansatz-based method and approximation of the invariant solutions by means of infinite series of exponential functions.
The structure of the article is as follows. In Section 2 we present a certain modification of the direct algebraic method, put forward in papers [7,8,9]. Employing this modification we obtain the kink-like and soliton-like TW solutions for the hyperbolic modification of Burgers equation, telegraph equation, nonlinear d'Alembert equation and the somitogenesis model. In Section 3 we present the results of a qualitative investigations of a set of TW solutions for the hyperbolic modification of Burgers equation. It demonstrates the existence of the TW solutionlike solutions, yet these solutions cannot be described by the finite algebraic combination of elementary functions. Therefore we propose to approximate the invariant solitary wave solutions by an infinite series of the exponential functions.
2 Exact soliton-like and kink-like solutions for transport equations

Modification of the direct algebraic method and its applications
The essence of the unified direct algebraic method in E. Fan's paper [7], is based on the observation that particular solutions of any system of PDE's which do not depend on (t, x) coordinates in explicit form, can be presented as a polynomial expansion where a µ are unknown parameters, while function φ(ξ) satisfies the equatioṅ Note that in order to obtain non-trivial solutions, one should "balance" the numbers r and s (for details see [7]). Depending on the conditions posed on the parameters c ν , solutions to the equation (3) are expressed by elliptic Jacobi functions, and in particular cases by hyperbolic or trigonometric functions [7]. The properties of these functions are inherited by the solution of the initial system which can be presented in the form (2). In fact this methodology is constructive and algorithmic if the functions H ν (·, ·, . . . ) arising in (1) are algebraic or rational ones. With this assumption, for any set of functions u i m i=1 that can be represented as (2), the LHS of the equations (1) can be presented in the form of finite linear combinations of functions φ k (ξ) and φ l (ξ)φ(ξ), k, l = 0, 1, 2, . . . : When equating to zero all the coefficients of these decompositions, one obtains a system of nonlinear algebraic equations, that determine the particular solutions of the initial system of PDE's (for more details see [7]). It must be stressed that the method does not guarantee obtaining of the general solution, even within the set (2), for its effectiveness depends on one's ability to solve a system of nonlinear algebraic equations. On the other hand, it enables calculation of particular solutions possessing some predictable properties. An enforced version of Fan's method was put forward recently in [8,9]. It was used for determining of exact solutions of the equation More precisely, the following ansatz was proposed: where z(ξ) = k µ=0 a µ φ µ (ξ), and φ(ξ) is a function satisfying equation (3). Due to this combination, the multi-parameter families of the exact solutions were obtained in the situation when the pure form of Fan's methodology does not work [8].
On analyzing different versions of the ansatz-based method, one can conclude that effectiveness of their employment is based on the mere fact that the family of functions φ µ (ξ) and φ νφ (ξ) is closed with respect to the algebraic operations and differentiation. In view of this, a quite natural generalization of the already mentioned ansatzes can be as follows: where the function φ(ξ) still satisfies the equation (3), but, in contrast to (4), dependence between the functions f and g is not assumed from the very beginning. Effectiveness of the ansatz (5) is demonstrated by applying it to some transport equations given below.

Soliton-like and kink-like TW solutions for the generalized transport equations
Let us consider the following equations: where τ , A, B, κ are non-negative constants. For A = 0 equation (6) coincides with the nonlinear telegraph equation; for A = 0 it coincides with the hyperbolic generalization of Burgers equation, while for A = B = 0 -with the nonlinear d'Alembert equation. The hyperbolic modifications of nonlinear transport equations arise naturally when the effects of temporal non-locality are taken into account [11].
Assuming that the solitons and kinks can be expressed by powers of function sech(ξ), which is the particular solution of equation (3) and function sinh(ξ) which appears in the odd derivatives of the function sech(ξ), we use the following ansatz: or, equivalently, Inserting the ansatz (7) ( or (8)) into (6) and executing all the necessary operations, we obtain an algebraic equation containing, respectively, functions sech µ (αξ), sech ν (αξ) sinh(αξ) or exp [µαξ]. Deeming them functionally independent and equating to zero the corresponding coefficients, we go to the system of algebraic equations. We do not demonstrate the details of these calculations since they are simple but cumbersome. The calculations were performed through the aid of a software application "Mathematica". The results obtained are presented below.
I. For arbitrary A, B and f (u satisfies equation (6) if the following conditions hold: Here and henceforth we use the notation (9) defines a kink-like regime when b 0 b 2 > 0 and a 0 /b 0 = a 2 /b 2 . Using the conditions (10), we can express the unknown parameters from formula (9) by the parameters characterizing equation (6), yet, in the general case it is too cumbersome. It is much easier to From these conditions we obtain the solution: For κ = 1, τ = 0, A = 0, B = 1 function (9) coincides with the solution obtained in [8].
Another example of the kink-like solution defined by the formulae (9), (10) is as follows: .

II. For
satisfies the equation (6) providing that the following conditions hold: This solution defines the solitary wave regime satisfies (6), when λ 1 , λ 3 are positive and the parameters are as follows: This solution is always singular, because for arbitrary values of the parameters the expression in the denominator of the formula (11) vanishes for some ξ ∈ R.
IV. Now let us consider the case satisfies equation (6) when the following conditions hold: For a 0 = 0, b 0 = 0, |a 1 | + |b 1 | = 0 equation (12) defines the soliton-like solution. One of the parameters contained in (12) can be chosen arbitrarily, while the rest can be expressed, using (13), as the functions of the parameters, defining equation (6). We omit doing this in the general case, but present one particular example. Thus, for b 1 = 0, b 0 = 1 and arbitrary α we have the solution and the additional condition defines a solution of (6) if the following conditions hold: VId. For f (u) = λ 1 u(t, x) + λ 3/2 u(t, x) αh defines a soliton-like solution of the equation (6).

A model of mathematical biology
Now we consider the following system [12]: where f (u, z), g(u, z) are given polynomials, and A is a positive constant. Due to the fact that the system is autonomous we apply variables in a travelling wave form: where ξ = x + µt, and we reduce system (16) to the following system of ODE: To simplify the computations, we restrict the family of admissible polynomials in RHS as follows: In order to obtain soliton-and kink-like solutions we proposed the following ansatz that shows the proper asymptotic behavior.
3 Approximated soliton-like solutions of the generalized Burgers equation

Existence of the soliton-like regimes within the set of TW solutions
We showed in the previous section that the ansatz (8) is quite effective in obtaining travelling wave solutions with given properties and asymptotics. Yet recently we realized that it is quite impossible to describe soliton-like solutions of the following generalization of the Burgers equation (GBE from now on): using the ansatz (8), which is more general than the generalized Hirota-like ansatzes employed in [8,9]. The first thing we are going to do in the situation when the ansatz-based methods fail, is to perform a qualitative study of the ODE system describing the TW solutions of the GBE and answer the question whether (and when) soliton-like solutions exist. Below we analyze in detail the whole set of TW solutions of equation (17). The results of this analysis are used in the following section in constructing the approximated soliton-like solutions.
Inserting the ansatz u = U (ξ) = U (x + µt) into the equations (17) we obtain the second order ODE leads us to the equation where A = 1/ hγ 2 andμ = µ/z 0 (we omit bars over the variables in the following formulae). Let us rewrite (19) as the following system of the first order ODEs: System (20) has three stationary points: B 0 (0, 0), B 1 = (−1, 0) and B 2 = (1, 0). Analysis shows [13] that the critical point B 0 is always a saddle. Depending on the values of the parameters, two other critical points are either nodes or foci. As was declared earlier, we are looking for the homoclinic trajectory and since the dynamical system (20) is of dissipative type, this trajectory can exist only for special values of the parameters. In our attempts to capture the homoclinic loop we begin with studying the condition of the limit cycle creation in proximity of the critical point B 2 . If there exists a limit cycle lying in the right half plane and surrounding the critical point B 2 , then its evolution in the vicinity of the unmovable saddle point B 0 , caused by the change of the driven parameter µ, most probably will lead to the homoclinic bifurcation.
It can be shown by inspection of the Andronov-Hopf theorem conditions and their fulfillment [14,13] that in vicinity of the critical value of the parameter µ = −1 an unstable limit cycle is created with zero radius. A further evolution of system (20) was undertaken by means of the numerical simulation in which we put A = 1. It shows that the radius of the periodic trajectory grows with the parameter µ, until the homoclinic bifurcation takes place at µ ≈ −0.836. Note that the homoclinic solution intersects the horizontal axis at the point U = x * ≈ 1.426095. On further growth of the parameter µ the homoclinic loop is destroyed and one of the phase trajectories starting at the unstable saddle B 1 reaches the critical point B 2 . Patterns of the phase trajectories of system (20) obtained by means of numerical simulations are shown in Fig. 2.
Thus, the homoclinic loop exists in system (20) merely at the specific values of the parameters. The construction of the approximated soliton-like solution, undertaken below is based in essential way upon the preliminary information obtained in this section.

Approximation of the soliton-like TW solution to the GBE
The results of qualitative analysis tell us that equation (17) possesses the SW solution for the specified values of the parameters. Let us assume that it can be described by means of the formula (8). Using the fact that the SW solution is represented by the homoclinic trajectory of system (19), that is bi-asymptotic to the origin, we can use for its description the following representation Proof . Since there is a one-to-one correspondence between the solutions of equations (18) and (19), we consider the latter one, trying to represent its homoclinic solution by means of the formula Let us denote the numerator and denominator of the above expression by F (T ), G(T ) respectively. Inserting the ansatz (22) into the equation (19) and multiplying the resulting equation by G 3 , we get the expression Now, gathering the coefficients standing at functions exp (αT ), exp [(3m + 2)αT ] and assuming that the coefficients a 1 , a m , b m+1 are nonzero, we obtain the following system of algebraic equations: that leads to a contradiction when A is nonzero.
Of course, the formula (22) does not give the most general expression for a solution, having the proper asymptotics at T → ±∞ and belonging to the family (8). For example, we can take the formula (23) instead of (22) for any natural numbers n 1 , n 2 , n 3 and m. But our experiments with this formula in the cases when the numerator and denominator contained, respectively, two and three terms, gave a result which is identical to that stated in Lemma 1. We do not mention it here because it is too cumbersome. Concerning an expression containing more terms the result is unclear as yet, but it is rather doubtful that we would obtain a positive answer by incorporating extra terms, for a number of extra algebraic equations accompanying the "inflation" of the formula (23) exceeds in an essential way a number of extra coefficients a k , b r to be determined.
So putting aside attempts to discover the analytical expression for soliton-like solutions of the GBE, we are searching for approximated solutions to the equation (19), represented by the series (24), with the additional condition where x * is a point of intersection of the homoclinic trajectory with the horizontal axis, which, generally speaking, can be obtained from the numerical experiment. With this conditions the two branches of the homoclinic trajectory are "tied up" in the point (x * , 0). There are two ways of determining the unknown constants a k and α. Firstly, it is possible to insert the series into (17) and equate to zero the coefficients standing at different powers of the function exp [αT ]. This way we obtain the recurrence a k+1 = ψ k+1 (a 1 , a 2 , . . . , a k ). Another possibility is to use the initial conditions (25) and their differential consequences. As in the previous case, this approach is quite algorithmic for any second order ODE that is linear with respect to the highest derivatives. The approach is based on the following statement, which can be easily proved by means of mathematical induction.

Lemma 2. Consider the Cauchy problem
Suppose we are looking for the approximate solution to (27), (28) in the form Then the coefficients {a k } N k=1 can be obtained from the conditions (28) and the N − 2 equations where D ξ denotes the total derivative. The functions Φ j (U (0), . . . , U (j+1) (0)) are recurrently expressed merely by the initial conditions x 0 , x 1 and do not depend upon a k .
So according to Lemma 2, using Cauchy data (28) and their differential consequences we get the linear system The matrixM N standing at the RHS of equation (30) is of the Vandermonde type [15]. Its determinant is given by the formula so it is always nonsingular and the system has a unique solution for each N .
Remark. The above lemma says nothing about the convergence of the series (29) to the real solution of the initial value problem (27), (28) as N → ∞. This problem needs a special investigation.
Our analysis suggests that the upper and lower branches of the homoclinic trajectory are so strongly asymmetric that the convergency of the approximated solution based on the first method of determining the coefficients of series expansions (24) takes place merely for the lower one, whereas the upper branch should be approximated with the second method only.
Thus, we approximate the lower branch of homoclinic solution by the series We omit the presentation of this recurrence since it is irregular and cannot be expressed it in the analytical form. Therefore all the calculations were performed with the help of the "Mathematica" package. They were based in an essential way upon the assumption that the phase trajectory approaches the point U = x * = 1.426095 of the horizontal axis as T → 0.
Approximate solutions corresponding to different values of the number N are shown in Fig. 3. Unfortunately we fail to employ a similar series decomposition and the method of finding the unknown coefficients for the upper branch, since the series do not converge to the numerical solution in proximity of the far end of the homoclinic trajectory, as it is seen in Fig. 4. In fact, the series is convergent yet the limiting function does not coincide with the solution, obtained by means of the highly precise numerical scheme. So in order to determine the series coefficients,    we turn to the second scheme. Keeping in the series (26) N elements and substituting the truncated series into equation (19), we obtain from (25) two conditions for the coefficients a k : The first differential consequence of these conditions gives us equation (19) itself. Taking it into the point T = 0 and employing equations (31), we obtain: In order to define all the coefficients of the series truncated on N -th element we need N − 3 extra conditions. We can obtain these missing conditions from further differential consequences, as it is stated in Lemma 2. Coefficients {a k } N k=1 satisfy the linear system (30) that delivers a unique solution, providing the RHS and the parameter α are determined. The remaining N − 3 components of the vector, positioned at the RHS of system (30), were obtained by means of some recurrent procedure, written in "Mathematica" package. In order to obtain the parameter α, we insert the solution (26) into (19) and equate to zero the coefficient α 2 +Aµα−1, at the lowest power of the exponential function. For A = 1 and µ = −0.836 this gives us the value α = 1.5018.
The above procedure was used to obtain the upper branch of the homoclinic curve of system (20) (or what is the same, the left "wing" of the SW solution of the equation (19)). The results obtained are shown in Fig. 5. The upper and lower branches sewed together are shown in Fig. 6, from which we can see that suggested approach quite well approximates the homoclinic solution of equation (17).

Concluding remarks
We have reviewed in this study the methods and tools that allow one to describe the solutions of the nonlinear PDEs with the required features and asymptotic behavior. The exact solutions presented here are obtained with the help of the ansatz (8) 1 , which proves to be effective when the methods described in [7,8] fail. We concentrated mainly upon the soliton-like TW solutions and their description, yet similar tools can be applied to the description of kink-like, periodic and multi-periodic solutions as well.
In contrast to the classical self-similarity methods, ansatz-based methods [7,8,9,10] enable to obtain the exact solution with required properties an asymptotics (such e.g. as the periodic, soliton-like and kink-like TW solutions). This is due to the fact that within these methods solutions are defined as algebraic combinations of certain set of functions (elementary of special ones), closed with respect to differentiation and algebraic operations, being represented in a form, enabling to control some features and the asymptotic behavior of the solutions to be described. Yet efficient employment of the direct algebraic balance methods is related to one's ability to solve large systems of nonlinear algebraic equations. In addition to this purely technical obstacle, there is a number of evidence that none of the already known version of the balance algebraic method is fully universal. For this reason the authors put forward the method of the asymptotic description of the soliton-like TW solutions by an infinite series of exponential functions. The effectiveness of this method, connected with some general properties of the set of exponential functions [17], is demonstrated by its application to the GBE. Previous qualitative investigations showed the existence of such solutions for the specific values of the parameters and further on it was shown that these solutions cannot be described as finite algebraic combination of exponential functions.
Let us say a few words about the problem of the truncated series (24) convergence to the exact solution. A uniform convergency of a series of exponential functions to an exact solution can be easily stated e.g. for the nonlinear wave equation, because the expressions for the coefficients a n and b n are known in this case [13]. Unfortunately, for GBE the recurrent expressions for indices a n , b n of the decomposition (24) are not known and the problem of the series convergence remains still open.
The method of approximation of solitary wave regimes, described in Section 3, is relatively easy to perform in case when the factorized system is two-dimensional. In the less trivial cases, when the factorized system is multidimensional, the homoclinic trajectories representing the soliton-like solutions become much more complicated. As an example let us mention a Shil'nikov homoclinic trajectory associated with the saddle-focus of a three-dimensional dynamical system [18,6]. In order to deliver an approximated description of such trajectory following the methodology outlined in Section 3 of this study, it is necessary to incorporate wider class of functions. In earlier attempts at describing the Shil'nikov homoclinics [19], a combination of the exponential and trigonometric functions was used. The approximation employed in [19] was not very effective because of the great technical difficulties encountered. Since then a number of powerful programs for the symbolic calculations were developed, and it is our hope that approximation similar to that used in [19] actually can be performed technically much more effectively.