kill(all)$

maxima: graphics :programming : 'phase plane' from RK(=rk) output ,version 1.0
peter.vlasschaert@gmail.com,10/08/2019

load packages 
============

load(draw)$
load(lrats)$

→  how draw double y-axis
etc

ref :https://themaximalist.org/2016/03/10/numerical-solutions-of-2d-ode-systems-and-phase-plots/
=>  use research.

conceptual : get phase plane output from rungekutta(numerical methods)
========================================================
rk : maxima rungekutta method

1e) step : wxphaseplot2d
-----------------------------------

function : see ref : above.

/* wxphaseplot2d takes the output of s: rk([rhs1,rhs2]) 
   and  plots in phase space */
p1:wxphaseplot2d(s):=wxplot2d([discrete,makelist([p[2],p[3]],p,s)],[xlabel,""],[ylabel,""]);

info
change:wxphaseplot2d(s):=
makelist([p[2],p[3]] = > makelist([p[1],p[2]] ,'same plot
rem : see below example (/Brusselator)
rk: [p[1],p[2]] => time ,x
rk: [p[1],p[3]] => time ,y
rk: [p[2],p[3]] => x ,y

2e)  How get curve from rk data
--------------------------------------------
ref :https://en.wikipedia.org/wiki/Brusselator

(a:1,b:4);

p1a:'diff(x,t) = a+x^2*y-b*x-x;
p1b:'diff(y,t) =b*x-x^2*y;

how to use rk : maxima.
==================
syntax :          rk (<ODE>, <var>, <initial>, <domain>)
                       rk ([<ODE1>, ..., <ODEm>], [<v1>, ..., <vm>], [<init1>, ...,
                            <initm>], <domain>)

1e)  <ODE1> = x^2*y-5*x+1,<ODE1> =4*x-x^2*y
2e)  [<v1>,<v2>] =[x,y]
3e)  [<init>]=[1,1]
4e)  [ <domain>] =[t,0,200,0.1],simulation : t,start,end,step

info :rk

     To solve numerically the system:
                  dx/dt = 4-x^2-4*y^2     dy/dt = y^2-x^2+1
     for t between 0 and 4, and with values of -1.25 and 0.75 for x and
     y at t=0:
          (%i1) sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
          (%i2) plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$
     The plot will show the solution for variable <y> as a function of
     <t>.
  There are also some inexact matches for `rk'.

3e)  use rk
--------------

p2:sol : rk([rhs(p1a),rhs(p1b)],[x,y],[0,0],[t,0,50,.1])$

rem :rk output

 p2a:rk([rhs(p1a),rhs(p1b)],[x,y],[0,0],[t,0,5,0.7]);

4e) use function :wxphaseplot2d(s)
------------------------------------------------

p3:wxphaseplot2d(sol);

intermezzo:

p5:[exp(-t),-exp(-t),exp(-2*t),-2*exp(-2*t),exp(-t)-exp(-2*t),-exp(-t)+exp(-2*t)];

ev(p5[1],t=0);

m:makelist([ev(p5[1],nouns),ev(p5[2],nouns)],t,0,10,1),float;

ma:float(m);

m1:makelist([ev(p5[1],nouns),ev(p5[2],nouns)],t,0,100,.1),float$
m2:makelist([ev(p5[3],nouns),ev(p5[4],nouns)],t,0,100,.1),float$
m3:makelist([ev(p5[5],nouns),ev(p5[6],nouns)],t,0,100,.1),float$

p6:wxdraw2d(point_size=2,
color =red,
points(m1),
color = blue,
points(m2),
color = green,
points(m3));

closed : intermezzo

we use :intermezzo
===============
init :[0,0] => red
init  [1,1] => blue

rem : rk ≠ initial conditions
====================
p7a:   [0,0]
p7b,c:[1,1]

p7a:sol1:rk([rhs(p1a),rhs(p1b)],[x,y],[0,0],[t,0,50,.1])$
p7b:sol2:rk([rhs(p1a),rhs(p1b)],[x,y],[1,1],[t,0,50,.1])$
p7c:sol3:rk([rhs(p1a),rhs(p1b)],[x,y],[1,1],[t,0,50,.1])$

p7c1:sol3a:rk([rhs(p1a),rhs(p1b)],[x,y],[1,1],[t,0,10,.5]);

plot : time versus x
==============

p8:wxdraw2d(point_size=.5,    
color =red,
points(sol1),
color = blue,
points(sol2));

definition :time versus x
==================

plottimex(s):=wxplot2d([discrete,makelist([p[1],p[2]],p,s)],[xlabel,""],[ylabel,""]);

p9a:plottimex(sol1);

definition :time versus  y
==================

plottimey(s):=wxplot2d([discrete,makelist([p[1],p[3]],p,s)],[xlabel,""],[ylabel,""]);

p9b:plottimey(sol1);

combiplot 
=======

p10a:makelist([p[1],p[2]],p,sol)$

p10b:makelist([p[1],p[3]],p,sol)$

1e) method:combiplot

p10:wxplot2d([[discrete,p10a],[discrete,
        p10b]],[style,[lines,1,1],[points,0.3]],[color,red,blue],
    [xlabel,"time"],[ylabel,"data"],[legend,"time_xplot","time_yplot"]);

;

definition :x versus  y  (  =wxphaseplot2d(s) ,see above 'p1')
=============================================

p11:plotxy(s):=wxplot2d([discrete,makelist([p[2],p[3]],p,s)],[xlabel,""],[ylabel,""]);

p9c:plotxy(sol1);

2e) method:combiplot

how to do : combiplot with use 'wxdraw2d'
================================
we use : p10a,p10b , 'see above'

p12:wxdraw2d(point_size=.5,   
color =red,
            points_joined =true,   
points(p10a),          
color = blue,
            points_joined = false,
points(p10b));

rem: wx before will be a inline plot: 'wxdraw2d' before the draw2d.

how to plot legend to 'p12'
====================

p12a:wxdraw2d(point_size=.5,   
color =red,
 key = "red:time_xplot",
points(p10a),
color = blue,
 key = "blue:time_yplot",   
points(p10b)
  
);

rem :double y axis
==============
ref :https://themaximalist.org/tag/wxdraw2d/

p12b:wxdraw2d(
    /* general : option */
                    grid =[1,1],
                    line_width = 3,
    /*************************/
                    color = red,                                   
                    key = "cos(2*x)",
                              explicit(cos(2*x),x,-2*π,+2*π),
                   color = blue,
       /* double y axis */
    yaxis_secondary = true,
    ytics_secondary = true,
                  key="cos(2*x)*exp(-2*x)",
                               explicit(cos(2*x)*exp(-2*x),x,-2*π,+2*π),
    /* axis labels */
xlabel = "time",
ylabel = " function"
                           );

rem : to extract data from output  Runge kutta ( Maxima = rk( ..,..,..,..))
=====================================================
use : list instruction , 'maxima'
command : first ,second ,last ,third

how to do this : rk  ( use : p7c1, 'see above')
==================================

1e) step (output : simulation with rk)
===========================

data structure (list), output  => rk
----------------------------------------------

see : info,Brusselator,data ( odes),'see above'

p13:rk([rhs(p1a),rhs(p1b)],[x,y],[1,1],[t,0,5,.5]);

rem :[[time1,x1,y1],[time2,x2,y2],...................................]

2e)  combination : first ( =time) & second ( =x).
===================================

p13a:makelist([first(p13[i]),second(p13[i])],i,1,length(p13));

3e)  combination : first ( =time) & third ( =y).
===================================

p13b:makelist([first(p13[i]),third(p13[i])],i,1,length(p13));

4e)  combination : second ( =x) & thrid ( =y).
====================================

1e) method: ( third)
--------------------------

p13c1:makelist([second(p13[i]),third(p13[i])],i,1,length(p13));

2e) method (last)
------------------------

p13c2:makelist([second(p13[i]),last(p13[i])],i,1,length(p13));