Discussion Paper No.269
Love Dynamics Between Cautious Romeo and Juliet
With Delays in Self and Cross Reactions
Akio Matsumoto
Chuo University
Ferenc Szidarovszky
University of Pécs
July 2016

INSTITUTE OF ECONOMIC RESEARCH
Chuo University
Tokyo, Japan

Love Dynamics Between Cautious Romeo andJuliet with Delays in Self and Cross Reactions
Akio Matsumotoy

Ferenc Szidarovszkyz

Abstract
We investigate love dynamics of two individuals in a delay Romeo andJuliet model in which both are assumed to be cautious, the most naturalromantic style. The local stability analysis proves …rst that the steadystate is fairly stable when there are no delays and second that solving
the characteristic equation generates a set of positive delays for which thesteady state loses stability. Through numerical analysis, we con…rm thefollowing three main results: (1) cyclic oscillations of love feeling emerge
via Hopf bifurcation; (2) multiple delays cause the double edge e¤ect
implying that alternation of stability and instability repeatedly occurs;(3) complicated dynamics involving chaotic oscillations emerges and thenmerges to a limit cycle as the length of one delay increases with …xed
values of the other delay.

Keywords: Love dynamics, Delay di¤erential equations, Cautious
lovers, Hopf bifurcation, Double edge e¤ect

The authors highly appreciate the …nanical supports from the MEXT-Supported Program
for the Strategic Research Foundation at Private Universities 2013-2017, the Japan Society
for the Promotion of Science (Grant-in-Aid for Scienti…c Research (C), 25380238, 26380316,
16K03556) and Chuo University (Joint Research Grant). The usual disclaimer applies.
y Professor, Department of Economics, Senior Researcher, International Center for further
Development of Dynamic Economic Research, Chuo University, 742-1, Higashi-Nakano, Hachioji, Tokyo, 192-0393, Japan. akiom@tamacc.chuo-u.ac.jp
z Professor, Department of Applied Mathematics, University of Pécs, Ifjúság u. 6., H-7624,
Pécs, Hungary. szidarka@gmail.com

1

1

Introduction

This paper considers love a¤air dynamics with two distinct delays and extends
the study of Matsumoto and Szidarovszky (2016) that deals with love dynamics
with one delay. It presents a new characterization of love evolution by whicha wide spectrum of love dynamics can be constructed, ranging from monotonic
convergence to complicated oscillations including chaotic phenomena.
Di¤erential equations have been used to describe dynamic phenomena in various …elds of science; Navier-Stockes equations in physics, Schrödinger equationin quantum chemistry, Lotka-Voltera equations in biology, price dynamic system following the law of supply and demand in economics, to name only a few.
In social psychology, Strogatz (1988) could be the …rst to model the dynamics
of romantic feelings of two individuals named Romeo and Juliet with two di¤erential equations. It is designed to evolve a love a¤air between Romeo and Juliet
through two routes. It may be very natural that one’s love grows more strongly
when the other expresses his/her love more passionately. It is, on the contrary,also possible that one’s feeling gets gradually down when this individual is a…ckle lover. So the …rst route describes a change in one’s own level of love feeling a¤ected by the partner’s level. We say that time evolution of the love a¤air
through this route is nurtured in a cross-reaction process. Second, one’s own
level can change even in absence of the partner. If Juliet disappears by somereasons, Romeo gradually loses his a¤ections as time goes on or he could live inold days memories, preserving his a¤ections to Juliet. We say that love a¤ection
through this route evolves in a self-reaction process. Further, how love grows
or decays depends on the romantic styles of Romeo and Juliet. According toStrogatz (1988), the styles are classi…ed into four speci…cations, "eager beaver"
if both of the self and cross reactions are positive, "cautious lover" if the self
reaction is negative and the cross-reaction is positive, "narcissistic nerd" if theself-reaction is positive and the cross-reaction is negative and "hermit" or "stoic
lover" if both are negative.
Taking into account of the real-life fact that the majority of the populationis the cautious lover, we con…ne our main analysis to a case in which Romeoand Juliet are cautious lovers. Further, real life observations often tell us two
distinguished evolutions of love after they meet for the …rst time. One is thatlove stories develop monotonically and, sooner or later, stay at their plateau fora long time (i.e., convergent dynamics). The other is that love stories arrive at aoscillating regime in which the feelings of the individuals sometimes fall into thebottomless pit, some other times rise up to the peak of happiness and these ups
and downs repeat themselves (i.e., cyclic dynamics). In the existing literature,
two types of models are known to explain these observations. On one hand, it
has been demonstrated that a variant of Romeo and Juliet model can interpret
various types of love stories described in "Pride and Prejudice," (Rinaldi et al.
(2014)) "Gone with the Wind" (Rinaldi at al. (2013a)) and "Beauty and theBeast" (Rinaldi at al. (2013b)) by focusing on the appeals of the individuals.
On the other hand, following natural phenomena that interactions between theindividuals could be delayed, Liao and Ran (2007), Son and Park (2011) and2

Bielczyk et al. (2013) focus on a delay in a nonlinear cross-reaction process
and show the birth of cyclic dynamics through Hopf bifurcation. This …nding
indicates that the nonlinear delay model of Romeo and Juliet explains various
oscillatory dynamics of two individuals’ romantic feelings. In the existing literature, however, not much has been revealed with respect to multiple delays
in reactions of the individuals to stimuli. In this paper we construct a class oflove dynamic models of cautious individuals in which there are delays in bothself-reaction and cross-reaction of love accumulation process and study the delaye¤ects on love dynamics analytically as well as numerically.
In what follows, Section 2 presents the basic love dynamic model with nodelays. In the …rst half of Section 3, we review Liao and Ran (2007) in which
two delays are assumed in the cross-reaction process. Then in its latter half,
we additionally introduce the self-reaction delay to their model and analyticallyderive a stability switching curve on which stability is switched to instability orvice versa. Section 4 conducts numerical simulations and demonstrates that thetwo delay model can generate rich dynamics to describe "many couples, many
ways to express their love and a¤ection". Section 5 concludes the paper andprovides directions of future research.

2

Basic Model

We now construct a Romeo and Juliet model without delays, which is called abasic model. Let x(t) and y(t) denote levels of the romantic feelings of Romeo
to Juliet and that of Juliet to Romeo at time t, respectively, if x(t) > 0 andy(t) > 0. The negativity of these state variables represents a level of the nonromantic or dislike feelings such as antagonism and disdain. Then the rates ofchange describe the feeling accumulation processes of the individuals and areassumed to have the following forms according to Rinaldi (1998):x(t)_= Ox (x(t)) + Rx (y(t)) + Ix
(1)y(t)_ = Oy (y(t)) + Ry (x(t)) + Iy
each of which is composed of three terms, oblivion denoted as Oz ; for z =x; y; return by Rz and instinct by Iz . First, Oz gives rise to a loss of interest
that describes the self-reaction process and depends on his/her own feeling level.
It characterizes decay of love at disappearance of the partner in the self-reaction
process. Second, Rz is a source of interest and describes the reaction of individual z to the partner’s love in the cross-reaction process. Lastly, Iz is alsoa source of interest and describes the reaction of individual z to the partner’s
appeal re‡ecting physical, …nancial, educational, intellectual, well-born properties, the family background, etc. We adopt the following forms of these reaction
functions:
Assumption 1: Ox (x) =Assumption 2: Rx (y) =x x;x > 0 and Oy (y) =x tanh(y) and Ry (x) =

3

y y;

y tanh(x):y > 0:

Assumption 3: Ix =x Ay ; Ay

0 and Iy =y Ax ;Ax

0:

Assumption 1 con…nes attention to the case where the a¤ection vanishes
exponentially in the absence of the partner. Assumption 2 implies that thereaction is determined by a product of two terms, the reaction coe¢ cient andthe hyperbolic reaction that is positive, increasing, concave and bounded fromabove for positive values and is negative, increasing, convex and bounded frombelow for negative values. The love a¤ection of individual z is encouraged or
discouraged by the partner according to whether z > 0 or z < 0. Assumption
3 implies that individuals have time-invariant appeal. Non-zero Az a¤ects not
only the location of the steady state but also love dynamics as mentioned in
the Introduction. In this study, since we con…ne attention to the delay e¤ects
on love evolution, Az = 0 is assumed to simplify the analysis.1 Under theseassumptions, equations (1) are reduced to a more speci…c system,
x(t)
_
=

x x(t) +x tanh[y(t)];

y(t)
_ =

y y(t) +y tanh[x(t)];(2)
where z is the self-reaction or forgetting coe¢ cient and z is the cross-reaction
or return coe¢ cient. The basic structure of system (2) is the same as that of
the model Strogatz (1988) proposes. The minor di¤erence is that a linear return
function of Strogatz is replaced with a nonlinear hyperbolic function. We callit a basic Romeo and Juliet model. Matsumoto and Szidarovszky (2016) havealready demonstrate that it has a unique zero steady state if x yx y and
two more non- zero steady states if x y > x y : Accordingly, let (x0 ; y0 ) be
the zero steady state and (xk ; yk ) the positive steady state if k = 1 and thenegative steady state if k = 2:
To examine stability of the steady states, system (2) is linearized,
x(t)
_
=

x x(t) +kx dy y(t);ky dx x(t)y y(t)

(3)
y(t)
_ =
where
dkx =

d tanh(y)d tanh(x)and dky =
:
dx
dy
x=x
y=y
k

k

t

t

Substituting exponential solutions x(t) = e u and y(t) = e v into system (3)yields the corresponding characteristic equation
1
0
k
+ x
x dy
A=0
det @
k
d
+
y
y x
1 On the other hand a number of papers of Rinaldi treats the appeal as an ingredient factor

for love evolution. Since our approach could complement Rinaldi’s approach, a delay modelwith the appeal could bring about more fruitful results.

4

or
2

+( x+

y)

+

x y

k k
x y dx dy = 0;(4)

where x + y > 0 always by Assumption 1. It is also demonstrated in Matk k
sumoto and Szidarovszky (2016) that x yx y dx dy > 0 for k = 1; 2 andkkdx = dy = 1 for k = 0: Hence we have the following result, which is a summaryof Theorem 2 of Matsumoto and Szidarovszky (2016);
Theorem 1 The non-zero steady states are stable nodes while the zero-steadystate is a saddle if x y > x y and stable otherwise.
The following numerical simulations con…rm Theorem 1 and are done with
y = 1, x = y = 3=2 in Figure 1(A) and x = y = 1=2 in Figure 1(B).
It is clearly seen that (x0 ; y0 ) is a saddle and (xk ; yk ) for k = 1; 2 are stable in
Figure 1(A) whereas the unique steady state (x0 ; y0 ) in Figure 1(B) is a stablenode.2
x =

(A)

x y >

x y >0(B)x y >x y >0Figure 1. Stability properties of the basic model3

Two Delay Model

We introduce time delays in the interactions of the individuals in the basicmodel and consider how the delays a¤ect love evolutions just described aboveby the nonlinear Romeo and Juliet model without any delays. Matsumoto
and Szidarovszky (2016) have already con…rmed the appearance of cyclic love
dynamics in the model with one delay. Henceforth we investigate the e¤ects
caused by two distinct delays in this study.
2 In the same way, it can be checked that the zero-steady state is stable in the cases of

0>

x y >

x x and 0 >

x x >

x y:

5

3.1

Liao and Ran ModelWe …rst review the dynamic results provided by Liao and Ran (2007) that
introduce two delays x and y into the cross-reaction processes of the basicmodel,x(t)_=x x(t) + x tanh[y(t
y )];
y(t)
_ =

y tanh[x(tx )]y y(t):The steady states of this model are the same as the ones of the basic model.
To examine the stability of the steady states, the model is linearized in theneighborhood of the steady state (xk ; yk )
x(t)
_
=

k
x dy y(tx x(t) +

y );

(5)
y(t)
_ =

k
y dx x(tx )]y y(t):

Substituting exponential solutions x(t) = e t u and y(t) = e t v into system (5)yields, after arranging the terms, the corresponding characteristic equation
2

+( x+

y)

+

x y

k k
x y dx dy e

=0

(6)

with = x + y : Although Liao and Ran (2007) introduce two delays into thecross-reaction processes, dynamics generated by their model is essentially thesame as the one by the model with one delay in the cross-reaction process since
only the sum of two delays plays a crucial role in determining dynamic behavior
in their model. Hence, following Matsumoto and Szidarovszky (2016), it canbe shown …rst that all pure complex eigenvalues of equation (6) are simple andsecond that the non-zero steady states are locally asymptotical stable. Liao andRao (2007), on the other hand, concern with dynamics of the zero-steady stateand show the following result, which restates their Theorem 1:
Theorem 2 Concerning the zero-steady state, (x0 ; y0 ); no stability switch occurs if x yx y and x y >x y while a stability switch occurs only onceat = 0 if 0 >x y > x y , implying that the zero stationary state is locallyasymptotically stable if < 0 , loses stability if = 0 and bifurcates to a limit
cycle if > 0 where the threshold value of is de…ned as
0 =

with
1
!0 = p
2

rq

1
sin 1
!0

( 2x

( x+

y )! 0x y2 )2 + 4
y

2
x y

2 +
x

2
y

:

A Hopf bifurcation resulting in cyclic dynamics occurs under two conditions,(i) the product of the cross reaction coe¢ cients is larger than the product of theself reaction coe¢ cients and (ii) one of the two individuals is a hermit and the6other is cautious. The …rst one could be possible, however, the second might beunusual, as already been pointed out by Son and Park (2011). To remedy thisde…cit, we will limit the styles of Romeo and Juliet to cautious lovers henceforthand then pursue possible emergencies of cyclic oscillations:
Assumption 2’.

3.2

x > 0 andy > 0:

Delay Romeo and Juliet Model

We introduce multiple delays into the self- and cross-reaction processes in sucha way that each delay has an independent role:
x(t)
_
=

x x(t

s
x) +

c
y )];

x tanh[y(t

(7)
y(t)
_ =

c
x )]

y tanh[x(ty y(t)

where sx > 0 is a self-reaction delay of Romeo3 while cx
0 and cy0cccross-reaction delays of Romeo and Juliet satisfying x + y > 0. The linearizedversion of model (7) is
x(t)
_
=

x x(t

s
x) +

k
y dx x(t

y(t)
_ =

c
x)

k
x dy y(tcy );y y(t);and the analysis of the Jacobian as before reveals that the characteristic equation
is
k k
1
2
( + y ) + x ( + y )e=0(8)
x y dx dy ewhere 1 = sx and 2 = cx + cy : Following the method provided by Gu et al.
(2005), we detect the location of the eigenvalues of this characteristic equationin the same way as in Matsumoto and Szidarovszky (2015).
k k
Since x y
x y dx dy > 0 for k = 1; 2 is shown in Matsumoto and Szidarovszky (2016), = 0 is not a solution of equation (8). Then dividing bothsides of the characteristic equation by ( + y )(6= 0) reduces the left hand sideto1
2
+ a2 ( )ea( ) = 1 + a1 ( )ewhere
a1 ( ) =x

and a2 ( ) =k k
x y dx dy

( +

y)

:

Theorem 1 ensures the stability of the non-zero steady state when there are
no delays (i.e., 1 = 2 = 0). We now examine whether the switching fromstability to instability or vice verse can occur at some positive values of thedelays. To this end, we check if the characteristic equation (8) has a pair ofpure imaginary solutions at the threshold values of the delays. Let us suppose
3 The qualitatively same results will be obtained even if the self-reaction delay of Romeo is
replaced with a self-reaction delay of Juliet.

7

that = i!; ! > 0 where we have the same results under ! < 0 as the solutions
are conjugate. We then havea1 (i!) =
and

x

i

!

k k
k k
y x y dx dy
x y dx dy !
+i:!(! 2 + 2y )!(! 2 + y )a2 (i!) =Their absolute values arex

ja1 (i!)j =andja2 (i!)j =Furthermore,arg [a1 (i!)] =!

k kx y dx dy
q:! ! 2 + 2yh i
3yand arg [a2 (i!)] = tan 1:
2!

Solving a(i!) = 0 directly is not an easy job. However, if the three termsin a(i!) are considered to be three vectors in the complex plane, then solving
it is equivalent to constructing a triangle from the three terms satisfying thefollowing (triangle) conditions,(1) 1

j 1 (i!)j + j 2 (i!)j(2) j 1 (i!)j

1 + j 2 (i!)j(3) j 2 (i!)j

1 + j 1 (i!)jConditions (1) and (2) are written asf (!) = (!

2
x)

!2 +

2
y

k k 2
x y dx dy

2
x)

!2 +

2
y

k k 2
:
x y dx dy

and condition (3) asg(!) = (! +Hence, the triangle conditions are satis…ed if
f (!)

k k 2
x y dx dyg(!):

It is to be noticed thatf (0) = ( x y )2 ; f ( x ) = 0; f ( 1) = 1 and f (!) < g(!):
Di¤erentiating f (!) yields
f 0 (!) = (!

x )(4!

8

2

2! x + 2 2y ):The …rst factor is negative if ! < x and positive if ! >the second factor seems to be ambiguous. Let'(!) = 4! 2x while the sign of

2! x + 2 2yand the roots of '(!) = 0 arex! 1;2 =

q

2
x

4

8 2y

:

For simplicity, we make the following4 :Assumption 4:x 'y such that22
x < 8 y always.

The negative discriminant implies that '(!) 0 for all !: Hence we have8
< < 0 if ! < x ;f 0 (!):> 0 if ! > x :So the domain for ! in which the triangle conditions are satis…ed is de…ned asfollows.
2

Theorem 3 Let B = x y dkx dky . Then the domain is between roots of f (!) =B if B ( x y )2 and between roots of g(!) = B and f (!) = B if B > ( x y )2 :Let the three terms in a(i!) be three vectors forming a triangle. We supposethat j1j ; the absolute value of vector 1, is its base and denote the angle between
j1j and j 1 (i!)j by 1 and the angle between j1j and j 2 (i!)j by 2 : Then by
the law of cosine, we havecos 1 =
andcos 2 =(! 2 +

2
2x )(! +

2 x
(! 2

2
y)
!(! 2 +

( x y dkx dky )22)y2) + ( x y dkx dky )2qy:
2 x 2 ! ! 2 + 2y22x )(! +

4 If 28 2y holds, then the discriminant is nonnegative. Both roots ! 1 and ! 2 are real,x
positive and less than x : So we have8
< > 0 if ! < ! 1 or ! > ! 2 ;'(!):< 0 if ! 1 < ! < ! 2 :Hence the sign of f 0 (!) is determined such as8< < 0 if ! < ! 1 or ! 2 < ! <f 0 (!):> 0 if ! 1 < ! < ! 2 or ! >x;
x.

So domain of ! can be de…ned in various ways. However, it is numerically checked thatdynamics obtained under 2x 8 2y is essentially the same as the one under 2x < 8 2y :9

Solving these two equations for 1 and 2 gives
"
#
2
2
2
2
k k 2
(!
+
)(!
+
)
(
d
d
)
x
y
x
y
x
y
1
1 (!) = cos
2 1 !(! 2 + 2y )and2 (!) = cos

2

14

(! 2

2
2
x )(! +

2 1 2!

2y)( x y dkx dky )2

q
!2 +

2
y

3

5:

Since the triangle may be located above and below the horizontal axis in thecomplex plane, we have two possibilities,arg

1 (i!)e

i! 1

+ 2k

1 (!) =arg

2 (i!)e

i! 2

+ 2n

2 (!) =

and
:

i! k

Using the formula arg[ k (i!)e] = arg[ k (i!)] + arg[e i! k ] and solving
these equations for 1 and 2 yield the threshold values of the delays,
1
farg [ 1 (i!)] + 2(k 1)
1 (!)g1 (!; k) =!
and1farg [ 2 (i!)] + 2(n 1)
2 (!)g2 (!; n) =!
for (k; n) = 0; 1; 2; ::: We can …nd the pairs of ( 1 ; 2 ) constructing the partition
curves consisting of two sets of parametric segments for k; n 0,
L1 (k; n) =+1 (!; k);2 (!; n)for ! 2 [! s ; ! e ]L2 (k; n) =1 (!; k);+2 (!; n)for ! 2 [! s ; ! e ]:andHere [! s ; ! e ] denotes the domain of ! in which the triangle conditions hold.
According to Theorem 2, the left hand extreme value or the starting value ! ssolves f (!) = B if B ( x y )2 and g(!) = B if B > ( x y )2 while the righthand extreme value or the ending value ! e always solves f (!) = B: The next
result con…rms that the segments of L1 (k; n) and L2 (k; n) with …xed k andvarying n smoothly connected to form one continuous curve.5
Theorem 4 With a …xed value of k; the segments of L1 (k; n) and L2 (k; n) form
a continuous curve as n increases.
Notice two issues, one is that characteristic equation (8) has a pair of pure
imaginary roots on these partition curves and the other is that given k and n;
the partition curve divides the ( 1 ; 2 ) region into subregions according to thenumber of the eigenvalues whose real parts are positive. In consequence, sincecharacteristic equation (8) is reduced to equation (4) as 1 = 2 = 0, the steadystate is locally stable for any pair of ( 1 ; 2 ) in the separated region including
the origin.
5 Applying Theorem 1 of Matsumoto and Szidarovszky (2015), we can demonstrate thatthese segments from a continuous curve10

4

Numerical Simulations

In this section we specify the parameter values and perform numerical simulations …rst to con…rm the conditions under which stability switching occurs andthen to detect what kind of dynamics emerges when stability is lost. Under
Assumptions 1, 2 and 2’, we can identify two cases, x y > x y > 0 andx y > x y > 0: These cases are successively considered in the sequel.

4.1

x y >

x

y > 0

Under this parametric conditions, the following two issues have been alreadyshown in the basic or non-delay model:(1) three steady states exist,(2) the non-zero steady states are locally asymptotically stable and the zerosteady state is a saddle.
In this section we are concerned with the delay e¤ects caused by changing
the values of 1 and 2 on dynamics. More precisely, specifying the parametervalues, we numerically examine whether the delays can destabilize the non-zerosteady states and what dynamics can emerge when the stability is lost. As in
Figure 1(A), we make the following:
Assumption 4:

1 =2 = 1 and1 =2 = 3=2

Under Assumption 4, Figure 2 illustrates the segments of L1 (k; n) in red andL2 (k; n) in blue with the values of ! varying from ! s ' 0:883 to ! e ' 1:104
for k = 0 and n = 0; 1; 2: The green and purple dots are connecting pointsof the red and blue segments, the red segments start at the green points andend at the purple points while the blue segments end at the green points andstart at the purple points.6 It is to be noticed that the red and blue segments
shift upward when n increases and rightward when k increases. The L1 (0; 0)segment is located below the horizontal axis so it is not depicted. In order tokeep 2 non-negative, the lower most blue segment of L2 (0; 0) is illustrated only
0
7
for 1
1 ' 1:617: Increasing the value of 2 along the vertical dotted line0
at 1 = 1 intersects the partition curve three times ata2 ' 2:618;b
2 ' 6:837 andc
2 ' 8:697

which are denoted by the black dots.
Concerning stability, we …rst examine dynamics in the following two regionsof Figure 2, the LHS (left hand side) region shaded by positive-sloping lines in
6 In the case of n = 1; segment L (0; 1) starts at the end point (i.e., the lowest green point)1of L2 (0; 0) and ends at the starting point (i.e., the lower purple point) of L2 (0; 1):The sameconnection is repeated as k increases.
7 Taking k = n = 0 and solving + (!; 0) = 0 yields ! = !: Substituting into2
1 (!; 0)presents approximately this threshold value.

11mwhich 00 and the RHS (right hand side) region shaded
1 < 1 and 2
M
by the negative-sloping lines in which 1
0. Here m
1 and 2
1 is theMminimum 1 -value of the segments L2 (0; n) and 1 is the maximum 1 -valueof the segments L1 (0; n):8 Since the origin, 1 = 2 = 0 is in the LHS region,
any combination of the delays in the LHS region do not a¤ect stability of the
non-zero steady states and such a delay is called harmless. In order to arriveat the RHS region for any value of 2 , increasing 1 from zero must cross thestability switching curve and then the vertical dotted line at 1 = M
1 : This
implies that at least one of the eigenvalues must have a positive real part forany pair ( 1 ; 2 ) in the RHS region. In other words, the stability is lost in the
RHS region. Roughly speaking, for any value of 2 ; the stability of the non-zerosteady states is preserved for smaller values of 1 and lost for larger values of
1 . These results are summarized as follows:

Proposition 5 The non-zero steady states of system (7) are locally asymptotm
M
ically stable if 0
0 and always unstable if 1
1 < 1 and 2
1 and
0:
2

Figure 2. The stability switching curvefor k = 0 and n = 0; 1; 2.
M
We now turn attention to the remaining area in which m
1
1
1 and
0: The red-blue connected curve in Figure 2 is called a stability switching2
curve since dynamics is switched from stability to instability (i.e., stability loss)
or from instability to stability (i.e., stability gain) on this curve. We numerically
8 Solving d

into
way.

1 =d! = 0 for ! to obtain the minimizer ! m ' 1:084 and then substituting it
m ' 1:393: The other value M ' 1:813 is determined in the similar
11

1 (!; 0) yields

12con…rm how the stability switch takes place on the stability switching curve andwhat dynamics arises when a stability loss takes place.9 We …rst consider thee¤ects caused by the change in 1 with a …xed value of 2 and then proceed tothe e¤ect caused by the change in 2 ; with a …xed value of 1 .
Simulation 1: 1

1

3:5;

2 = 2:The dotted horizontal line at 1 = 2 crosses the red segment L1 (0; 1) atthe dotted point in Figure 2. Although it is not labelled, the 1 -value of theintersection will be denoted by B
1 ' 1:776 in which B implies bifurcation.
Figure 3 illustrates two bifurcation diagrams with respect to 1 in which delaysystem (7) runs for 0 t 1000 and initial functions de…ned for t 2 [ 1 ; 0] areselected to be constant, x0 (t) = x1 + 0:1 and y0 (t) = y1 + 0:2. In Figure 3(A),1 is increased from 1 to 3:5 with an increment of 1=400 along the dotted line at
t 1000 are
2 = 2. The local maximum and minimum values of y(t) for 950
plotted again each value of 1 to take away initial disturbances. The horizontal
Bpart of the diagram for y = y1 and 1
1 < 1 implies that the positive steadystate is locally asymptotically stable. The stability is lost at 1 = B
1 and alimit cycle oscillating around y1 emerges for 1 > B
whenthediagramhas two1
branches. It then gets complicated more through a period-doubling-like cascadeas 1 further increases. As 1 becomes closer to 3:5, the complicated dynamics
converges to a big limit cycle including two steady states inside. If we takethe values of the constant initial functions in the neighborhood of the negative
steady state, y2 ; we then have exactly the same dynamics whose bifurcation
diagram is symmetric to the one in Figure 3(A) with respect to the horizontal
axis, which is not illustrated to avoid confusion.
Simulation 2: 1

1

4:5;

2 =

a
2:

In Figure 3(B), the …xed value of 2 is increased to a2 and the same procedure
generates a di¤erent shape of the bifurcation diagram showing larger and more0complicated oscillations when stability is lost at 1 = B
1 (= 1 ) at which theahorizontal line at 2 = 2 crosses the L1 (0; 1) segment in Figure 2. It is to be0
Bnoticed that B
1 = 1 ' 1:617 in Figure 3(A) and 1 ' 1:776 in Figure 3(B)9 It is possible to analytically con…rm the direction of the stabililty switch. See Theorems
3 and 4 in Matsumoto and Szidarovszky (2015) and Proposition 6.1 in Gu et al. (2005).

13

though the same notation is used.

2 =

a
2

Figure 3. Bifurcation diagrams with respect to

1(A)

Simulation 3:

2 =2

1 = 2:5 and

(B)

2 =2

In the next simulation, we keep 2 = 2 and …x the values of 1 at 2:5. By
doing so, we will see dynamics scattered along the dotted vertical line at 1 = 2:5
seen in Figure 3(A) from a di¤erent point of view. With these parameter values,
the dynamic system runs for 0
t1000 with two di¤erent initial functions
de…ned for t 0, one with x0 (t) = x1 + 0:1 and y0 (t) = y1 + 0:2 as before andthe other with x0 (t) = x2 0:1 and y0 (t) = y2 0:2: Two phase diagrams areplotted in the (x; y) plane in Figure 4(A) in which a trajectory with the positive
initial functions converges to an upper attractor surrounding the positive steadystate y1 and a trajectory with the negative initial functions approaches thelower attractor surrounding the lower steady state y2 . The attractors are point
symmetric with respect to the origin. The corresponding time trajectories ofy(t) for 900 t 1000 are depicted in Figure 4(B) in which both are symmetricwith respect to the horizontal axis.
These simulations explain well our two typical love experiences. The …rst
one is that love developments strongly depends on the …rst feelings that theindividuals obtain at their …rst meeting. When trajectories start in the …rst
quadrant of the (x; y) plane, both individuals are thought to have "good feeling"
and can develop their romantic feelings to more favorable situation. On thecontrary when both have "not good feelings," the developments grow into the
opposite direction. The second one is that love feelings are often ‡uctuating ina wavering manner around their "true love" levels even if they fell in love. As
is seen in Figure 4(B), the romantic feeling of Rome and Juliet are alternating

14

between staying and going.

(A) Phase diagrams(B) Time trajectories

Figure 4. Dynamics with
Simulation 4:

1 = 3 and1 = 2:5 and

2 =2

2 =2

We increase 1 to 3 from 2:5; taking 2 = 2 in the next example. We again
observe love dynamics on the vertical line at 1 = 3 in Figure 3(A). Resultant
dynamics with the positive initial functions is depicted in red and dynamics
with the negative initial functions in black.10 A parametric di¤erence seems tobe minor, however, major di¤erent dynamics emerges. Comparing Figure 5(A)with Figure 4(A), we can see that increasing the value of 1 makes two attractors
merge into one large attractor including the two steady states inside by enlarging
the size of each attractor. In Figure 5(A) two large attractors are depicted in
red and black. Due to the point symmetry, the attractors are not identical,
however, generated dynamics are essentially the same and line-symmetric withrespect to the horizontal axis as seen in Figure 5(B). Although love and hate areopposite emotions, our daily-life experiences indicate that one can often turninto the other and vice versa. These examples describe these alternative feelings
between love and hate. These results remind us of cyclical emotional feelings ofPetrach, an Italian poet of the 14th century to Laura, a beautiful but married
woman described by Jones (1995). Stimulated by this work, Rinaldi (1998) has
successfully described the poet’s regular cyclic behavior ranging from ecstasy todespair with a dynamic system of three ordinary di¤erential equations without
delays but with instincts (i.e., appeals) which are assumed away in our model
1 0 The initial functions x (t) = 0:1 and y (t) = 0:2 for t < 0 are selected for the red curve00
and x0 (t) = 0:1 and y0 (t) = 0:2 for the black curve.

15

by Assumption 3.

(A) Phase diagrarams(B) Time trajectories

Figure 5. Dynamics with
Simulation 5:

1 = 3 and2 =1 = 3 and

2 =2

a
1

We now turn attention to the bifurcation diagram in Figure 3(B) with= a2 : Taking 1 = 3; Figure 6(A) depicts a phase diagram with the initial functions x0 (t) = 0:4 and y0 (t) = y2 0:1 for t 0: Notice that dynamics
starts with the point (x0 (0); y0 (0)) denoted by the black dot in the fourth quadrant of the (x; y) plane, oscillates around the negative steady state for a while
and then converges to a limit cycle surrounding the positive steady state. This
numerical example provides a green light to the Beast, a witched young handsome prince, in "Beauty and the Beast." Even starting with a unfavorable point
where the Beast has good feelings to Beauty named Belle but Beauty has negative feelings due to the terrifying and ugly appearance of the Beast, their lovestory described by delay system (7) with suitable parametric values arriving at
a regime where Belle is happy and the Beast can be transformed to a prince.11
It is also numerically con…rmed that a similar cyclic oscillation around the negative steady state is obtained when the initial point is selected in the second
quadrant where Romeo has negative feelings and Juliet has positive feeling when2

1 1 Rinaldi et al.(2013b) interprete this love story with non-delay model (1) without Assumption 3. Positivity of the appeal Az is a key factor that brings Beauty and the Beast to theirplateau.

16

they …rst meet.

(A) Phase diagram(B) Time trajectory

Figure 6. Dynamics with
Simulation 6:

1 = 3:5 and2 =

1 = 3 and

2 =

a
2

a
2

In the last simulation with Figure 3(B), the value of 1 is increased to 3:5.
More realistic love-hate alternations are described in Figures 7(A) and 7(B).
Emotions oscillate around the positive steady state (i.e., love regime) for awhile and then suddenly change their directions, going toward the negative
steady state, around which emotions oscillate in the neighborhood sooner orlatter, a sudden change occurs again to bring back the negative oscillations to
the positive oscillations. This alternative oscillation continues. Contrary to aregular cyclic pattern in Figure 5(A), a quite irregular pattern appears ranging
from one extreme to the other.

(A) Phase diagram(B) Time trajectory

Figure 7. Dynamics with17

1 = 3:5 and

2 =

a
2

Simulation 7: 1

1

2;

2 =

b
2:

Now the value of 2 is increased to b2 from a2 : In Figure 8(A) two blue
curves, one located slightly above the y = y1 locus and the other below they = y2 locus, present bifurcation diagrams with respect to 1 : With this, it is…rst observed that no more complicated dynamics appears but only a simple
limit cycle emerges when the stability of the non-zero steady states is lost forB
0 12This result implies two issues. One is that given 1 , a larger
1 > 1 (= 1 ):value of 2 generates simpli…ed dynamics. The other is that multi-stability canoccur as two branches are illustrated for a small deviation of 1 from B
1 : A basinof attraction is presented in Figure 8(B) with 1 = 3=2 and 2 = b2 : The regionof the initial points, (x0 (0); y0 (0)); is divided into three subregions, the light-red
region including the positive steady state denoted by the green dot, the lightblue region including the negative steady state by the red dot and the white
region including the zero-steady state denoted by the yellow dot. This divisionindicates an initial point dependency of dynamics, namely, dynamics converges
to the positive steady state, to the negative steady state or to a limit cycle
according to whether the point (x0 (0); y0 (0)) is selected from the red region, theblue region or the white region of Figure 8(B). This dependency also explains
something that happens often in a love story, that is, even though Romeo andJuliet are into something goods for example, a small misunderstanding describedby a small change in the initial functions can lead them to a di¤erent worlds inwhich one’s love feeling ‡uctuates from one extrem to the other.

(A) Bifurcation diagrams(B) Basin of attraction

Figure 8. Dynamics with
Simulation 8: 1
to

1

2;

2 =

1 = 3=2 and

2 =

a
2

c
2:

In the last simulation in this section, the value of 2 is further increasedc
2 . It is con…rmed that multi-stability occurs in the neighborhood of the

1 2 Although it is not illustrated in Figure 8(A), the simple shape of the bifurcation diagrams
is numerically con…rmed for larger values of 1 :18

0bifurcation value of B
1 (= 1 ): As can be seen, small red limit cycles surrounding
the non-zero steady state coexist with a large blue limit cycle surrounding thezero steady state and including the non-zero steady states inside. Concerning
the red cycle, we notice the followings: the positive steady state y1 is stable forBB
1 < 1 and bifurcates to a limit cycle for 1 > 1 ; a radius of the limit cycle
is sensitive to the value of 1 ; the red cycle suddenly disappears when the valueof 1 is larger than about roughly 1:9. Concerning the blue cycle, we also notice
the followings: the limit cycle can emerge for 1 < B
1 ; for 1 in the vicinityof 1:8, the bifurcation diagram has a collection of the local maximum value,indicating that the cycle repeats small ‡uctuations in the local maximum andminimum, although ‡udtuations in the local minimum are not seen in Figure
9(A). These numerical simulations imply the inital point dependency with whicha small change in the initial point might result in largely di¤erent dynamics.

(A) Bifurcation diagrams(B) Time trajectories

Figure 9. Dynamics with

4.2

x

y >

1 = 1:8 andc2x y > 0

We now focus on the case in which the self reaction dominates the cross reaction.
Keeping the values of x and y …xed, we increase the values of x and y to
the extent that the direction of the inequality is reversed:
Assumption 5:x =y = 8=5.

Under these parametric conditions, two issues can be con…rmed: (1) the basicmodel has the zero solution as the unique steady state and (2) the eigenvaluesof the characteristic equation (4) are real and negative indicating that the zerosteady state is a stable node, as shown in Figure 1(B). Our concerns are, as
before, upon the delay e¤ects on dynamics, in particular we are interested inwhether larger delays can destabilize the zero steady state and what dynamics
could emerge when the stability is lost.
19

Figure 10 illustrates the parts of the segments of L1 (k; n) in red and L2 (k; n)in blue for k = 0 and n = 0; 1; 2 under Assumption 5 and x = y = 3=2. The
stability switching curve is the outer envelope of the partition curves for various
values of n and its shape di¤ers from the one in Figure 2.13 Since the steadystate is stable without delays (i.e., 1 = 2 = 0) and the origin of Figure 10
is located left to the stability switching curve, the steady state with positive
delays is stable in the left part to the curve and unstable in the right. The valueBof B
1 (= 0:6) is selected in such a way that the dotted vertical line at 1 = 1
intersects the stability switching curve three times ata2 ' 1:771;b
2 ' 3:536 andc
2 ' 4:478:

Multiple intersections suggest that delay 2 has the double edge e¤ect forB
1 , which is numerically con…rmed below.

1 =Figure 10. Stability switching curves
Simulation 9:

1 =

B
1 ;

0

2

6:

Figure 11(A) presents a bifurcation diagram with respect to 2 in which avalue of 2 is increased along the dotted vertical line at 1 = B
1 in Figure 10.
As is seen both in Figure 10 and Figure 11(A), the steady state loses stability
at 2 = a2 and 2 = c2 and regains it at 2 = b2 : Further, Figure 11(A) depicts
time trajectories of x(t) in blue and y(t) in red. It is seen that the steadystate bifurcates to a limit cycle surrounding the zero steady state whenever thestability is lost and the cycle merges to the zero steady state when the stability
is regained. Provided a positive delay in 1 ; some positive value of 2 keeps both
1 3 Under these speci…cations, the end points of the domain for the stability switching curves
are ! s ' 0:805 and ! e ' 1:163:

20

Romeo and Juliet in suspense resulting in a cyclic behavior and some other value
makes them stay at the zero steady state.
Simulation 10: 0:5

1

1;

2 =

a
2

Figure 11(B) gives bifurcation diagrams of x(t) in blue and y(t) in red withrespect to 1 : Two issues are observed: (i) Given the …xed value of 2 ; di¤erent
values of 1 do not qualitatively a¤ect dynamics as the time trajectories in
Figure 11(B) exhibit similar behavior (i.e., limit cycles) and (ii) time trajectoriesof x(t) exhibit huge oscillations for 1 -values close to unity, indicating that themodel may be inappropriate to describe love evolutions of Romeo and Juliet.

(A) With respect to(B) With respect to2

1

Figure 11. Bifurcation diagrams

5

Concluding Remarks

We have investigated a delay Romeo and Juliet model in which both are supposed to be cautious lovers, the most natural romantic style between two individuals. After checking that the no-delay version of the model exhibits fairly
stable dynamics, we theoretically derive the instability conditions of the delay
version by analyzing the characteristic equations and then use them to illustrate a stability switching curve on which stability of the delay version is lost.
Concerning global behavior, we perform numerical simulations to examine thatthe unstable steady state turns to cyclic oscillations through Hopf bifurcation.
In comparison with a one-delay model, the stability switching curve has more
complicated shape and dynamics become more complicated in the two-delay
model. In consequence, it is con…rmed that the multiple delays have the double edge e¤ect, implying the occurrence of multiple alternations of stability andinstability with increasing value of one delay and …xing the value of the otherdelay at some level. Furthermore, various love evolutions ranging from simple

21

to complex dynamics involving chaos can explain various types of love stories
between two individuals.
For further investigation, we …rst plan to introduce delays in the two selfreaction processes, with which the delay model could be more realistic. Secondly,
it may be interesting to add one more individual to the model and consider atriangle relation whose dynamics could be much more complicated and therefore
much more interesting.

22References
Bielczyk, N., U. Foryś and T. Platkowski, Dynamical Models of Dyadic Interactions with Delay, Journal of Mathematical Sociology 37, 223-249, 2013.
Gu, K., S-I., Nicolescu and J. Chen, On Stability Crossing Curves for General
Systems with Two Delays, Journal of Mathematical Analysis and Applications 311, 231-253, 2005.
Liao, X. and J. Ran, Hopf Bifurcation in Love Dynamical Models with Nonlinear Couples and Time Delays, Chaos Solitions and Fractals 31, 853-865,
2007.
Jones, F., J. The Structure of Petrarchi’s Canzoniere, Brewer, Cambridge, UK,
1995.
Matsumoto, A. and F. Szidarovszky, Love A¤airs Dynamics with One Delayin Losing Memory or Gaining A¤ection, IERCU DP #260, Institute of
Economic Research, Chuo University, 2016.
Matsumoto, A. and F. Szidarovszky, Learning Monopolies with Delay Feedbac on Price Expectations, Communications in Nonlinear Science andNumerical Simulation 28, 151-165, 2015.
Rinaldi, S., F. D. Rossa and P. Landi, A Mathematical Model of "Pride andPrejudice," Nonlinear Dynamics, Psychology, and Life Science 18, 199211, 2014.
Rinaldi, S., F. D. Rossa and P. Landi, A Mathematical Model of "Gone with
the Wind," Physica A 392, 3231-3239, 2013a.
Rinaldi, S., P. Landi and F. D. Rossa, Small Discoveries Can Have Great Consequences in Love A¤airs: The Case of Beauty and Beast, International
Journal of Bifurcation and Chaos 23, doi: 10.1142/S0218127413300383,
2013b.
Rinaldi, S. and A. Gragnani, Love Dynamics Between Secure Individuals: A
Modeling Approach, Nonlinear Dynamics, Psychology and Life Science 2,283-301, 1998.
Rinaldi, S., Laura and Petrarch: an Intriguing Case of Cyclical Love Dynamics,
SIAM Journal of Applied Mathematics 58, 1205-1221, 1998.
Son, W-S. and Y-J. Park, Time Delay E¤ect on the Love Dynamical Model,
Journal of the Korean Physical Society 59, 2197-2204, 2011.
Strogatz, S., Love A¤airs and Di¤erential Equations, Mathematics Magazine
61, p. 35, 1988.

23