<< DiscreteMath`RSolve`

Discrete Systems and Eigenvalues

Mathematica is able to solve simple discrete linear systems using RSolve:

Clear[x, y]
req1  =  
      x[n + 1]  ==   -4 
            x[n] + 2  y[n] 
req2  =  
      y[n + 1]  ==   x[n]  +  2  y[n]
RSolve[{req1, req2}, {x[n], y[n]}, n]

However, more complicated nonlinear systems will usually be too difficult for RSolve to solve:

req1  =  
      x[n + 1]  ==   -4 x[n] + 2  y[n]  +  
          x[n] y[n]
req2  =  
      y[n + 1]  ==   x[n]  +  2  y[n] - 
          4  x[n] y[n]
RSolve[{req1, req2}, {x[n], y[n]}, n]

This is the discrete dynamical system which corresponds to the one we looked at in Lecture 10 for continuous dynamical systems.

We are forced to resort to constructing sequences of solutions and analysing those. We can make linear approximations near the fixed points, and try to determine stability based on the eigenvalues of a matrix.

Fixed Points

The fixed points for the system will satisfy

req1  =  
      x  ==   -4 x + 2  y  +  x y
req2  =  
      y  ==   x  + 
          2  y - 4  x y 
Solve[{req1, req2}, {x, y}]

Which, not surprisingly, are the same fixed points we found for the continuous system.

Solution

Clear[x, y]
num = 10;
x[0]  :=  7/19;
y[0]  :=  7/9 + 0.001;
x[n_]  :=  x[n]  = -4 x[n - 1] + 2  y[n - 1]  +  
          x[n - 1] y[n - 1]
y[n_]  :=  y[n]  = 
        x[n - 1]  +  2  y[n - 1] - 
          4  x[n - 1] y[n - 1]
data  =  
        Table[{n, x[n], y[n]}, {n, 0, num}];

    
phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
datax1  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 2]}, {n, 0, 
            num}];
datax2  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 3]}, {n, 0, 
            num}];
plotphase  =  
      ListPlot[phase, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"x1", "x2"}]
plotphase  =  
      ListPlot[datax1, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"n", "x1"}]
plotphase  =  
      ListPlot[datax2, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"n", "x2"}]

matrix  =  {{-4, 2}, {1, 2}}
MatrixForm[%]
Eigenvalues[matrix]
N[%]

matrix  =  {{-29/9, 405/171}, {-361/171, 
          10/19}}
MatrixForm[%]
Eigenvalues[matrix]
N[%]

Sqrt[(-1.35)^2 + (1.22)^2]

matrix  =  {{-0.4, 0.2}, {0.1, 0.2}}
MatrixForm[%]
Eigenvalues[matrix]
N[%]

Clear[x, y]
num = 20;
x[0]  :=  4/9;
y[0]  :=  2/3;
x[n_]  :=  x[n]  = -0.4 x[n - 1] + 0.2  y[n - 1]  +  
          1  x[n - 1] y[n - 1]
y[n_]  :=  y[n]  = 
        0.1  x[n - 1]  +  0.2  y[n - 1] - 
          4  x[n - 1] y[n - 1]
data  =  
        Table[{n, x[n], y[n]}, {n, 0, 
            num}];

Solve[{x  == -0.4 x + 0.2  y   +  
            1  x y
, 
        y == 0.1  x  +  0.2  y - 4  x y}, {x, 
        y}]

phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
datax1  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 2]}, {n, 0, 
            num}];
datax2  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 3]}, {n, 0, 
            num}];
plotphase  =  
      ListPlot[phase, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"x1", "x2"}]
plotx1  =  
      ListPlot[datax1, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"n", "x1"}]
plotx2  =  
      ListPlot[datax2, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"n", "x2"}]

<< Graphics`PlotField`

f1  =  -0.4 x + 0.2  y  +  1  x y  - 
        x
f2  = 0.1  x  +  0.2  y - 4  x y  -  y

field = 
      PlotVectorField[{f1, f2}, {x, -1, 1}, {y, -1, 3}, 
        Frame -> True, 
        FrameLabel -> {"x1", "x2"}]
Show[field, plotphase, AspectRatio -> 1]

Clear[x, y]
num = 6;
x[0]  :=  5/9;
y[0]  :=  2/3;
x[n_]  :=  x[n]  = -0.4 x[n - 1] + 0.2  y[n - 1]  +  
          1  x[n - 1] y[n - 1]
y[n_]  :=  y[n]  = 
        0.1  x[n - 1]  +  0.2  y[n - 1] - 
          4  x[n - 1] y[n - 1]
data  =  
        Table[{n, x[n], y[n]}, {n, 0, 
            num}];

Solve[{x  == -0.4 x + 0.2  y   +  
            1  x y
, 
        y == 0.1  x  +  0.2  y - 4  x y}, {x, 
        y}]

phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
datax1  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 2]}, {n, 0, 
            num}];
datax2  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 3]}, {n, 0, 
            num}];
plotphase  =  
      ListPlot[phase, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"x1", "x2"}]
plotphase  =  
      ListPlot[datax1, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"n", "x1"}]
plotphase  =  
      ListPlot[datax2, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel -> {"n", "x2"}]

<< Graphics`PlotField`

f1  =  -0.4 x + 0.2  y  +  1  x y  - 
        x
f2  = 0.1  x  +  0.2  y - 4  x y  -  y

field = 
      PlotVectorField[{f1, f2}, {x, -1, 1}, {y, -1, 1}, 
        Frame -> True, 
        FrameLabel -> {"x1", "x2"}]

-722/171 + 1

-81/171 + 1

Clear[x]
RSolve[{x[n] ==  λ x[n - 1]  +  x[n - 1]}, x[n], 
      n]

matrix  =  {{-0.4, 0.2}, {α, 
          0.2}}
MatrixForm[%]
Eigenvalues[matrix]
N[%]
TeXForm[
      (1/2) (-0.2 - 
            0.89 (Sqrt[0.45 +  α]))]

Plot[{(1/2) (-0.2 - 
              0.89 (Sqrt[0.45 +  α])), 
        (1/2) (-0.2 + 
              0.89 (Sqrt[0.45 +  α])), 
        1, -1}, {α, -0.45, 10}, Frame -> True, 
      FrameLabel -> {"alpha", "eigenvalues"}]

Solve[
      (1/2) (-0.2 - 
              0.89 (Sqrt[0.45 +  α])) == -1, α]