Analysis of the dynamics of new models of nonlinear systems with state variable damping and elastic coefficients.

This paper, is an analysis of the dynamics of new models of nonlinear systems in which the state damping variables with elastic coefficients, given by functions c cos ⁡ ( p x ) , c sin ⁡ ( p x ) , c cos ⁡ ( p x ˙ ) and c sin ⁡ ( p x ˙ ) are investigated in their autonomous and excited states. They exhibit periodic regions of stability and instability in their autonomous states and a rich dynamic behavior. The analysis of limit cycles shows the presence of isolated curves around the origin (0.0), which explains the presence of periodic solutions (limit cycles). The dynamics obtained allows to describe qualitatively the cardiac activity (artificial pacemaker). A chaos analysis shows the appearance of regular and chaotic behaviors. These studies allowed us to show the effect of the damping of the state variable and the elastic coefficients on the dynamics of these models. The presence of analog functions makes the experimental study complex. An implementation based on microcontroller simulation technology has been proposed. The microcontroller results are consistent with the numerical results.

Introduction

Owing to their broad importance in many fields and their numerous engineering applications, nonlinear oscillators have recently attracted the attention of a large number of researchers due to their extremely rich dynamics (Joshi, 2021; Dashkovskiy and Pavlichkov, 2020; Ahmed et al., 2017; Cheng and Zhan, 2020; Kudryashov, 2018; Tang et al., 2020; Dashkovskiy and Pavlichkov, 2020; Fonkou et al., 2022; Han et al., 2019; Ramirez et al., 2020; FitzHugh, 1961; Rahman et al., 2021a, Rahman et al., 2021b, Rahman et al., 2021c, Rahman et al., 2021d). Their fields of application are between others: seismology, communication and neurophysiology (Rahman et al., 2021a, Rahman et al., 2021b, Rahman et al., 2021c, Rahman et al., 2021d; Lucero and Schoentgen, 2013; Rowat and Selverston, 1993; Balachandran and Kandiban, 2009). These oscillators exhibit rich dynamics among which limit cycle oscillations of sinusoidal and relaxation nature, since one of their important characteristic is their capacity to present limit cycle behaviors which is an important criterion in the characterization of the artificial pacemaker (Steeb, 1977; Hochstadt and Stephan, 1967; D'Heedene, 1996; Steeb and Kunick, 1987; Steeb et al., 1983). When subjected to an external periodic excitation, numerical studies and singular point analysis have revealed chaotic behaviors, allowing the analysis of phenomena such as control and cardiac activity with numerous technological applications. (Steeb and Kunick, 1982; Steeb and Kunick, 1983; Forger, 1999; Enrique et al., 2020; Rahman et al., 2019; Kai and Tomita, 1979; Rahman et al., 2021a, Rahman et al., 2021b, Rahman et al., 2021c, Rahman et al., 2021d; Van der Pol and Van der Mark, 1926; Van der Pol and Van der Mark, 1928; Alhasnawi et al., 2021).

The Van der Pol oscillator is one of these systems, it is nonlinear and its differential equation is the second order. It is a self-sustained oscillator because it maintains its oscillations by itself. Under certain conditions, it also exhibits the very rich dynamical behaviors (Van der Pol and Van der Mark, 1928; Alhasnawi et al., 2021; Chedjou et al., 1997; Makouo and Woafo, 2017; Han et al., 2018; Simo and Woafo, 2012; Bao et al., 2018; Jasim et al., 2020; Han and Bi, 2012; Ma et al., 2021; Grudzinski and Zebrowski, 2004; Magnitskii and Sidorov, 2004). The refs. (Van der Pol and Van der Mark, 1926; Van der Pol and Van der Mark, 1928) have shown that the qualitative characteristics of the heart action potential are closely related to the dynamic response of the classical Van der Pol oscillator (VdP), because it is the starting point for modeling heartbeats.

In the work of Refs. (Steeb and Kunick, 1987; Steeb et al., 1983), the authors study the dynamic behavior of anharmonic systems with limit cycles in which the potentials are only dependent variables. The models used are modifications of the Van der Pol oscillator, the appearance of chaotic dynamics and limit cycle dynamics is due solely to the effect of potentials and the damping coefficient.

In the present work, an analysis of the dynamics of four new models of nonlinear systems under external periodic disturbance is proposed: analytical study and experimental simulation by microcontroller. These models have state variable damping and elastic coefficients. They are described by equations (1) and (2).

However, periodic sine or cosine function of system variable may introduce some coexisting attractors based on offset boosting and therefore self-reproducing chaotic system may be devised.

Faced with this, the following question was asked: how can sinusoidal functions show monotonic solutions of limit cycles and the onset of chaotic dynamics under external excitation? To answer this question, the present work has been organized as follows: Section 2 presents a description of the models used. Section 3 presents an analysis of their dynamics in their autonomous states. Section 4 is devoted to the study of the systems in the excited state. Section 5 shows the results through the practical implementation using an Arduino Uno board. Section 6 ends with the conclusion.

Description of the models

Four new nonlinear oscillator models in which the variables have been substituted by the terms with state variable damping and elastic coefficients given by functions , , and are used. Eq. (1) and Eq. (2) present their dynamics With

With , and p, the characteristic parameters of systems. To maintain the self-oscillating character of the systems, these coefficients must be positive. For an expansion of the functions and up to the order two, the so-called Van der Pol equation is obtained, which is important for triode oscillators necessary for the characterization of the heart (Van der Pol and Van der Mark, 1928). Qualitatively, the quantity represents the friction, if this one is very , one obtains an autonomous oscillatory character and the fixed point unstable tends with the complex eigenvalues of the Jacobean matrix towards a repulsive focus. On the other hand, if it is very strong the oscillatory character is not autonomous and the fixed point stable tends towards an attractive focus with the complex eigenvalues of the Jacobean matrix. Thus influence on the more or less sinusoidal character of these systems.

Quantitatively, the frequency character of the system is linked to the quantity b. Thus, when b increases, the frequency of the oscillations increases and their amplitude decreases. On the other hand, the decrease in b leads to a decrease in the frequency of the oscillations but an increase in their amplitude. Being the natural frequency of the systems, b is chosen so that the systems evolve with a frequency of which is the mean frequency of a healthy heart (Van der Pol and Van der Mark, 1926). Where p is a constant. Using these very rich dynamics, the artificial pacemaker can be a very important application of these models.

Analysis of dynamics in their autonomous state

Analysis of stability

Equation (1) can be rewritten as follows: Obtaining the equilibrium points be done by solving the system below:

For the resulting fixed point is . The Jacobean matrix obtained from this point . By solving , the resulting characteristic equation is and the discriminant is .

However, for , the fixed point is given by . Its Jacobean matrix is . The characteristic equation obtained by solving , is and its discriminant .

It can be seen that: By considering the Eq. (2), it can still be written by:

For and , unstable and is stable.

For and , if , is unstable. However, if , is stable. For , is stable however if , is unstable.

For and , if , is unstable, if , is stable. However, is stable.

From equation (4), we obtain for , the fixed point . Its Jacobean matrix is . By solving , its characteristic equation gives and its discriminant .

However, for , the fixed point is . The Jacobean matrix is . The characteristic equation obtained by solving is given by and the discriminant by . is the identity matrix.

Thus we have:

For and , and are stable.

For and , if , is stable. However, if , is unstable. For , is stable and unstable if .

Be it Eq. (1) where Eq. (2) we find that when, the equilibrium point is stable and the models converge. On the other hand when , the equilibrium point is unstable and the models diverge.

The only critical point being the origin (zero coordinates), a solution of period one exists. In the (x, y)-plane, this solution can describe limit cycle dynamics

Analysis of the dynamics of limit cycles

Recent studies on the nature of integral curves have shown the presence of singular points in systems of differential equations (Chaté, 1994; Bendixson, 1901). Generally, they are observed in conservative systems. However, for some non-conservative systems, closed trajectories or limit cycles could be approached. Eqs. (6) and (7) describe the dynamics of a nonlinear system of order 2 And

Then a limit cycle is necessarily formed in a domain whose limits are respectively the innermost and most external circles of radii and which touch the contact curve (Hayashi, 1964).

By considering the models (3) and (5), since the origin is an only singular point, a family of concentric circles around this origin is considered, that is to say:

In this paragraph, it will be demonstrated that the models described by systems (3) and (5) exhibit limit cycle oscillations under origin (0, 0).

By considering Eq. (1), for , starting from the given Eq. (8), we obtain:

When the polar coordinates given by Eq. (10) are introduced in Eq. (9), and apart from , the expressions of the radii and are presented in Eq. (11):

Still in considering Eq. (1), but this time with , the expressions of the radii and are given by:

When the dynamics of system are described by Eq. (2) with , by substituting Eq. (10) in Eq. (9) and apart from , the expressions of the radii and are given by:

Still by taking Eq. (2) with , the expressions of the radii and are given by:

Consequently from expressions given in (11), (12), (13), and (14), a ring domain centered at the origin exists and contains all possible limit cycles; these boundaries are the radii and . They represent the boundaries of the innermost and outermost circles, respectively, that touch the contact curve. These expressions show that if a limit cycle occurs around the origin, it must necessarily be in a ring domain whose center is at the singular point.

In Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, the temporal traces and phase portraits have been plotted for the initial conditions and . This allows us to observe that these systems are nonlinear and self-excited and show limit cycle behaviors of sinusoidal and relaxation nature. The limit cycles are almost isolated curves so the trajectories converge around the origin.

Figure 1

Time series and phase portrait obtained from Eq. (1) with a = 0.5, b = 39.4384, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

Figure 2

Time plots and phase portrait obtained from Eq. (1) with a = 0.5, b = 39.4384, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

Figure 3

Time plots and phase portrait obtained from Eq. (1) with a = 0.001, b = 39.4384, c = 5, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

Figure 4

Time plots and phase portrait obtained from Eq. (2) with a = 0.001, b = 39.4384, c = 5, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

Figure 5

Time plots and phase portrait obtained from Eq. (2) with a = 0.001, b = 39.4384, c = 5, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

The limit cycles obtained from Eq. (1) (case of ) are presented by figure.

For , the temporal trace and phase portrait are given in Fig. 1.

On the other hand, for , we have Fig. 2.

Still by considering Eq. (1) (case of ), the time serie and the phase portrait are given in Fig. 3.

By considering the Eq. (2) (case of ), the limit cycles are given in Fig. 4.

When we consider Eq. (2) (case of ), the time plot and phase portrait are given in Fig. 5.

It can be seen in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 that the models used are nonlinear and self-excited and exhibit limit cycle dynamics of sinusoidal and relaxation nature qualitatively comparable with those exhibited by artificial pacemakers (Van der Pol and Van der Mark, 1928). In addition, the sinusoidal oscillations evolve with a frequency of 1 Hz for where and .

Analysis of dynamics in their excited state

Chaos analysis

In the anharmonic systems with limit cycles, it has been shown that the variation of their parameters as well as those of the external force can lead to chaotic behaviors. The Refs. (Steeb and Kunick, 1987; Steeb et al., 1983) have demonstrated the influence of potentials and damped coefficients on the appearance of chaotic dynamics in nonlinear systems. In this study, our objective is to see the effect of state variable damping and elastic coefficients on the appearance of chaotic dynamics in nonlinear systems under the effect of an external excitation of sinusoidal nature . This is done by plotting the bifurcation curves each associated with its maximum Lyapunov exponent with the amplitude of the external excitation as parameter. These curves (bifurcation and maximum Lyapunov exponent) are presented in Figure 6, Figure 7, Figure 8, Figure 9. On these figures, we observed the presence of chaotic and multiperiodic dynamics. However, for small values of the parameter , the Lyapunov exponent tends towards 0. This behavior reflects the presence of quasiperiodic dynamics. A few are shown in Fig. 6. In Fig. 11 it only about the presentation of some phase portraits illustrating the chaotic states.

Figure 6

Bifurcation tree and Lyapunov exponent obtained when Eq. (1) is excited (case of ), with ω = 4.8 rad.s−1, a = 0.5, b = 0.25, c = 2 and p = 5.

Figure 7

Bifurcation tree and Lyapunov exponent obtained when Eq. (1) is excited (case of ), with ω = 4.8 rad.s−1, a = 0.5, b = 0.25, c = 8 and p = 5.

Figure 8

Bifurcation tree and Lyapunov exponent obtained when Eq. (1) is excited (case of ) with ω = 4.8rad.s−1, a = 0.001, b = 50.25, c = 5, and p = 5.

Figure 9

Bifurcation tree and Lyapunov exponent obtained when Eq. (2) is excited (case of ) with ω = 4.8rad.s−1, a = 0.001, b = 50.25, c = 5, and p = 5.

Figure 11

Chaotic phase portraits for x(0)=0.01 and dx/dt(0)=0.04.

When Eq. (1) is excited (case of ), the bifurcation trees and maximal Lyapunov exponents are given in Figure 6, Figure 7.

When , the information given by Fig. 6 are: a period-2T orbit is obtained, this behavior disappears and gives way to the quasiperiodic behavior for . A second orbit of period-13T appears for , then disappears to give a quasiperiodic orbit for . For , a new orbit of period-11T is observed, the latter disappears and gives a quasiperiodic behavior for . When , the orbit of period-1T settles.

On the other hand, when , on the Fig. 7, an orbit of period-1T appears for then disappears and gives way to the quasiperiodic orbit for . For the first orbit of period-34T appears. When the value of goes from 1.625V to 2.5V, a second quasiperiodic orbit appears and disappears, giving way to the period-32T orbit for . For , quasiperiodic orbits appear, a second orbit of period-30T appears for . This orbit disappears and the first chaotic behavior appears for . This first chaotic behavior is replaced by a second orbit of period-27T when the value of increases from 3.776V to 4.385V. Then a second chaotic orbit appears for . Subsequently, the alternations between periodic and chaotic orbits are observed. For an orbit of period 25T is observed. When a third chaotic orbit is obtained. Then a new orbit of period-23T for . For a fourth chaotic orbit appears. Between 5.975V and 6.642V, an orbit of period-21T is also observed. A fifth chaotic orbit between 6.642V and 6.696V. For , a new orbit of period-19T appears. When a sixth chaotic orbit appears. Between 7.396V and 8.07V an orbit of period-17T appears. Between 8.07V and 8.14V a seventh chaotic orbit is observed. Between 8.14V and 8.769V an orbit of period-15Tis obtained. For an eighth chaotic orbit appears. For an orbit of period-14T appears. Between 9.46V and 9.53V, a ninth chaotic behavior is obtained. Between 9.53V and 10.18V the orbit of period-13T appears. For , a tenth chaotic orbit is obtained. Between 10.27V and 10.89V, an orbit of period-11T appears and between 10.89 V and 11.02V an eleventh chaotic behavior is obtained. For , the orbit of period-9T is observed and for , a twelfth chaotic behavior. When varies between 11.88V and 12.4V, there is the appearance of a period-7T orbit. However, between 12.4V and 13.35V, the quasiperiodic orbit reappears. From , only the orbit of period 1T evolves.

When we excite Eq. (1) (case of ), we have in Fig. 8 the bifurcation tree and its associated maximal Lyapunov exponent.

For this case, a chaotic orbit is obtained for , it then disappears and gives way to the period-3T orbit for . For a new chaotic orbit appears. When the value of goes from 36.5V to 50V, a second orbit of period-2T appears, then disappears, giving way to a new chaotic orbit for . For , A new appearance of the period-3T orbit. A third chaotic orbit appears for . This orbit disappears and the first quasiperiodic orbit appears for . Then, we observe a succession between chaotic and quasiperiodic orbits for . Then reappearance of the period-2T orbit for . When , a new succession between chaotic and quasiperiodic orbits appears.

The bifurcation curves and maximal Lyapunov exponent are given respectively in Fig. 9 for and Fig. 10 for .

Figure 10

Bifurcation tree and Lyapunov exponent obtained when Eq. (2) is excited (case of ) with ω = 4.8rad.s−1, a = 0.001, b = 50.25, c = 5, and p = 5.

For these cases (Fig. 9 and Fig. 10), an alternation between chaotic orbits and periodic orbits is observed.

To illustrate the chaotic behaviors presented by the bifurcation curves and Lyapunov exponents of Figure 7, Figure 8, Figure 9, Figure 10, some phase portraits are given in Fig. 11.

Simulation based on microcontroller technology

Microcontroller simulation is one of the simplest because it is performed using a computer in which suitable simulation software is installed. Recent work has shown that it can be used to do electronics using just a programming language (Fonkou et al., 2021a, Fonkou et al., 2021b, Fonkou et al., 2021c). Experimentally, it is less cumbersome and requires significant period savings for component wiring, a higher degree of integration, and lower power consumption.

A microcontroller is an integrated circuit that gathers the essential elements of a computer. Recent work has shown that microcontrollers can be used to simulate nonlinear differential equations (Fonkou et al., 2021a, Fonkou et al., 2021b, Fonkou et al., 2021c) by producing simple and complex electrical signals. Compared to conventional electronic systems based on separate components, they can reduce the size of the components, as well as the cost of the products. In this work, we will use an Arduino microcontroller from the ATMEGA family (Arduino Card). Its structure is shown in Fig. 12.

Figure 12

Device mounted on Proteus allowing the observation of the signal on its oscilloscope.

To produce the electrical signal described by Eq. (1) and Eq. (2), these equations are discretized. Subsequently, a simple digital code using appropriate software (here the Arduino 1.8.9) is written and loaded into the microcontroller. The signal obtained from Proteus software with its own oscilloscope is shown in Fig. 12.

To produce this signal experimentally from a real oscilloscope, the practical configuration presented in Fig. 13 is used. It consists of a computer for inserting the digital code into the microcontroller and the oscilloscope for viewing the signal delivered by the latter. The real electrical signals obtained are given in Fig. 14, Fig. 15, and Fig. 16.

Figure 13

Experimental device allowing the visualization of the signal.

Figure 14

Experimental results obtained from Eq. (1), with a = 0.5, b = 0.25, c = 2, p = 5 (Fig. 14a) and a = 0.5, b = 0.25, c = 8, p = 5 (Fig. 14b) for x(0)=0.01 and dx/dt(0)=0.04.

Figure 15

Experimental results obtained from Eq. (1), with a = 0.001, b = 50.25, c = 5, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

Figure 16

Experimental results obtained from Eq. (2) (Fig. 16a: case of ) and (Fig. 16b: case of ) with a = 0.001, b = 50.25, c = 5, p = 5, x(0)=0.01, and dx/dt(0)=0.04.

For , we have Fig. 14.

For , we obtain in figure.

These curves (Fig. 14 and Fig. 15) are in agreement with those obtained respectively in the Fig. 2 and Fig. 3.

On the other hand, when we use the Eq. (2), the dynamics are given by figure.

The curves presented in the figures above (Fig. 16a and Fig. 16b) show a good agreement with those obtained numerically (Fig. 4 and Fig. 5).

For the non-autonomous case, an experimental investigation similar to the autonomous case was carried out using the Runge Kutta algorithm of order 4. Two Arduino Uno boards are used practically. The first board is used to produce the electrical excitation signal . The Arduino code used to produce this signal and the complete flowchart used to program this signal are given in Table 1.

Table 1

The Arduino code and the complete flowchart for programming the signal

The Arduno code used to produce the signal E_0sin⁡(ωt)	The complete flowchart used to program the signal E_0sin⁡(ωt)
// constants	Image 1
float ts=(float)0.00016;
float E0=(float)4.37;
float w=(float)4.8;
float t =(float)0.0;
// vriables
float E0, ts;
void setup () {
pinMode (0, OUTPUT);
pinMode (1, OUTPUT);
pinMode (2, OUTPUT);
pinMode (3, OUTPUT);
pinMode (4, OUTPUT);
pinMode (5, OUTPUT);
pinMode (6, OUTPUT);
pinMode (7, OUTPUT); }
void loop() {
t=t+ts;
E=E0*sin(w*t);
PORTD =(E-20)*45.5;
t = t; }

The second card is used for programming the self-excited oscillator. The Arduino code and the complete flowchart used to program this oscillator are presented in Table 2.

Table 2

The Arduino code and the complete flowchart for programming the self-excited oscillator.

The Arduino code	The complete flowchart
//constants	Image 2
float h=(float)0.1;
float a=(float)0.5;
float b=(float)0.25;
float c=(float)8.;
float p=(float)5.0;
//initial conditions
float x1 =(float)0.01;
float y1 =(float)0.04;
float t1 =(float)0.0;
// vriables
double atan( double t1);
double sin( double t1);
double cos( double t1);
float x,x2,y2,t,d;
float y,z,r,tt1,xx1;
float L1,L2,L3,L4;
float M1,M2,M3,M4;
void setup( ) {
pinMode (0, OUTPUT);
pinMode (1, OUTPUT);
pinMode (2, OUTPUT);
pinMode (3, OUTPUT);
pinMode (4, OUTPUT);
pinMode (5, OUTPUT);
pinMode (6, OUTPUT);
pinMode (7, OUTPUT);
pinMode (A0, INPUT); }
void loop() {
L1=y1;
M1=-(a-c*cos(p*x1))*y1-b*x1;
z=x1+(h/2)*L1;
r=y1+(h/2)*M1;
t=t1+(h/2);
L2=r;
M2=-(a-c*cos(p*z))*r-b*z;
z=x1+(h/2)*L2;
r=y1+(h/2)*M2;
t=t1+(h/2);
L3=r;
M3=-(a-c*cos(p*z))*r-b*z;
z=x1+h*L3;
r=y1+h*M3;
t=t1+h;
L4=r;
M4=-(a-c*cos(p*z))*r-b*z;
x1=x1+(h/6)*(L1+2*L2+2*L3+L4);
y1=y1+(h/6)*(M1+2*M2+2*M3+M4);
t1=t1+h;
PORTD=(x1+20)*55.5;
x1 = x1;
y1 = y1;
t1 = t1; }

To obtain the result, the signal produced by the first Arduino board is inserted into an input port of the second Arduino board to excite the self-excited oscillator (Eq. (1) or Eq. (2)) programmed in this second card, using an appropriate line of code.

The device cabled on Proteus is shown on Fig. 17 and its practical device is presented in Fig. 18. As shown in Figure 19, Figure 20, only chaotic states have been presented.

Figure 17

Microcontroller device representing our exciter signal coupled contains Eq. (1) or Eq. (2).

Figure 18

Experimental device allowing the visualization of the signal.

Figure 19

Chaotic dynamics obtained numerically and by microcontroller, with x(0)=0.01 and dx/dt(0)=0.04.

Figure 20

Chaotic dynamics obtained numerically and by microcontroller, with x(0)=0.01 and dx/dt(0)=0.04.

The chaotic dynamics obtained for and are presented in Fig. 19. These curves show a good agreement between the experimental results and the results obtained numerically.

For and the chaotic dynamics obtained experimentally are presented in Fig. 20, and a good agreement between the experimental results and those obtained numerically is observed.

Conclusion

This work was devoted to the study of four new models of nonlinear systems with state variables damping and elastic coefficients defined by functions , , , and . The study of limit cycles has shown that these systems, to some degree of approximation, can describe the electrical activity of the heart. The study of their dynamics in their autonomous and non-autonomous states has shown that the values of the coefficients influence the nature of the stability/instability of these systems. The production of real electrical signals by discretizing differential equations and writing digital codes in the microcontroller show that the latter can be used for the design of nonlinear oscillator models. In doing so, it appeared that when these nonlinear systems are subjected to external and periodic sinusoidal excitation, they presented periodic, quasiperiodic, and chaotic behaviors. Since these models are essentially autonomous oscillators, they continue to deliver a signal if the action of ceases.

Declarations

Author contribution statement

Rodrigue Fatou Fonkou, Dr: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper. Patrick Louodop, Associate Professor: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper. Pierre Kisito Talla, Professor; Paul Woafo, Professor: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

No data was used for the research described in the article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.