{"id":2811,"date":"2018-04-27T02:45:33","date_gmt":"2018-04-27T07:45:33","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=2811"},"modified":"2018-04-25T13:01:11","modified_gmt":"2018-04-25T18:01:11","slug":"of-love-music-shakespeare-and-dynamic-systems","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2018\/04\/27\/of-love-music-shakespeare-and-dynamic-systems\/","title":{"rendered":"Of Love, Music, Shakespeare and Dynamic Systems"},"content":{"rendered":"\r\n<div class=\"content\"><!--introduction--><p>For any human being love is one of the biggest source of joy, happiness...problems and puzzles. Today's guest blogger, Aldo Caraceto, one of my fellow Application Engineers from Italy, is going to convince you to approach it from a different angle.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#0ce8abba-1b14-4184-aa10-e5317781e5f0\">From Love to MATLAB<\/a><\/li><li><a href=\"#e4fcdebb-baa3-469c-ab0b-0c0660162cb0\">A very simple emotional mechanism<\/a><\/li><li><a href=\"#98f17dce-170c-4cc9-bea2-1bfc12784d3d\">When things get complicated...<\/a><\/li><li><a href=\"#9cd6c125-d452-4d95-b445-d7248e766a48\">Phase portrait and ODE solutions<\/a><\/li><li><a href=\"#4f31094f-ff5e-4e04-a5c2-6d880c558ea1\">Relationship Models<\/a><\/li><li><a href=\"#03450f15-bf0e-4a87-a588-7064121174a8\">Conclusion<\/a><\/li><\/ul><\/div><h4>From Love to MATLAB<a name=\"0ce8abba-1b14-4184-aa10-e5317781e5f0\"><\/a><\/h4><p>I love music. One of my favourite songs is the popular Toto track <i>\"Hold The Line\"<\/i>. Very briefly, this song is about love and timing, as you may infer from the main chorus of the song: \"Hold the line, love isn't always on time\". Ok, it's not always on time, but can I predict when it will be? To answer this question, I've found a couple of insightful works, <i>\"Strogatz, S.H. (1988), Love affairs and differential equations, Math. Magazine 61, 35\"<\/i> and <i>\"Eva M. Strawbridge, Romeo and Juliet, a Dynamical Love Affair\"<\/i> which I've translated in MATLAB code. The mentioned papers answer a very simple question: how to express mathematically love dynamics, taking inspiration from <i>\"Romeo &amp; Juliet\"<\/i> by W.Shakespeare. MATLAB may be a good friend to help us in examining the evolution of this relationship. An <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/ode23.html\"><tt>ode23<\/tt><\/a> solver has been used to solve the system of equations representing the story. The results are quite funny.<\/p><h4>A very simple emotional mechanism<a name=\"e4fcdebb-baa3-469c-ab0b-0c0660162cb0\"><\/a><\/h4><p>This first example examines the condition where each person&#8217;s feelings are only affected by the other. Besides, we want our Romeo not following the expected behaviour: in this case, he'll be a confused and confusing lover: the more Juliet loves him, the more he begins to dislike her, generating surprising reactions with regards to the original tragedy.<\/p><p>$\\frac{dr}{dt} = - a \\cdot j,$<\/p><p>$\\frac{dj}{dt} = b  \\cdot r$<\/p><p>In the system of differential equations representing the phenomenon, <tt>a<\/tt> and <tt>b<\/tt> are real, positive numbers, and measured in 1\/time, hence are frequencies: the bigger they are, the shorter in time will be the swing of emotions for the two lovers. In addition, just a quick summary of the variables represented in the plots:<\/p><div><ul><li><tt>R(t)<\/tt> = Romeo&#8217;s feelings for Juliet at time t .<\/li><li><tt>J(t)<\/tt> = Juliet&#8217;s feelings for Romeo at time t.<\/li><li><tt>R(t), J(t) &gt; 0<\/tt> signify love, passion, and attraction.<\/li><li><tt>R(t), J(t) &lt; 0<\/tt> values signify dislike.<\/li><li><tt>R(t),J(t) == 0<\/tt> signifies indifference.<\/li><\/ul><\/div><pre class=\"codeinput\">a = 2e-1;\r\nb = 5e-1;\r\nloveRJ = @(t,y) [-a*y(2); b*y(1)];\r\n<\/pre><p>Now, we can set the problem up, collect data, and plot results, adding some annotations:<\/p><pre class=\"codeinput\">tstart = 0;\r\ntfinal = 50;\r\ntspan = [tstart tfinal];\r\ny0 = [-1;2]; <span class=\"comment\">% initial conditions<\/span>\r\n[t,y] = ode23(@(t,y) loveRJ(t,y),tspan,y0);\r\nplot(t,y')\r\nax = gca;\r\nax.XLim = [0 50];\r\nax.YLim = [-3 3];\r\nax.XLabel.String = <span class=\"string\">'Time'<\/span>;\r\nax.YLabel.String = <span class=\"string\">'Emotions'<\/span>;\r\nax.Title.String = <span class=\"string\">'Romeo &amp; Juliet''s relationship'<\/span>;\r\nlegend(<span class=\"string\">'Romeo'<\/span>,<span class=\"string\">'Juliet'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2018\/RJ_1_1_LSnew_Final_01.png\" alt=\"\"> <h4>When things get complicated...<a name=\"98f17dce-170c-4cc9-bea2-1bfc12784d3d\"><\/a><\/h4><p>In our first example, we have described the story with the changing emotions felt by Romeo and Juliet as they feed on each other's passion and indifference while the time is passing by. As you know, things - and people - are never so simple. Let's try to describe the relationship between Romeo and Juliet a bit differently, trying to compare the first, simple representation of the relationship with a more complex version. Then, we'll plot them together.<\/p><p>The new mathematical model can be expressed as:<\/p><p>$\\frac{dr}{dt} = a \\cdot r + b \\cdot j ,$<\/p><p>$\\frac{dj}{dt} = c \\cdot r + d \\cdot j$<\/p><p>As you may notice, this time there are  four parameters, because we'd like to consider also how Romeo and Juliet are influenced by their own feelings:<\/p><pre class=\"codeinput\">a = -0.15;\r\nd = 0.17;\r\nb = 0.9;\r\nc = -0.9;\r\n<\/pre><p>We calculate the solutions for the two relationship model (you may find the model definitions at the end of the article):<\/p><pre class=\"codeinput\">tstart = 0;\r\ntfinal = 50;\r\ntspan = [tstart tfinal];\r\ny0 = [1.00;-0.5]; <span class=\"comment\">% initial conditions<\/span>\r\n\r\n[t2,y2] = ode23(@(t,y) loveRJ_simple(t,y,a,d) ,tspan,y0);\r\n[t1,y1] = ode23(@(t,y) loveRJ_ownEffect(t,y,a,b,c,d) ,tspan,y0);\r\n<\/pre><p>and we plot them all on the same axis<\/p><pre class=\"codeinput\">figure\r\nax1 = subplot(2,1,1); <span class=\"comment\">% simple plot<\/span>\r\nplot(ax1, t1,y1')\r\nax1.YLim = [-2 2];\r\nax1.XLabel.String = <span class=\"string\">'Time'<\/span>;\r\nax1.YLabel.String = <span class=\"string\">'Emotions'<\/span>;\r\nax1.Title.String = <span class=\"string\">'Romeo &amp; Juliet''s relationship - simple'<\/span>;\r\nlegend(<span class=\"string\">'Romeo'<\/span>,<span class=\"string\">'Juliet'<\/span>)\r\n\r\nax2 = subplot(2,1,2); <span class=\"comment\">% complex plot<\/span>\r\nplot(ax2, t2,y2')\r\nax2.YLim = [-1.5 1.5];\r\nax2.XLabel.String = <span class=\"string\">'Time'<\/span>;\r\nax2.YLabel.String = <span class=\"string\">'Emotions'<\/span>;\r\nax2.Title.String = <span class=\"string\">'Romeo &amp; Juliet''s relationship - self-effect(complex)'<\/span>;\r\nlegend(<span class=\"string\">'Romeo'<\/span>,<span class=\"string\">'Juliet'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2018\/RJ_1_1_LSnew_Final_02.png\" alt=\"\"> <p>Here we've seen the impact of different models on the evolution of emotions; it looks like the chosen values tend to reduce the frequency of changes in emotions. Let's call it self-control (or selfishness).<\/p><h4>Phase portrait and ODE solutions<a name=\"9cd6c125-d452-4d95-b445-d7248e766a48\"><\/a><\/h4><p>Phase plane analysis is a very important technique to study the behavior of dynamic systems; it covers a particularly relevant role in the nonlinear case, where widely applicable methods for computing analytical solutions are not available. In a nutshell, it basically consist of drawing the derivatives of solutions against the solutions in the phase plane. The derivatives of solutions are usually drawn in form of vector fields, to emphasize how large are changes in the solutions at a specific point in the phase plane and show the trajectories of the solutions, given specific initial conditions. Therefore, by superposing the two plots, it is possible to infer how the solutions might evolve, for the purposes to build our understanding under which conditions a system might stable or not.<\/p><p>We start calculating the derivatives <tt>y1'<\/tt> and <tt>y2'<\/tt> for each point of the phase plane. We create a grid of points where we want to draw out plots:<\/p><pre class=\"codeinput\">y1 = linspace(-10,10,20);\r\ny2 = linspace(-10,10,20);\r\n<\/pre><p><tt>meshgrid<\/tt> creates two matrices: one for all the uu-values of the grid, and one for all the vv-values of the grid. Then, we consider <tt>x1<\/tt> and <tt>x2<\/tt> matrices: <tt>x1<\/tt> will contain the value of <tt>y1'<\/tt> at each uu and vv position, while <tt>x2<\/tt> will contain the value of <tt>y2'<\/tt> at each uu and vv position of our grid.<\/p><pre class=\"codeinput\">[uu,vv] = meshgrid(y2,y1);\r\nx1 = zeros(size(uu));\r\nx2 = zeros(size(vv));\r\n<\/pre><p>Now we compute the vector field and plot the phase portrait. Our derivatives are computed for each point (y1, y2) at <tt>init_time = 0<\/tt>, through a <tt>for<\/tt> loop. We have obtained the phase portrait.<\/p><pre class=\"codeinput\">a = -2e-1;\r\nb = 5e-1;\r\ninit_time=0;\r\nloveRJ = @(t,y) [a*y(2); b*y(1)];\r\n\r\n<span class=\"keyword\">for<\/span> i = 1:numel(uu)\r\n    Yder = loveRJ(init_time,[uu(i); vv(i)]);\r\n    x1(i) = Yder(1);\r\n    x2(i) = Yder(2);\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>Finally we compute a couple of solutions and plot them, along with the phase portrait, on the phase plane.<\/p><pre class=\"codeinput\">figure\r\nquiver(gca,uu,vv,x1,x2,<span class=\"string\">'r'<\/span>);\r\nxlabel(<span class=\"string\">'Juliet Emotions'<\/span>);\r\nylabel(<span class=\"string\">'Romeo Emotions'<\/span>);\r\naxis <span class=\"string\">tight<\/span> <span class=\"string\">equal<\/span>;\r\n\r\n<span class=\"comment\">% Calculate and plot 1st Solution<\/span>\r\ntstart = 0;\r\ntfinal = 50;\r\ntspan = [tstart tfinal];\r\n\r\ny0_1 = [2;-1]; <span class=\"comment\">% initial conditions<\/span>\r\n[t,y1] = ode23(@(t,y) loveRJ(t,y),tspan,y0_1);\r\nfigure(gcf)\r\nhold <span class=\"string\">on<\/span>\r\nplot(y1(:,2),y1(:,1),<span class=\"string\">'b'<\/span>)\r\nplot(y1(1,2),y1(1,1),<span class=\"string\">'mo'<\/span>,<span class=\"keyword\">...<\/span><span class=\"comment\"> % starting point<\/span>\r\n<span class=\"string\">'MarkerEdgeColor'<\/span>,<span class=\"string\">'k'<\/span>,<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerFaceColor'<\/span>,[.49 1 .63],<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerSize'<\/span>,10)\r\nplot(y1(end,2),y1(end,1),<span class=\"string\">'bs'<\/span>,<span class=\"keyword\">...<\/span><span class=\"comment\"> % ending point<\/span>\r\n<span class=\"string\">'MarkerEdgeColor'<\/span>,<span class=\"string\">'k'<\/span>,<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerFaceColor'<\/span>,[.49 .63 1],<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerSize'<\/span>,10)\r\n\r\n<span class=\"comment\">% Calculate 2nd Solution<\/span>\r\ny0_2 = [4;1]; <span class=\"comment\">% initial conditions<\/span>\r\n[t,y2] = ode23(@(t,y) loveRJ(t,y),tspan,y0_2);figure(gcf)\r\nplot(y2(:,2),y2(:,1),<span class=\"string\">'c'<\/span>)\r\nplot(y2(1,2),y2(1,1),<span class=\"string\">'ko'<\/span>,<span class=\"keyword\">...<\/span><span class=\"comment\"> % starting point<\/span>\r\n<span class=\"string\">'MarkerEdgeColor'<\/span>,<span class=\"string\">'k'<\/span>,<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerFaceColor'<\/span>,[.49 1 .63],<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerSize'<\/span>,10)\r\nplot(y2(end,2),y2(end,1),<span class=\"string\">'bs'<\/span>,<span class=\"keyword\">...<\/span><span class=\"comment\"> % ending point<\/span>\r\n<span class=\"string\">'MarkerEdgeColor'<\/span>,<span class=\"string\">'k'<\/span>,<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerFaceColor'<\/span>,[.49 .63 1],<span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'MarkerSize'<\/span>,10)\r\nhold <span class=\"string\">off<\/span>\r\ntitle(<span class=\"string\">\"Two solutions plotted on vector field\"<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2018\/RJ_1_1_LSnew_Final_03.png\" alt=\"\"> <h4>Relationship Models<a name=\"4f31094f-ff5e-4e04-a5c2-6d880c558ea1\"><\/a><\/h4><p>This is the first model of their relationship, just dependent on the <tt>a<\/tt> \/ <tt>b<\/tt> parameters, defining Romeo's \/ Juliet's attraction for her\/his counterpart.<\/p><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> dy = loveRJ_simple(t,y,a,b)\r\ndR = a*y(2);\r\ndJ = b*y(1);\r\ndy = [dR;dJ];\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>In the second model of relationship, two additional parameters appear:<\/p><div><ul><li><tt>a,d<\/tt> = how Romeo, Juliet are influenced by their own feelings;<\/li><li><tt>b,c<\/tt> = Romeo's, Juliet's attraction for her\/his counterpart.<\/li><\/ul><\/div><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> dy = loveRJ_ownEffect(t,y,a,b,c,d)\r\ndR  = a*y(1) + b*y(2);\r\ndJ = c*y(1) + d*y(2);\r\ndy = [dR;dJ];\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4>Conclusion<a name=\"03450f15-bf0e-4a87-a588-7064121174a8\"><\/a><\/h4><p>Today we have taken some steps in the analysis of this ''special'' dynamic system. Others can be done, exploiting tools already available in MATLAB. For example, making the system even more realistic, using the same ODE solver would have been a good deal or would you have chosen another one and why? Do you think calculating eigenvalues would have shed some more light on your understanding of the system? How would you do it in MATLAB? Try to answer these questions and ask yourself some new ones! Let us know what you think <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=2811#respond\">here<\/a>.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_99d96fece93b48068ffbe5afcccb6e5e() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='99d96fece93b48068ffbe5afcccb6e5e ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 99d96fece93b48068ffbe5afcccb6e5e';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2018 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_99d96fece93b48068ffbe5afcccb6e5e()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2018a<br><\/p><\/div><!--\r\n99d96fece93b48068ffbe5afcccb6e5e ##### SOURCE BEGIN #####\r\n%% Of Love, Music, Shakespeare and Dynamic Systems\r\n% For any human being love is one of the biggest source of joy,\r\n% happiness...problems and puzzles. Today's guest blogger, Aldo Caraceto,\r\n% one of my fellow Application Engineers from Italy, is going to convince\r\n% you to approach it from a different angle.\r\n%% From Love to MATLAB  \r\n% I love music. One of my favourite songs is the popular Toto track _\"Hold\r\n% The Line\"_. Very briefly, this song is about love and timing, as you may\r\n% infer from the main chorus of the song: \"Hold the line, love isn't always\r\n% on time\". Ok, it's not always on time, but can I predict when it will be?\r\n% To answer this question, I've found a couple of insightful works, _\"Strogatz, S.H.\r\n% (1988), Love affairs and differential equations, Math. Magazine 61, 35\"_\r\n% and _\"Eva M. Strawbridge, Romeo and Juliet, a Dynamical Love Affair\"_\r\n% which I've translated in MATLAB code. The mentioned papers answer a very\r\n% simple question: how to express mathematically love dynamics, taking\r\n% inspiration from _\"Romeo & Juliet\"_ by W.Shakespeare. MATLAB may be a\r\n% good friend to help us in examining the evolution of this relationship.\r\n% An <https:\/\/www.mathworks.com\/help\/matlab\/ref\/ode23.html |ode23|> solver\r\n% has been used to solve the system of equations representing the story.\r\n% The results are quite funny.\r\n%% A very simple emotional mechanism\r\n% This first example examines the condition where each person\u00e2\u20ac&#x2122;s feelings\r\n% are only affected by the other. Besides, we want our Romeo not following\r\n% the expected behaviour: in this case, he'll be a confused and confusing\r\n% lover: the more Juliet loves him, the more he begins to dislike her,\r\n% generating surprising reactions with regards to the original tragedy.\r\n%\r\n% $\\frac{dr}{dt} = - a \\cdot j,$\r\n%\r\n% $\\frac{dj}{dt} = b  \\cdot r$\r\n%\r\n% In the system of differential equations representing the phenomenon, |a|\r\n% and |b| are real, positive numbers, and measured in 1\/time, hence are\r\n% frequencies: the bigger they are, the shorter in time will be the swing\r\n% of emotions for the two lovers. In addition, just a quick summary of the\r\n% variables represented in the plots:\r\n%%\r\n%\r\n% * |R(t)| = Romeo\u00e2\u20ac&#x2122;s feelings for Juliet at time t .\r\n% * |J(t)| = Juliet\u00e2\u20ac&#x2122;s feelings for Romeo at time t.\r\n% * |R(t), J(t) > 0| signify love, passion, and attraction.\r\n% * |R(t), J(t) < 0| values signify dislike.\r\n% * |R(t),J(t) == 0| signifies indifference.\r\n%%\r\na = 2e-1; \r\nb = 5e-1;\r\nloveRJ = @(t,y) [-a*y(2); b*y(1)];\r\n%%\r\n% Now, we can set the problem up, collect data, and plot results, adding\r\n% some annotations:\r\ntstart = 0;\r\ntfinal = 50;\r\ntspan = [tstart tfinal];\r\ny0 = [-1;2]; % initial conditions\r\n[t,y] = ode23(@(t,y) loveRJ(t,y),tspan,y0);\r\nplot(t,y')\r\nax = gca;\r\nax.XLim = [0 50];\r\nax.YLim = [-3 3];\r\nax.XLabel.String = 'Time';\r\nax.YLabel.String = 'Emotions';\r\nax.Title.String = 'Romeo & Juliet''s relationship';\r\nlegend('Romeo','Juliet')\r\n%% When things get complicated...\r\n% In our first example, we have described the story with the changing\r\n% emotions felt by Romeo and Juliet as they feed on each other's passion and\r\n% indifference while the time is passing by. As you know, things - and\r\n% people - are never so simple. Let's try to describe the relationship\r\n% between Romeo and Juliet a bit differently, trying to compare the first,\r\n% simple representation of the relationship with a more complex version.\r\n% Then, we'll plot them together.\r\n% \r\n% The new mathematical model can be expressed as: \r\n%\r\n% $\\frac{dr}{dt} = a \\cdot r + b \\cdot j ,$\r\n%\r\n% $\\frac{dj}{dt} = c \\cdot r + d \\cdot j$\r\n%\r\n% As you may notice, this time there are  four parameters, because we'd\r\n% like to consider also how Romeo and Juliet are influenced by their own\r\n% feelings:\r\n%%\r\na = -0.15;\r\nd = 0.17;\r\nb = 0.9; \r\nc = -0.9; \r\n%%\r\n% We calculate the solutions for the two relationship model (you may find\r\n% the model definitions at the end of the article):\r\ntstart = 0;\r\ntfinal = 50;\r\ntspan = [tstart tfinal];\r\ny0 = [1.00;-0.5]; % initial conditions\r\n\r\n[t2,y2] = ode23(@(t,y) loveRJ_simple(t,y,a,d) ,tspan,y0);\r\n[t1,y1] = ode23(@(t,y) loveRJ_ownEffect(t,y,a,b,c,d) ,tspan,y0);\r\n%% \r\n% and we plot them all on the same axis\r\n%\r\nfigure\r\nax1 = subplot(2,1,1); % simple plot\r\nplot(ax1, t1,y1')\r\nax1.YLim = [-2 2];\r\nax1.XLabel.String = 'Time';\r\nax1.YLabel.String = 'Emotions';\r\nax1.Title.String = 'Romeo & Juliet''s relationship - simple';\r\nlegend('Romeo','Juliet')\r\n\r\nax2 = subplot(2,1,2); % complex plot\r\nplot(ax2, t2,y2')\r\nax2.YLim = [-1.5 1.5];\r\nax2.XLabel.String = 'Time';\r\nax2.YLabel.String = 'Emotions';\r\nax2.Title.String = 'Romeo & Juliet''s relationship - self-effect(complex)';\r\nlegend('Romeo','Juliet')\r\n%%\r\n% Here we've seen the impact of different models on the evolution of\r\n% emotions; it looks like the chosen values tend to reduce the frequency of\r\n% changes in emotions. Let's call it self-control (or selfishness).\r\n%% Phase portrait and ODE solutions \r\n%\r\n% Phase plane analysis is a very important technique to study the behavior\r\n% of dynamic systems; it covers a particularly relevant role in the\r\n% nonlinear case, where widely applicable methods for computing analytical\r\n% solutions are not available. In a nutshell, it basically consist of\r\n% drawing the derivatives of solutions against the solutions in the phase\r\n% plane. The derivatives of solutions are usually drawn in form of vector\r\n% fields, to emphasize how large are changes in the solutions at a specific\r\n% point in the phase plane and show the trajectories of the solutions,\r\n% given specific initial conditions. Therefore, by superposing the two\r\n% plots, it is possible to infer how the solutions might evolve, for the\r\n% purposes to build our understanding under which conditions a system\r\n% might stable or not.\r\n%\r\n% We start calculating the derivatives |y1'| and |y2'| for each point of the phase plane.\r\n% We create a grid of points where we want to draw out plots:\r\ny1 = linspace(-10,10,20);\r\ny2 = linspace(-10,10,20);\r\n%%\r\n% |meshgrid| creates two matrices: one for all the uu-values of the grid, and\r\n% one for all the vv-values of the grid. Then, we consider |x1| and |x2|\r\n% matrices: |x1| will contain the value of |y1'| at each uu and vv position,\r\n% while |x2| will contain the value of |y2'| at each uu and vv position of\r\n% our grid.\r\n[uu,vv] = meshgrid(y2,y1);\r\nx1 = zeros(size(uu));\r\nx2 = zeros(size(vv));\r\n%%\r\n% Now we compute the vector field and plot the phase\r\n% portrait. Our derivatives are computed for each point (y1, y2) at\r\n% |init_time = 0|, through a |for| loop. We have obtained the phase\r\n% portrait.\r\na = -2e-1; \r\nb = 5e-1;\r\ninit_time=0;  \r\nloveRJ = @(t,y) [a*y(2); b*y(1)];\r\n\r\nfor i = 1:numel(uu)\r\n    Yder = loveRJ(init_time,[uu(i); vv(i)]);\r\n    x1(i) = Yder(1);\r\n    x2(i) = Yder(2);\r\nend\r\n%%  \r\n% Finally we compute a couple of solutions and plot them, along with the phase\r\n% portrait, on the phase plane.\r\nfigure\r\nquiver(gca,uu,vv,x1,x2,'r'); \r\nxlabel('Juliet Emotions');\r\nylabel('Romeo Emotions');\r\naxis tight equal;\r\n\r\n% Calculate and plot 1st Solution\r\ntstart = 0;\r\ntfinal = 50;\r\ntspan = [tstart tfinal];\r\n\r\ny0_1 = [2;-1]; % initial conditions\r\n[t,y1] = ode23(@(t,y) loveRJ(t,y),tspan,y0_1);\r\nfigure(gcf)\r\nhold on\r\nplot(y1(:,2),y1(:,1),'b')\r\nplot(y1(1,2),y1(1,1),'mo',... % starting point\r\n'MarkerEdgeColor','k',...\r\n    'MarkerFaceColor',[.49 1 .63],...\r\n    'MarkerSize',10)\r\nplot(y1(end,2),y1(end,1),'bs',... % ending point\r\n'MarkerEdgeColor','k',...\r\n    'MarkerFaceColor',[.49 .63 1],...\r\n    'MarkerSize',10)\r\n\r\n% Calculate 2nd Solution\r\ny0_2 = [4;1]; % initial conditions\r\n[t,y2] = ode23(@(t,y) loveRJ(t,y),tspan,y0_2);figure(gcf)\r\nplot(y2(:,2),y2(:,1),'c')\r\nplot(y2(1,2),y2(1,1),'ko',... % starting point\r\n'MarkerEdgeColor','k',...\r\n    'MarkerFaceColor',[.49 1 .63],...\r\n    'MarkerSize',10)\r\nplot(y2(end,2),y2(end,1),'bs',... % ending point\r\n'MarkerEdgeColor','k',...\r\n    'MarkerFaceColor',[.49 .63 1],...\r\n    'MarkerSize',10)\r\nhold off\r\ntitle(\"Two solutions plotted on vector field\")\r\n%% Relationship Models\r\n%\r\n% This is the first model of their relationship, just dependent on the |a|\r\n% \/ |b| parameters, defining Romeo's \/ Juliet's attraction for her\/his\r\n% counterpart.\r\n%%\r\nfunction dy = loveRJ_simple(t,y,a,b)\r\ndR = a*y(2);\r\ndJ = b*y(1);\r\ndy = [dR;dJ];\r\nend\r\n%%\r\n% In the second model of relationship, two additional parameters appear:\r\n%%\r\n% \r\n% * |a,d| = how Romeo, Juliet are influenced by their own feelings;\r\n% * |b,c| = Romeo's, Juliet's attraction for her\/his counterpart.\r\n%%\r\nfunction dy = loveRJ_ownEffect(t,y,a,b,c,d)\r\ndR  = a*y(1) + b*y(2);\r\ndJ = c*y(1) + d*y(2);\r\ndy = [dR;dJ];\r\nend\r\n%% Conclusion\r\n% Today we have taken some steps in the analysis of this ''special'' dynamic\r\n% system. Others can be done, exploiting tools already available in MATLAB.\r\n% For example, making the system even more realistic, using the same ODE solver would \r\n% have been a good deal or would you have chosen another one and why?\r\n% Do you think calculating eigenvalues would have shed some more light on \r\n% your understanding of the system? How would you do it in MATLAB?\r\n% Try to answer these questions and ask yourself some new ones! Let us know\r\n% what you think <https:\/\/blogs.mathworks.com\/loren\/?p=2811#respond here>.\r\n\r\n\r\n\r\n##### SOURCE END ##### 99d96fece93b48068ffbe5afcccb6e5e\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2018\/RJ_1_1_LSnew_Final_03.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>For any human being love is one of the biggest source of joy, happiness...problems and puzzles. Today's guest blogger, Aldo Caraceto, one of my fellow Application Engineers from Italy, is going to convince you to approach it from a different angle.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2018\/04\/27\/of-love-music-shakespeare-and-dynamic-systems\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[33],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2811"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=2811"}],"version-history":[{"count":3,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2811\/revisions"}],"predecessor-version":[{"id":2821,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2811\/revisions\/2821"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=2811"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=2811"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=2811"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}