[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

49. dynamics


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

49.1 Introduction to dynamics

The additional package dynamics includes several functions to create various graphical representations of discrete dynamical systems and fractals, and an implementation of the Runge-Kutta 4th-order numerical method for solving systems of differential equations.

To use the functions in this package you must first load it with load("dynamics").


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

49.2 Definitions for dynamics

Function: chaosgame ([[x1, y1]...[xm, ym]], [x0, y0], b, n, ...options...);

Implements the so-called chaos game: the initial point (x0, y0) is plotted and then one of the m points [x1, y1]...[xm, ym] will be selected at random. The next point plotted will be in the segment from the previous point to the point chosen randomly, at a fraction b of the distance from the random point. The procedure is repeated n times.

Function: evolution (F, y0, n,...options...);

Draws n+1 points in a two-dimensional graph, where the horizontal coordinates of the points are the integers 0, 1, 2, ..., n, and the vertical coordinates are the corresponding values y(n) of the sequence defined by the recurrence relation

 
        y(n+1) = F(y(n))

With initial value y(0) equal to y0. F must be an expression that depends only on the variable y (and not on n), y0 must be a real number and n must be a positive integer.

Function: evolution2d ([F, G], [x0, y0], n, ...options...);

Shows, in a two-dimensional plot, the first n+1 points in the sequence of points defined by the two-dimensional discrete dynamical system with recurrence relations

 
        x(n+1) = F(x(n), y(n))    y(n+1) = G(x(n), y(n))

With initial values x0 and y0. F and G must be two expressions that depend just on x and y.

Function: ifs ([r1,...,rm],[A1,...,Am], [[x1,y1]...[xm, ym]], [x0,y0], n, ...options...);

Implements the Iterated Function System method. This method is similar to the method described in the function chaosgame, but instead of shrinking the segment from the current point to the randomly chosen point, the 2 components of that segment will be multiplied by the 2 by 2 matrix Ai that corresponds to the point chosen randomly.

The random choice of one of the m attractive points can be made with a non-uniform probability distribution defined by the weights r1,...,rm. Those weights are given in cumulative form; for instance if there are 3 points with probabilities 0.2, 0.5 and 0.3, the weights r1, r2 and r3 could be 2, 7 and 10.

Function: orbits (F, y0, n1, n2, [x, x0, xf, xstep], ...options...);

Draws the orbits diagram for a family of one-dimensional discrete dynamical systems, with one parameter x; that kind of diagram is used to study the bifurcations of a one-dimensional discrete system.

The function F(y) defines a sequence with a starting value of y0, as in the case of the function evolution, but in this case that function will also depend on a parameter x that will take values in the interval from x0 to xf with increments of xstep. Each value used for the parameter x is shown on the horizontal axis. The vertical axis will show the n2 values of the sequence y(n1+1),..., y(n1+n2+1) obtained after letting the sequence evolve n1 iterations.

Function: rk (ODE, var, initial, domain)
Function: rk ([ODE1,...,ODEm], [v1,...,vm], [init1,...,initm], domain)

The first form solves numerically one first-order ordinary differential equation, and the second form solves a system of m of those equations, using the 4th order Runge-Kutta method. var represents the dependent variable. ODE must be an expression that depends only on the independent and dependent variables and represents the derivative of the dependent variable with respect to the independent variable.

The independent variable is specified with domain, which must be a list of four elements such as:

 
[t, 0, 10, 0.1]

the first element of the list identifies the independent variable, the second and third elements are the initial and final values for that variable, and the last element sets the increments that should be used within that interval.

If m equations are going to be solved, there should be m dependent variables v1, v2, ..., vm. The initial values for those variables will be init1, init2, ..., initm. There will still be just one independent variable defined by domain, as in the previous case. ODE1, ..., ODEm are the expressions for the derivatives of each dependent variable in terms of the independent variable. The only variables that may appear in those expressions are the independent variable and any of the dependent variables.

The result will be a list of lists with m+1 elements. Those m+1 elements will be the value of the independent variable, followed by the values of the dependent variables corresponding to that point in the interval of integration. If at some point one of the variables becomes too large, the list will stop there. Otherwise, the list will extend until the last value of the independent variable specified by domain.

Function: staircase (F, y0, n, ...options...);

Draws a staircase diagram for the sequence defined by the recurrence relation

 
        y(n+1) = F(y(n))

The interpretation and allowed values of the input parameters is the same as for the function evolution. A staircase diagram consists of a plot of the function F(y), together with the line G(y) = y. A vertical segment is drawn from the point (y0, y0) on that line until the point where it intersects the function F. From that point a horizontal segment is drawn until it reaches the point (y1, y1) on the line, and the procedure is repeated n times until the point (yn, yn) is reached.

Options

The options accepted by the functions that plot graphs are:

Examples

Graphical representation and staircase diagram for the sequence: 2, cos(2), cos(cos(2)),...

 
(%i1) load("dynamics")$
(%i2) evolution(cos(y), 2, 11, [yaxislabel, "y"], [xaxislabel,"n"]);
(%i3) staircase(cos(y), 1, 11, [domain, 0, 1.2]);

figures/dynamics1
figures/dynamics2

If your system is slow, you'll have to reduce the number of iterations in the following examples. And the pointsize that gives the best results depends on the monitor and the resolution being used.

Orbits diagram for the quadratic map

 
        y(n+1) = x + y(n)^2
 
(%i4) orbits(y^2+x, 0, 50, 200, [x, -2, 0.25, 0.01], [pointsize, 0.9]);

figures/dynamics3

To enlarge the region around the lower bifurcation near x = -1.25 use:

 
(%i5) orbits(x+y^2, 0, 100, 400, [x,-1,-1.53,-0.001], [pointsize,0.9],
             [ycenter,-1.2], [yradius,0.4]);

figures/dynamics4

Evolution of a two-dimensional system that leads to a fractal:

 
(%i6) f: 0.6*x*(1+2*x)+0.8*y*(x-1)-y^2-0.9$
(%i7) g: 0.1*x*(1-6*x+4*y)+0.1*y*(1+9*y)-0.4$
(%i8) evolution2d([f,g],[-0.5,0],50000,[pointsize,0.7]);

figures/dynamics5

And an enlargement of a small region in that fractal:

 
(%i9) evolution2d([f,g],[-0.5,0],300000,[pointsize,0.7], [xcenter,-0.7],
                  [ycenter,-0.3],[xradius,0.1],[yradius,0.1]);

figures/dynamics6

A plot of Sierpinsky's triangle, obtained with the chaos game:

 
(%i9) chaosgame([[0, 0], [1, 0], [0.5, sqrt(3)/2]], [0.1, 0.1], 1/2,
                 30000, [pointsize,0.7]);

figures/dynamics7

Barnsley's fern, obtained with an Iterated Function System:

 
(%i10) a1: matrix([0.85,0.04],[-0.04,0.85])$
(%i11) a2: matrix([0.2,-0.26],[0.23,0.22])$
(%i12) a3: matrix([-0.15,0.28],[0.26,0.24])$
(%i13) a4: matrix([0,0],[0,0.16])$
(%i14) p1: [0,1.6]$
(%i15) p2: [0,1.6]$
(%i16) p3: [0,0.44]$
(%i17) p4: [0,0]$
(%i18) w: [85,92,99,100]$
(%i19) ifs(w,[a1,a2,a3,a4],[p1,p2,p3,p4],[5,0],50000,[pointsize,0.9]);

figures/dynamics8

To solve numerically the differential equation

 
          dx/dt = t - x^2

With initial value x(t=0) = 1, in the interval of t from 0 to 8 and with increments of 0.1 for t, use:

 
(%i20) results: rk(t-x^2,x,1,[t,0,8,0.1])$

the results will be saved in the list results.

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

 
(%i21) sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$

[ << ] [ >> ]           [Top] [Contents] [Index] [ ? ]

This document was generated by root on October, 18 2006 using texi2html 1.76.