share/dynamics/dynamics.mac

   1 /************************************************************************
   2     Copyright (C) 2003, 2004, 2006, 2007 Jaime E. Villate <villate@fe.up.pt>
   3
   4     This program is free software; you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation; either version 2 of the License, or
   7     (at your option) any later version.
   8
   9     This program is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with this program; if not, write to the Free Software
  16     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  17
  18 $Id: dynamics.mac,v 1.6 2010-04-20 12:55:46 villate Exp $
  19 ***********************************************************************/
  20
  21 load("complex_dynamics.lisp")$
  22 load("plotdf.lisp")$
  23
  24 /************************************************************************
  25  Function evolution
  26
  27     Description: Maxima program to draw orbits of one-dimensional mappings
  28     Usage:
  29       evolution(fun, initial, n, ..., [options], ...);
  30     "fun" must be an expression which depends only on one variable.
  31     "initial" is the initial value for that variable, and n is the number
  32     of steps. [options] can be any options accepted by plot2d.
  33 */
  34
  35 evolution(fun, initial, n, [options]) :=
  36   block
  37   ([z:initial, data: [[0, initial]], kwds: [], plot: numer: true, float: true],
  38     if length(listofvars(fun)) # 1 then
  39            error("fun should depend on one variable"),
  40     for i thru n do
  41       (z: ev(fun, listofvars(fun)[1]=z),
  42        if not numberp(z) then
  43                  error("The function gave a non-numerical value:", z),
  44        data: cons([i, z], data)),
  45     for i thru length(options) do
  46       (if not listp(options[i]) then error("Wrong option", options[i]),
  47        kwds: cons(options[i][1], kwds)),
  48     if not member('x,kwds) then options: cons(['x, -0.5, n+0.5],options),
  49     if not member('xlabel,kwds) then options: endcons(['xlabel, "n"],options),
  50     if not member('ylabel,kwds) then
  51           options: endcons(['ylabel,sconcat(listofvars(fun)[1],"(n)")],options),
  52     if not member('style,kwds) then options: endcons(['style,'points],options),
  53     options: cons(['discrete, data], options),
  54     apply(plot2d, options))$
  55
  56 /************************************************************************
  57  Function staircase
  58
  59     Description: Maxima program to draw staircase diagrams of
  60                  one-dimensional mappings
  61     Usage:
  62       staircase(fun, initial, n, [options]);
  63     "fun" must be an expression which depends only on one variable.
  64     "initial" is the initial value for that variable, and n is the number
  65     of steps. [options] can be any options accepted by plot2d.
  66 */
  67
  68 staircase(fun, initial, n, [options]) :=
  69   block
  70   ([zf, z:initial, z1:initial, z2:initial, stair:[[initial, initial]],
  71     kwds: [], numer: true, float: true],
  72     if length(listofvars(fun)) # 1 then
  73            error("fun should depend on one variable"),
  74     for i thru n do
  75       (zf: ev(fun, listofvars(fun)[1]=z),
  76        if not numberp(zf) then
  77                  error("The function gave a non-numerical value:", zf),
  78        stair: append(stair, [[z, zf], [zf, zf]]),
  79        z: zf,
  80        if z < z1 then z1: z,
  81        if z > z2 then z2: z),
  82     for i thru length(options) do
  83       (if not listp(options[i]) then error("Wrong option", options[i]),
  84        kwds: cons(options[i][1], kwds)),
  85     if not member(listofvars(fun)[1],kwds) then
  86        options: cons([listofvars(fun)[1],
  87                         4*(1.1*z1-.1*z2)/3,4*(1.1*z2-.1*z1)/3],options),
  88     if not member('y,kwds) then
  89        options: endcons(['y,1.1*z1-.1*z2,1.1*z2-.1*z1],options),
  90     if not member('xlabel,kwds) then
  91        options: endcons(['xlabel,sconcat(listofvars(fun)[1],"(n)")],options),
  92     if not member('ylabel,kwds) then
  93        options: endcons(['ylabel,sconcat(listofvars(fun)[1],"(n+1)")],options),
  94     if not member('legend,kwds) then options: endcons(['legend,false],options),
  95     options: cons([['discrete, stair], fun, listofvars(fun)[1]], options),
  96     apply(plot2d, options))$
  97
  98 /************************************************************************
  99  Function evolution2d
 100
 101     Description: Maxima program to show the evolution of a two-dimensional
 102                  discrete dynamical system in a 2-dimensional graph.
 103     Usage:
 104       evolution2d([fun1, fun2], [u, v], [u0, v0], n, ..., [options], ...);
 105
 106     fun1 and fun2 must be expressions which depend only on two variables
 107     (u and v) named on the second list. The third list must give two
 108     initial numerical values for the two variables. n is the number of
 109     points to be shown. The options can be any accepted by plot2d.
 110 */
 111
 112 evolution2d(fun, state, initial, n, [options]) :=
 113   block
 114   ([x: initial[1], y: initial[2], xaux, data: [],
 115     fx: false, fy: false, kwds: [], numer: true, float: true],
 116     if length(listofvars(fun)) >2 then
 117            error("fun should depend on two variables."),
 118     for i thru length(options) do
 119       (if not listp(options[i]) then error("Wrong option", options[i]),
 120        if options[i][1]=state[1] then
 121           (x1: options[i][2],
 122            x2: options[i][3],
 123            options[i][1]:'x,
 124            fx: true),
 125        if options[i][1]=state[2] then
 126           (y1: options[i][2],
 127            y2: options[i][3],
 128            options[i][1]: 'y,
 129            fy: true),
 130        kwds: cons(options[i][1], kwds)),
 131     for i thru n do
 132       (if (fx and fy)
 133         then
 134           (if ((y >= y1) and (y <= y2) and (x >= x1) and (x <= y2)) then
 135               data: cons([x, y], data))
 136         else
 137           (data: cons([x, y], data)),
 138        xaux: ev(fun[1], state[1]=x, state[2]=y),
 139        if not numberp(xaux) then
 140                  error("Function",fun[1],"gave a non-numerical value:", xaux),
 141        y: ev(fun[2], state[1]=x, state[2]=y),
 142        if not numberp(y) then
 143                  error("The function",fun[2],"gave a non-numerical value:", y),
 144        x: xaux),
 145     if not member('xlabel,kwds) then
 146        options: endcons(['xlabel,string(state[1])],options),
 147     if not member('ylabel,kwds) then
 148        options: endcons(['ylabel,string(state[2])],options),
 149     if not member('style,kwds) then options: endcons(['style,'points],options),
 150     options: cons(['discrete, data], options),
 151     print("Graph passed to plot2d..."),
 152     apply(plot2d, options))$
 153
 154 /************************************************************************
 155  Function chaosgame
 156
 157     Description: Maxima program to play the "chaos game".
 158     Usage:
 159       chaosgame([[x1,y1]...[xm,ym]], [x0,y0], beta, n, [options]);
 160
 161     where [x0,y0] are the coordinates of the initial point, [x1,y1]
 162     ... [xm,ym] are the attracting points, beta is the ratio of the
 163     final to initial distance from the attracting point, and n the
 164     number of points to show.
 165 */
 166
 167 chaosgame(point,p0,b,n,[options]) :=
 168 block
 169   ([p:ev(p0,numer), j, m:length(point), data: [p0], kwds: [],
 170     xlabel, ylabel, numer:true, float:true],
 171     if length(p0) # 2 then
 172         error("Initial point",p0,"should be a list with 2 components"),
 173     if not (numberp(p[1]) and numberp(p[2])) then
 174         error("The components of the initial point", p,"should be numbers"),
 175     for i thru m do
 176         if length(point[i]) # 2 then
 177          error("Point",point[i],"should be a list with 2 components"),
 178     for i thru n do
 179        (j: random(m) + 1,
 180         p: ev(point[j] + b*(p-point[j]), numer),
 181         if not (numberp(p[1]) and numberp(p[2])) then
 182               error("Non-numerical coordinates were obtained:", p),
 183        data: cons(p, data)),
 184     for i thru length(options) do
 185       (if not listp(options[i]) then error("Wrong option", options[i]),
 186        kwds: cons(options[i][1], kwds)),
 187     if not member('xlabel,kwds) then
 188        (xlabel: sconcat("The chaos game with ",m," points"),
 189         options: endcons(['xlabel,xlabel],options)),
 190     if not member('ylabel,kwds) then
 191        (ylabel: sconcat("contraction factor: ",b),
 192         options: endcons(['ylabel,ylabel],options)),
 193     if not member('style,kwds) then options: endcons(['style,'points],options),
 194     options: cons(['discrete, data], options),
 195     apply(plot2d, options))$
 196
 197 /************************************************************************
 198  Function ifs
 199
 200     Description: Maxima program to create 2-d fractals using Barnley's
 201                  IFS (Iterated Function Systems) method.
 202     Usage:
 203       ifs([r1,...,rm], [A1,...,Am], [p1,...,pm]], p0, n, [options]);
 204
 205     where r1,...,rm are cumulative weights for the attracting points.
 206     For instance if there are 3 points with probabilities 0.2, 0.5 and
 207     0.3, you can use [2,7,10].
 208     A1,...,Am, are the matrices for the m attracting points and
 209     p1,...,pm are the 2D coordinates of those points. p0 are the 2D
 210     coordinates of the initial point and n the number of points to be
 211     shown.
 212 */
 213
 214 ifs(prob, mat, point, p0, n, [options]) :=
 215 block
 216   ([p:ev(p0,numer), s, r, w:last(prob), data:[p0], xlabel, kwds: [],
 217     m:length(prob), plist, numer:true, float:true],
 218     if length(p0) # 2 then
 219         error("Initial point",p0,"should be a list with 2 components"),
 220     if not (numberp(p[1]) and numberp(p[2])) then
 221         error("The components of the initial point", p,"should be numbers"),
 222     if not ((length(mat)=m) and (length(point)=m)) then
 223         error("There should be the same number of probabilities, matrices and points"),
 224     for i thru m do
 225        (if length(point[i]) # 2 then
 226            error("Point",point[i],"should be a list with 2 components"),
 227         if not numberp(prob[i]) then
 228            error("Cumulative probability",prob[i],"should be a number"),
 229         if not ((length(mat[i])=2) and (length(mat[i][1])=2)) then
 230            error("Matrix",mat[i],"should be a 2x2 matrix")),
 231     for i thru n do
 232        (r: random(w) + 1,
 233         s: 0,
 234         for j while (r-prob[j]) > 0 do s:j,
 235         p: ev(mat[s+1].p + point[s+1], numer),
 236         plist: first(transpose(p)),
 237         if not (numberp(plist[1]) and numberp(plist[2])) then
 238            error("Non-numerical coordinates were obtained:", plist),
 239         data: cons(plist, data)),
 240     for i thru length(options) do
 241        (if not listp(options[i]) then error("Wrong option", options[i]),
 242         kwds: cons(options[i][1], kwds)),
 243     if not member('xlabel,kwds) then
 244        (xlabel: sconcat("Iterated Function System of ",m," transformations"),
 245         options: endcons(['xlabel,xlabel],options)),
 246     if not member('ylabel,kwds) then options: endcons(['ylabel,""],options),
 247     if not member('style,kwds) then options: endcons(['style,'points],options),
 248     options: cons(['discrete, data], options),
 249     apply(plot2d, options))$
 250
 251 /************************************************************************
 252  Function orbits
 253
 254     Description: Maxima program to draw orbit diagrams (also dubbed as
 255                  bifurcation diagrams) for first-order discrete
 256                  dynamical systems
 257     Usage:
 258       orbits(f, initial, n1, n2, [u, u0, uf], ..., [options], ...)
 259
 260     It will plot the orbits diagram for a family of first order discrete
 261     dynamical system, with a parameter.
 262     f must be an expression which depends only on one variable and
 263     the parameter. The parameter must be named in the first component of
 264     the list [u, u0, uf] (in this example it would be u) that sets up
 265     the range of the variation for the parameter, from u0, uf. The initial
 266     value of the parameter can be smaller than the final value, in the
 267     cases when the bifurcations appear as the parameter decreases.
 268     "initial" is an initial value for the state variable.
 269     The parameter will be represented in the horizontal axis; the range of
 270     the parameter will be divided by the value given by the nticks option,
 271     or by 200 if that option is not used. For each value of the parameter,
 272     the system will be left to evolve during n1 steps, and then the next
 273     n2 steps will be plotted in the vertical axis.
 274     The program tuns faster if a range is given to the vertical axis.
 275     By default, the vertical scale is divide in 600 pixels. That number
 276     can be changed with the option pixels, for instance, [pixels, 700].
 277 */
 278
 279 orbits(f, y0, n1, n2, domain, [options]) :=
 280 block
 281  ([x, y, data: [], y1, y2, yflag: false, kwds: [], var, points, plot: [],
 282    x0: domain[2], s, n, pixels: 600, numer:true, float: true],
 283    var: listofvars(f),
 284    var: delete(domain[1],var),
 285     if length(listofvars(fun)) >2 then
 286            error("f should depend on two variables."),
 287     for i thru length(options) do
 288       (if not listp(options[i]) then error("Wrong option", options[i]),
 289        if options[i][1]=var[1] then
 290           (y1: options[i][2],
 291            y2: options[i][3],
 292            if not y1 < y2 then
 293            error("The first number has to be smaller than the second in:",options[i]),
 294            options[i][1]: 'y,
 295            yflag: true),
 296        if options[i][1]='nticks then
 297           (if integerp(options[i][2]) then
 298               n: options[i][2]
 299            else
 300               error("Option nticks should have an integer value")),
 301        if options[i][1]='pixels then
 302           (pixels: options[i][2],
 303            delete(options[i],options))
 304        else
 305           plot: endcons(options[i], plot),
 306        kwds: cons(options[i][1], kwds)),
 307    if  not (member('nticks,kwds) and integerp(n)) then n: 200,
 308    s: (domain[3] - domain[2])/n,
 309    if yflag then
 310      (for i:0 thru n do
 311          (x: x0+i*s,
 312           y: y0,
 313           points: [],
 314           for j thru n1 do
 315              (y: ev(f, domain[1]=x, var[1]=y),
 316               if not numberp(y) then
 317                  error("The function gave a non-numerical value:", y)),
 318           for j thru n2 do
 319              (y: ev(f, domain[1]=x, var[1]=y),
 320               if ((y >= y1) and (y <= y2)) then
 321                  (k: entier(pixels*(y-y1)/(y2-y1))+1,
 322                   if k>pixels then k: pixels,
 323                   if not member(k, points) then points: cons(k, points))),
 324           for j thru length(points) do
 325               data: cons([x,(y2-y1)*(points[j]-0.5)/pixels+y1], data)))
 326    else
 327      (for i:0 thru n do
 328          (x: x0+i*s,
 329           y: y0,
 330           for j thru (n1+n2)/2 do
 331              (y: ev(f, domain[1]=x, var[1]=y),
 332               if not numberp(y) then
 333                  error("The function gave a non-numerical value:", y),
 334               if not yflag then
 335                  (y1:y, y2:y, yflag: true)
 336               else
 337                  (if y<y1 then y1: y,
 338                   if y>y2 then y2: y))),
 339       if not y1 < y2 then error("It was not possible to find any points"),
 340       for i:0 thru n do
 341          (x: x0+i*s,
 342           y: y0,
 343           points: [],
 344           for j thru n1 do y: ev(f, domain[1]=x, var[1]=y),
 345           for j thru n2 do
 346              (y: ev(f, domain[1]=x, var[1]=y),
 347               if ((y >= y1) and (y <= y2)) then
 348                  (k: entier(pixels*(y-y1)/(y2-y1))+1,
 349                   if k>pixels then k: pixels,
 350                   if not member(k, points) then points: cons(k, points))),
 351           for j thru length(points) do
 352               data: cons([x,(y2-y1)*(points[j]-0.5)/pixels+y1], data))),
 353    if (length(data) > 0) then
 354       (if not member('x,kwds) then
 355           (if domain[3] > domain[2] then
 356               plot: cons(['x,domain[2],domain[3]],plot)
 357            else
 358               plot: cons(['x,domain[3],domain[2]],plot)),
 359        if not member('xlabel,kwds) then
 360           plot: endcons(['xlabel,string(domain[1])],plot),
 361        if not member('ylabel,kwds) then
 362           plot: endcons(['ylabel,string(var[1])],plot),
 363        if not member('style,kwds) then plot:endcons(['style,'points],plot),
 364        plot: cons(['discrete, data], plot),
 365        print("Graph passed to plot2d..."),
 366        apply(plot2d, plot)))$
 367