A First Look 
at
Discrete Dynamical Systems

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

Special Mathematica Functions used in this NoteBook

Needs["Graphics`Colors`"];

Notebook Overview

This notebook will allow you to experiment with discrete dynamical systems.  The "change form" of such systems is:
	
					p[0] = P0
					p[t+1] = p[t] + f[p[t]]
					
where you define the change function f[p_] := ??? as in the samples.  The change in  p[.]  as  t  goes from  t  to  t+1  is p[t+1] - p[t], on one hand, and given by the function  f[.], on the other,

			p[t+1] - p[t] = f[ p[t] ]   <=>   		p[t+1] = p[t] + f[p[t]]
			
	The function iteration form that is used with "cobwebs" is:
	
					p[0] = P0
					p[t+1] = g[p[t]]
					
where g[p_] := p + f[p]  or  f[p_] := g[p] - p.  Both of these forms are useful.  The change form is more like Euer's approximation to solution of differential equations.

References to the Text

This NoteBook goes with Chapter 20 of the text.  Also see the CobWeb NoteBook.

Step-by-step: Simple Iteration of Numbers

Automatic Recursion with a Mathematica Rule

Suppose we have

p[0] = P0
	p[t+1] = p[t] + f[p[t]]

with

f[p] = - p/3

Solution of this equation can be done recursively.  The initial starting value is given, so we can find  p[1]  by substituting  p[0]  on the right side of the change formula.  Then we can find  p[2]  by substituting  p[1]  on the right hand side, etc.  Let's do a few steps by hand, so we can check the computer work to follow.  Suppose we begin with

p[0] = P0 = 30

p[1] = p[0] + f[p[0]] = 30 - 30/3 = 20

p[2] = p[1] + f[p[1]] = 20 - 20/3 = 40/3

p[3] = p[2] + f[p[2]] = 40/3 - 40/9 = 80/9

Mathematica is based on "re-writing rules," so if we simply tell Mathematica the rules (at  t  looking back at  t-1)  and ask for  p[3], it will look up p[3-1], if it knows it stop, otherwise, look up  p[2-1], etc. until it comes to  p[0] which it knows.  Enter the next two samples:

Clear[p, t, f]; p[0] = P0p[t_] := p[t - 1] + f[p[t - 1]]p[1]p[2]p[3]

f[p_] = -p/3;
P0 = 30
p[1]
p[2]
p[3]

Mathematica will automatically compute the recursive steps, for example, without computing  p[4]  explicit, ask for  p[5] next:

p[5]

The rule approach to discrete dynamical systems is very simple - it leaves all the work to Mathematica, but to plot the solution sequence, we need to create a list as follows:

ListPlot[ Table[p[t],{t,1,10}] ];

If you prefer floating point values instead of exact fractions, simply define  f[.]  numerically as in the next sample:

Clear[p,t,f];
p[0] = 30.0
p[t_] := p[t-1] + f[p[t-1]]
f[p_] := -0.25 p
p[1]
p[2]
p[3]
ListPlot[ Table[p[t],{t,1,10}] ];

Practice Exercise

Define and plot the sequence defined by

p[0] = 30
	p[t + 1] = p[t] + f[ p[t] ]

where

f[p] = +p/5

Compute the first few values by hand to be sure your program works properly.

NestList[.] on the Iteration function  g[p] = p + f[p]

Let's return to the example

p[0] = P0
	p[t+1] = p[t] + f[p[t]]
with
	f[p] = - p/3

which is solved recursively by

p[0] = P0 = 30

p[1] = p[0] + f[p[0]] = 30 - 30/3 = 20

p[2] = p[1] + f[p[1]] = 20 - 20/3 = 40/3

p[3] = p[2] + f[p[2]] = 40/3 - 40/9 = 80/9

Notice that each step applys the same PROCESS or function to the output of the previous step

p[t] = g[g[g[...t times ...g[P0] ...]]]
where
	g[p] = p + f[p] = p - p/3 = 2p/3 in the example above.

We may use the formula   p[t+1] = g[p[t]] = 2 p[t]/3 to compute p[1]

p[1] = g[p[0]] = 2p[0]/3 = 20

Use the formula   p[t+1] = g[p[t]] = p[t]/2  again to compute p[2]

p[2] = g[p[1]] = 2p[1]/3 = 40/3

Mathematica has a single command to make a list of values arising from successively applying the same function over and over.  <Enter> the computation in the next cell and check the iterations

f[p_] := -(p/3.0); g[p_] := p + f[p]; iterations = 10; P0 = 30iterates = NestList[g, P0, iterations]plot1 = 
      ListPlot[iterates, PlotJoined -> True, PlotStyle -> {Red}]; 
    P0 = 60iterates = NestList[g, P0, iterations]plot2 = 
      ListPlot[iterates, PlotJoined -> True, PlotStyle -> {Blue}]; 
    Show[plot1, plot2];

Practice Exercise

Use the NestList command as above to compute  10  iterates of the difference equation
	p[t+1] = p[t] - 1.5 p[t]
starting with the initial values
	p[0] = 30.0
and
	p[0] = 60.0

\(f[p_] := 
      \(?? 
g[p_]:=p+??
P0 = 
          \(?? 
iterates=NestList[??]

              \(ListPlot[iterates, PlotJoined -> True]; \)\)\)\)

A Nonlinear Example

The previous examples used very simple change functions.  Here is a more complicated example:

p[0] = P0
	p[t + 1] = p[t] + (1/3) p[t] (1 - p[t]^2)

so if we start at  P0 = 1/2,

p[1] = (1/2) + (1/2)(1 - 1/4)/3 = 5/8

p[2] = (5/8) + (5/8)(1 - 25/64)/3 = 385/512

p[3] = 115545045/134217728
			...

f[p_] := (1/3) p (1 - (p^2)); g[p_] := p + f[p]; iterations = 10; P0 = 0.5; iterates = NestList[g, P0, iterations]plot1 = 
      ListPlot[iterates, PlotJoined -> True, PlotRange -> {0, 2}, 
        PlotStyle -> {Red}]; P0 = 1.5; iterates = NestList[g, P0, iterations]plot2 = 
      ListPlot[iterates, PlotJoined -> True, PlotRange -> {0, 2}, 
        PlotStyle -> {Blue}]; Show[plot1, plot2];

Practice Exercise

Plot 15 iterations of the system
	p[0] = P0
	p[t + 1] = p[t] + (1/3) p[t] (1 - p[t]^2)
starting at  p[0] = 0.1 and p[0] = 1.9.

Exercises

Exercise CD20.1.1: A first look at stability of p[t+1] = p[t] + a p[t] + b

Experiment with your choice of calculations in this NoteBook to formulate a conjecture about the role of a and b in determining both the stability and limiting value of the general affine linear dynamical system.  Consider the cases:
	
	(1)  a = -0.2, 	b = 100,	P0 = 100, 500, 1000
	(2)  a = +0.2,	b = 100,	P0 = 100, 500, 1000
	(3)  a = -0.2,	b = 50,		P0 = 100, 500, 1000
	(4)  a = +0.2,	b = 50,		P0 = 100, 500, 1000
	(5)  a = -0.1,	b = 100,	P0 = 100, 500, 1000
	(6)  a = +0.1,	b = 100,	P0 = 100, 500, 1000
	(7)  a = -1.1,	b = 110,	P0 = 100, 500, 1000
	(8)  a = -2.0,	b = 200,	P0 = 100, 500, 1000
	(9)  a = -2.2,	b = 220,	P0 = 100, 500, 1000

Exercise CD20.1.2

Compute and plot  5  steps of the nonlinear system

p[0] = P0
	p[t + 1] = p[t] + p[t] (1 - p[t]/50)

by hand for  P0 = 25  and repeat the computation with Mathematica.  Then compute for  P0 = 25, 40, 60, 75  and plot 30 steps for each initial value.