︠3d5a2adc-356e-4b9b-8a4e-bd313ac565f2si︠ %hide %html

A seriously silly model of romantic love!

How do I love thee? Let me model the ways!
- My apologies to Elizabeth Barrett Browning for nerdifying her lovely poem!

(Adapted from various sources. But if you're seriously interested, check out this article!)

The scene

Juliet loves Romeo, and Romeo loves Juliet, at least at time t=0. The question is, how will their feelings evolve with time, and as they meet other, potentially fascinating, new people! Here is our modeling strategy:

Love is an ever changing function

 
To help you get comfortable with these love functions and their phase plots, let us begin with a warmup exercise.

Each graph below shows the time history of various hypothetical love affairs between Juliet and Romeo. Interpret what is happening in each graph.
Fig. 1 Fig. 2
Fig. 3 Fig. 4
Building a simple love-affair model
 
For simplicity, we will assume that there are only two factors that affect Juliet's future feelings for Romeo: (i) her own current feelings about him, and (ii) his current feelings towards her. And similarly, there are only two factors that affect Romeo's future feelings about Juliet.

As usual, we want to think in terms of the derivatives $\displaystyle \frac{dJ}{dt}$ and $\displaystyle \frac{dR}{dt}$, in order to model $J(t)$ and $R(t)$. Here are some question to explore:

A general formulation of the model that includes all the cases considered above is
                      $J^\prime(t) = a J + b R $
                      $R^\prime(t) = c J + d R $

where all the coefficients denote constant parameter values that model different personality types and/or behavior traits. For example, here is one classification scheme based on Romeo's personality that has been proposed in the literature:

  1. Eager beaver: $ c > 0$ and $ d > 0$ (Romeo is encouraged by his own feelings as well as by Juliet's.)
  2. Narcissistic nerd: $ c < 0$ and $ d > 0$ (Romeo wants more of what he feels, but retreats from Juliet's feelings.)
  3. Cautious lover: $ c > 0$ and $ d < 0$ (Romeo retreats from his own feelings, but is encouraged by Juliet's.)
  4. Hermit: $ c < 0$ and $ d < 0$ (Romeo retreats from his own feelings as well as from Juliet's.)
︡21de48c9-b601-46f2-9488-84fbeb8ecc91︡{"hide":"input"}︡{"html":"\n

\nA seriously silly model of romantic love!\n

\n\n\n\n \n \n \n
\n \n How do I love thee? Let me model the ways!\n
\n - My apologies to Elizabeth Barrett Browning\n for nerdifying her lovely poem!\n \n
\n\n

\n(Adapted from various sources. But if you're \nseriously interested, check out\nthis article!)\n\n\n\n\n

\n\n\nThe scene\n\n

\nJuliet loves Romeo, and Romeo loves Juliet, \nat least at time t=0. The question is, how will their \nfeelings evolve with time, and as they meet other, potentially \nfascinating, new people! Here is our modeling strategy:\n\n\n

\n\n\nLove is an ever changing function\n\n

 
\nTo help you get comfortable with these love functions and their \nphase plots, let us begin with a warmup exercise.\n

\n\nEach graph below shows the time history of various hypothetical \nlove affairs between Juliet and Romeo. Interpret what is \nhappening in each graph.\n\n\n \n \n \n \n \n \n \n \n\n\n \n \n \n \n \n \n \n \n
\n \n \n \n
\n Fig. 1\n \n Fig. 2\n
\n \n \n \n
\n Fig. 3\n \n Fig. 4\n
\n \n \n\nBuilding a simple love-affair model\n\n
 
\nFor simplicity, we will assume that there are only two \nfactors that affect Juliet's future feelings for Romeo: \n(i) her own current feelings about him, and (ii) his \ncurrent feelings towards her. And similarly, there\nare only two factors that affect Romeo's future\nfeelings about Juliet.

\n\nAs usual, we want to think in terms of the derivatives \n$\\displaystyle \\frac{dJ}{dt}$ and $\\displaystyle \\frac{dR}{dt}$, \nin order to model $J(t)$ and $R(t)$. Here are some question \nto explore:\n

\n \nA general formulation of the model \nthat includes all the cases considered above is
\n\n                 \n    $J^\\prime(t) = a J + b R $
\n                 \n    $R^\\prime(t) = c J + d R $

\nwhere all the coefficients denote constant parameter \nvalues that model different personality types and/or \nbehavior traits. For example, here is one classification \nscheme based on Romeo's personality that has been proposed \nin the literature: \n\n

    \n
  1. Eager beaver: $ c > 0$ and $ d > 0$ \n(Romeo is encouraged by his own feelings as well as by Juliet's.)\n
  2. Narcissistic nerd: $ c < 0$ and $ d > 0$ \n(Romeo wants more of what he feels, but retreats from Juliet's \nfeelings.)\n
  3. Cautious lover: $ c > 0$ and $ d < 0$ \n(Romeo retreats from his own feelings, but is encouraged by\nJuliet's.)\n
  4. Hermit: $ c < 0$ and $ d < 0$ \n(Romeo retreats from his own feelings as well as from \nJuliet's.)\n
\n"}︡{"done":true}︡ ︠82feebf6-5e07-40d9-9066-ae73b30ebe10si︠ %hide %html Explore the model's predictions
 
The four parameters $a$, $b$, $c$, $d$, in the general formulation above can be chosen to mimic different personality types. For example, an affair between two Eager beavers would correspond to positive values for all 4 parameters. Pick a few different personality combinations and explore what the model predicts about their love affair. Be sure to try different starting points for each situation; e.g., $J(0)$ , $R(0)$ need not always be positive. Use the interactive Sage script below to carry out your explorations. ︡48a187ce-5d08-4331-b311-b1a8f8451287︡{"hide":"input"}︡{"html":"\n\nExplore the model's predictions\n\n
 
\n\nThe four parameters $a$, $b$, $c$, $d$, in the general formulation \nabove can be chosen to mimic different personality types. For example, \nan affair between two Eager beavers would correspond to positive \nvalues for all 4 parameters. Pick a few different personality \ncombinations and explore what the model predicts about their \nlove affair. Be sure to try different starting points for each \nsituation; e.g., $J(0)$ , $R(0)$ need not always be positive. Use \nthe interactive Sage script below to carry out your explorations.\n\n"}︡{"done":true}︡ ︠a2987edf-49f6-4f58-9500-7dfda967dbb2si︠ %hide # System is: j' = a*j + b*r ; r' = c*j + d*r jl, rl, t = var('jl rl t') print "Interactive simulation of simple love-affair model." print "Enter/adjust parameter values for system: dJ/dt = a*J + b*R ; dR/dt = c*J + d*R" @interact def _( a=slider(-5, 5, step_size=0.5, default=1, label='a', width=20), b=slider(-5, 5, step_size=0.5, default=1, label='b', width=20), c=slider(-5, 5, step_size=0.5, default=1, label='c', width=20), d=slider(-5, 5, step_size=0.5, default=1, label='d', width=20), ic_j=input_box(0.5, 'J(0)=', width=20), ic_r=input_box(1, 'R(0)=', width=20), t_end=input_box(5, 'End time=', width=20) ): de1 = a*jl + b*rl de2 = c*jl + d*rl P = desolve_system_rk4 ([de1, de2], [jl, rl], ics=[0, ic_j, ic_r], ivar=t, end_points=[0,t_end] ) Q = [ [j,k] for i, j, k in P] #list_plot(Q, axes_labels=['$S$', '$B$']) show(line(Q, axes_labels=['$J$', '$R$'], color='green', thickness=2, title='Phase plot')) # The following creates plots of S vs t and B vs t: QS = [ [i,j] for i, j, k in P] QB = [ [i,k] for i, j, k in P] SS = line(QS, axes_labels=['$t$', 'J/R'], legend_label='J', color='blue', thickness=2) SB = line(QB, axes_labels=['$t$', 'J/R'], legend_label='R', color='green', thickness=2) show(SS+SB) ︡825cb5ca-6743-4dc2-bfaa-2e16267e4054︡{"hide":"input"}︡{"stdout":"Interactive simulation of simple love-affair model.\n"}︡{"stdout":"Enter/adjust parameter values for system: dJ/dt = a*J + b*R ; dR/dt = c*J + d*R\n"}︡{"interact":{"controls":[{"animate":true,"control_type":"slider","default":12,"display_value":true,"label":"a","vals":["-5.00000000000000","-4.50000000000000","-4.00000000000000","-3.50000000000000","-3.00000000000000","-2.50000000000000","-2.00000000000000","-1.50000000000000","-1.00000000000000","-0.500000000000000","0.000000000000000","0.500000000000000","1.00000000000000","1.50000000000000","2.00000000000000","2.50000000000000","3.00000000000000","3.50000000000000","4.00000000000000","4.50000000000000","5.00000000000000"],"var":"a","width":20},{"animate":true,"control_type":"slider","default":12,"display_value":true,"label":"b","vals":["-5.00000000000000","-4.50000000000000","-4.00000000000000","-3.50000000000000","-3.00000000000000","-2.50000000000000","-2.00000000000000","-1.50000000000000","-1.00000000000000","-0.500000000000000","0.000000000000000","0.500000000000000","1.00000000000000","1.50000000000000","2.00000000000000","2.50000000000000","3.00000000000000","3.50000000000000","4.00000000000000","4.50000000000000","5.00000000000000"],"var":"b","width":20},{"animate":true,"control_type":"slider","default":12,"display_value":true,"label":"c","vals":["-5.00000000000000","-4.50000000000000","-4.00000000000000","-3.50000000000000","-3.00000000000000","-2.50000000000000","-2.00000000000000","-1.50000000000000","-1.00000000000000","-0.500000000000000","0.000000000000000","0.500000000000000","1.00000000000000","1.50000000000000","2.00000000000000","2.50000000000000","3.00000000000000","3.50000000000000","4.00000000000000","4.50000000000000","5.00000000000000"],"var":"c","width":20},{"animate":true,"control_type":"slider","default":12,"display_value":true,"label":"d","vals":["-5.00000000000000","-4.50000000000000","-4.00000000000000","-3.50000000000000","-3.00000000000000","-2.50000000000000","-2.00000000000000","-1.50000000000000","-1.00000000000000","-0.500000000000000","0.000000000000000","0.500000000000000","1.00000000000000","1.50000000000000","2.00000000000000","2.50000000000000","3.00000000000000","3.50000000000000","4.00000000000000","4.50000000000000","5.00000000000000"],"var":"d","width":20},{"control_type":"input-box","default":"0.500000000000000","label":"J(0)=","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"ic_j","width":20},{"control_type":"input-box","default":1,"label":"R(0)=","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"ic_r","width":20},{"control_type":"input-box","default":5,"label":"End time=","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"t_end","width":20}],"flicker":false,"id":"aa0fed41-ce6a-4a91-8300-87d15ef6e4f0","layout":[[["a",12,null]],[["b",12,null]],[["c",12,null]],[["d",12,null]],[["ic_j",12,null]],[["ic_r",12,null]],[["t_end",12,null]],[["",12,null]]],"style":"None"}}︡{"done":true}︡ ︠43f881e8-957e-43a8-b816-35a210dd4cc9si︠ %hide %html A contrived love affair
 
It is not too hard to design a model to mimic any long-term behavior we want. For example, suppose we want both lovers to stay in love permanently, after an initial period of meandering and uncertainty! Consider the following model:


                      $J^\prime(t) = -0.2(J-3) - (R-4) $
                      $R^\prime(t) = (J-3) - 0.2(R-4) $


Pick any choice of initial values for $J$ and $R$ -- in fact, try various combinations such as hate-hate, love-hate, etc. Use the Sage script below to run these simulations and discuss your results: what are the values of $J(t)$ and $R(t)$ in the long-term? arguing from a calculus perspective, can you conjecture why? what are the personality types of Romeo and Juliet in this model (e.g., Eager beaver, Narcissistic nerd, etc.)? Can you think of a simple change to the model that would make them both hate each other in the long-term? ︡eff074c4-e9df-41eb-87c0-9e02ea618655︡{"hide":"input"}︡{"html":"\n\nA contrived love affair\n\n
 
\nIt is not too hard to design a model to mimic any \nlong-term behavior we want. For example, suppose we want both \nlovers to stay in love permanently, after an initial period of \nmeandering and uncertainty! Consider the following model:


\n                 \n    $J^\\prime(t) = -0.2(J-3) - (R-4) $
\n                 \n    $R^\\prime(t) = (J-3) - 0.2(R-4) $


\n\n\nPick any choice of initial values for $J$ and $R$ -- in fact, \ntry various combinations such as hate-hate, love-hate, etc. \nUse the Sage script below to run these simulations and discuss \nyour results: what are the values of $J(t)$ and $R(t)$ \nin the long-term? arguing from a calculus perspective, \ncan you conjecture why? what are the personality types of \nRomeo and Juliet in this model (e.g., Eager beaver, \nNarcissistic nerd, etc.)? Can you think of a simple \nchange to the model that would make them both hate each\nother in the long-term?\n"}︡{"done":true}︡ ︠d46f042e-e2a1-40ea-b240-7f237c270ed6si︠ %hide # System is: j' = -0.2*(j-3) - (r-4) ; r' = (j-3) - 0.2(r-4) j, r, t = var('j r t') print "Interactive simulation of model: j' = -0.2*(j-3) - (r-4) , r' = (j-3) - 0.2(r-4)" print "Enter initial values for J and R" @interact def _( ic_j=input_box(0.5, '$J(0)=$', width=20), ic_r=input_box(1, '$R(0)=$', width=20), t_end=input_box(10, 'End time=', width=20) ): de1 = -0.2*(j-3) - (r-4) de2 = (j-3) - 0.2*(r-4) P = desolve_system_rk4 ([de1, de2], [j, r], ics=[0, ic_j, ic_r], ivar=t, end_points=[0,t_end] ) Q = [ [jc,k] for i, jc, k in P] #list_plot(Q, axes_labels=['$S$', '$B$']) show(line(Q, axes_labels=['$J$', '$R$'], color='green', thickness=2, title='Phase plot')) # The following creates plots of S vs t and B vs t: QS = [ [i,jc] for i, jc, k in P] QB = [ [i,k] for i, jc, k in P] SS = line(QS, axes_labels=['$t$', 'J/R'], legend_label='J', color='blue', thickness=2) SB = line(QB, axes_labels=['$t$', 'J/R'], legend_label='R', color='green', thickness=2) show(SS+SB) ︡b0c9d035-2f9a-42af-8312-65d206b0d179︡{"hide":"input"}︡{"stdout":"Interactive simulation of model: j' = -0.2*(j-3) - (r-4) , r' = (j-3) - 0.2(r-4)\n"}︡{"stdout":"Enter initial values for J and R\n"}︡{"interact":{"controls":[{"control_type":"input-box","default":"0.500000000000000","label":"$J(0)=$","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"ic_j","width":20},{"control_type":"input-box","default":1,"label":"$R(0)=$","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"ic_r","width":20},{"control_type":"input-box","default":10,"label":"End time=","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"t_end","width":20}],"flicker":false,"id":"1eac2357-2586-4fdd-a0d4-e2afc76d6936","layout":[[["ic_j",12,null]],[["ic_r",12,null]],[["t_end",12,null]],[["",12,null]]],"style":"None"}}︡{"done":true}︡ ︠d4515c4a-e8f8-432b-8948-c75a1dd077a8si︠ %hide %html The mathematics behind it
 
First, note the fixed points of the system:


                      $J^\prime(t) = -0.2(J-3) - (R-4) $
                      $R^\prime(t) = (J-3) - 0.2(R-4) $

Next, let's look at the Jacobian matrix and is eigenvalues. Since the system is linear, the terms in the Jacobian (say, $\bf A$) are just the coefficients: $ {\bf A} = \left[\begin{array}{cc} -0.2 & -1 \\ 1 & -0.2 \end{array}\right]$

Its eigenvalues are found using the Sage segment below. ︡54d766a1-8406-4f74-9d20-ec7eddb7762d︡{"hide":"input"}︡{"html":"\n\nThe mathematics behind it\n\n
 
\nFirst, note the fixed points of the system:


\n\n                 \n    $J^\\prime(t) = -0.2(J-3) - (R-4) $
\n                 \n    $R^\\prime(t) = (J-3) - 0.2(R-4) $

\nNext, let's look at the Jacobian matrix and is eigenvalues. \nSince the system is linear, the terms in the Jacobian (say, $\\bf A$) are \njust the coefficients:\n$ {\\bf A} = \\left[\\begin{array}{cc} -0.2 & -1 \\\\ \n 1 & -0.2 \\end{array}\\right]$

\n\nIts eigenvalues are found using the Sage segment below.\n\n"}︡{"done":true}︡ ︠2b4348b7-1457-41f6-90ba-d6cb57214a4ds︠ # Define Jacobian matrix by rows in brackets [] separated by commas A = matrix( [ [-0.2, -1], [1, -0.2] ] ) # Find eigenvalues of matrix A.eigenvalues() ︡6014a66d-d461-4cd5-b0e9-2947d7f46bd7︡{"stdout":"[-0.200000000000000 - 1.00000000000000*I, -0.200000000000000 + 1.00000000000000*I]\n"}︡{"stderr":":1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.\n"}︡{"done":true}︡ ︠f08821ec-cddb-479b-884e-f5701f3eaf11is︠ %hide %html Notice the eigenvalues are complex conjugates, with negative real part. These types of eigenvalues always attract the solution to the fixed point in an oscillating pattern. Hence, the phase plot displays converging spiral solutions.


An easy way to change the problem and make the solution diverge would be to change the matrix so its eigenvalues are positive. For example, $ {\bf A} = \left[\begin{array}{cc} 0.2 & -1 \\ 1 & 0.2 \end{array}\right]$ has complex conjugate eigenvalues with positive real part. Thus, the following system will have only diverging spiral solutions

                      $J^\prime(t) = 0.2(J-3) - (R-4) $
                      $R^\prime(t) = (J-3) + 0.2(R-4) $ ︡96c5aa69-48d2-4821-97d3-668b828cf62a︡{"hide":"input"}︡{"html":"\nNotice the eigenvalues are complex conjugates, with negative \nreal part. These types of eigenvalues always attract the \nsolution to the fixed point in an oscillating pattern. Hence, \nthe phase plot displays converging spiral solutions.


\n\nAn easy way to change the problem and make the solution \ndiverge would be to change the matrix so \nits eigenvalues are positive. For example, \n$ {\\bf A} = \\left[\\begin{array}{cc} 0.2 & -1 \\\\ \n 1 & 0.2 \\end{array}\\right]$ has complex conjugate \neigenvalues with positive real part. Thus, the following system \nwill have only diverging spiral solutions

\n\n                 \n    $J^\\prime(t) = 0.2(J-3) - (R-4) $
\n                 \n    $R^\\prime(t) = (J-3) + 0.2(R-4) $ "}︡{"done":true}︡ ︠f50a231a-f746-4167-9823-2adfb0d475dbsi︠ %hide %html

More complex love affairs
The models we have seen so far are all linear, and their solutions exhibit a fairly limited range of behaviors. As you might suspect, things get much more interesting when we try to model non-linear personalities, or other more complicated situations such as 3-way affairs or polyamorous relationships.

  1. Nonlinear personalities
    Here is one simple way to make the Juliet-Romeo model nonlinear
                          J'(t) = a J + b R(1-|R|)
                          R'(t) = c J(1-|J|) + d R

    Can you conjecture what type of real-life behavior the nonlinear terms might capture (assume b and c are positive)? Depending on the choice of parameters, different solution behaviors can be seen. For example, try to run some simulations with a=b=-2 and c=d=1. To really understand what is going on, you will need to try a few different initial values for J(0) and R(0).

  2. A 3-way affair
    Suppose Juliet has an affair with someone (say, Sky). A model will be much easier to develop if she keeps the affair a secret from Romeo! In that case there will be 2 equations for Juliet, plus 1 each for Romeo and Sky. Here is one example of such a model

                          JR'(t) = a JR + b (R-S)
                          JS'(t) = a JS + b (S-R)
                          R'(t) = c JR + d R
                          S'(t) = e JS + f S

    where JR(t) and JS(t) are functions that represent Juliet's love for Romeo and Sky, respectively. As before, R(t) and S(t) represent the respective partners' love for Juliet. If you assume positive parameter values for a-f, can you conjecture meaningful interpretations for all the terms in this model?

    Despite the seeming complexity of this last model, it is still linear, and has a fairly limited range of solution behaviors. To see chaotic solution behaviors it is necessary to include nonlinear terms in the model. ︡ba2a08de-57f2-4489-af62-7ce8d358a22e︡{"hide":"input"}︡{"html":"\n\n

    More complex love affairs
    \n\n\nThe models we have seen so far are all linear, and their \nsolutions exhibit a fairly limited range of behaviors. \nAs you might suspect, things get much more interesting \nwhen we try to model non-linear personalities, or other \nmore complicated situations such as \n3-way affairs or polyamorous relationships.

    \n\n

      \n
    1. Nonlinear personalities
      \nHere is one simple way to make the Juliet-Romeo \nmodel nonlinear \n
      \n\n                 \n    J'(t) = a J + b R(1-|R|)
      \n                 \n    R'(t) = c J(1-|J|) + d R

      \nCan you conjecture what type of real-life behavior \nthe nonlinear terms might capture (assume b \nand c are positive)? \nDepending on the choice of parameters, different solution \nbehaviors can be seen. For example, try to run some \nsimulations with a=b=-2 \nand c=d=1. To really understand what is going \non, you will need to try a few different initial values for \nJ(0) and R(0).\n

      \n\n

    2. A 3-way affair
      \nSuppose Juliet has an affair with someone (say, Sky). \nA model will be much easier to develop if she keeps \nthe affair a secret from Romeo! In that case there \nwill be 2 equations for Juliet, plus 1 each for Romeo and \nSky. Here is one example of such a model

      \n                 \n    JR'(t) = \n a JR + b (R-S)
      \n                 \n    JS'(t) = \n a JS + b (S-R)
      \n                 \n    R'(t) = c JR + d R
      \n                 \n    S'(t) = e JS + f S

      \n\nwhere JR(t) and \nJS(t) are functions \nthat represent Juliet's love for Romeo and \nSky, respectively. As before, R(t) and \nS(t) represent the respective partners' \nlove for Juliet. If you assume positive parameter values \nfor a-f, can you conjecture meaningful \ninterpretations for all the terms in this model?

      \n\nDespite the seeming complexity of this last model, it is \nstill linear, and has a fairly limited range of solution \nbehaviors. To see chaotic solution behaviors \nit is necessary to include nonlinear terms in the model."}︡{"done":true}︡