Chapter 5

Question 1

Needs["DiscreteMath`RSolve`"]
Needs["Graphics`Graphics3D`"]
Needs["Graphics`PlotField3D`"]

The two equations are:
	(Subscript[x, t + 1]) = -8 - (Subscript[x, t]) + (Subscript[y, t])(Subscript[y, t + 1]) = 4 - 0.3  (Subscript[x, t]) + 0.9  (Subscript[y, t])

with difference equations:
	(Subscript[Δx, t])=-8 - 2  (Subscript[x, t]) + (Subscript[y, t])(Subscript[Δy, t])=4 - 0.3  (Subscript[x, t]) - 0.1  (Subscript[y, t])

Solve[{-8 - 2  xstar + ystar == 0, 
        4 - 0.3  xstar - 0.1  ystar == 0}, {xstar, ystar}]

matrixA=(({{-1, 1}, {-0.3, 0.9}}))

Eigenvalues[matrixA]

Eigenvectors[matrixA]

It will be noted that although the eigenvalues are those given in the text, it appears that the eigenvectors are different. But this is not the case. Consider the eigenvector (v^r) associated with the eignenvalue r = 0.726209. This is given in the programme by,
	(v^r)=(({{-0.501268}, {-0.865292}}))
But in the text we arbitrarily set (Subscript[v, 2]) = 1. We therefore must divide each element of the vector by -0.865292. Since,

-0.501268/-0.865292

then the eigenvector is,
	(v^r)=(({{0.579305}, {1}}))
which is the eigenvector derived in the text.

V=Transpose[Eigenvectors[matrixA]]

V1 = Inverse[V]

u0={-4.4, -12.8}

matrixD=(({{(-0.8262087348)^t, 0}, {0, 0.7262087348^t}}))

u = V . matrixD . V1 . u0

u1 = u[[1]]

u2 = u[[2]]

xt = 6.4 + u1

yt = 20.8 + u2

points = Table[{xt, yt}, {t, 0, 20}];

traj = 
        ListPlot[points,  PlotJoined -> True,  
          DisplayFunction -> Identity];

isoclines = 
        Plot[{8 + 2  x,  40 - 3  x}, {x, -2, 10},  
          DisplayFunction -> Identity];

Show[{traj, isoclines},  AspectRatio -> 1,  
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction];

This figure clearly shows the trajectory oscillating across the phase line, but converging on the fixed point, which confirms the stability established in Section 5.3.

Question 2

Example 5.10 is given by the following equations
	(Subscript[x, t + 1]) = -0.85078 (Subscript[x, t]) - (Subscript[y, t])
(Subscript[y, t + 1]) = (Subscript[x, t]) + 2.35078  (Subscript[y, t])
which has a fixed point at the origin. Subtracting (Subscript[x, t]) from both sides of the first equation and (Subscript[y, t]) from both sides of the second equation we obtain,
	(Subscript[Δx, t])=(Subscript[x, t + 1]) - (Subscript[x, t]) = -1.85078 (Subscript[x, t]) - (Subscript[y, t])
(Subscript[Δy, t])=(Subscript[y, t + 1]) - (Subscript[y, t]) = (Subscript[x, t]) + 1.35078  (Subscript[y, t])
The two phase lines are found by setting Δ(Subscript[x, t]) = 0 and Δ(Subscript[y, t]) = 0 respectively, i.e.,

Solve[-1.85078 x - y == 0, y]

Solve[x + 1.35078  y == 0, y]

It should be noted that both phase lines have a negative slope.

The vector forces can be established by noting:
	If (Subscript[Δx, t])>0 then (Subscript[y, t]) < -1.85078 (Subscript[x, t]), i.e. below Δ(Subscript[x, t]) = 0, x is rising
	If (Subscript[Δx, t])<0 then (Subscript[y, t]) > -1.85078 (Subscript[x, t]), i.e. above Δ(Subscript[x, t]) = 0, x is falling
and
	If (Subscript[Δy, t])>0 then (Subscript[y, t]) > -0.740313 (Subscript[x, t]), i.e. above Δ(Subscript[x, t]) = 0, y is rising
	If (Subscript[Δy, t])<0 then (Subscript[y, t]) < -0.740313 (Subscript[x, t]), i.e. below Δ(Subscript[x, t]) = 0, y is falling

isoclines = 
        Plot[{-1.85078 x, -0.740313 x},  {x, -3, 
            3}];

matrixA=(({{-0.85078, -1}, {1, 2.35078}}))

Eigenvalues[matrixA]

Eigenvectors[matrixA]

To find the saddle-path equations

Solve[0.331007  x == -0.943628 y, y]

Solve[-0.943628 x == 0.331007  y, y]

saddles = 
        Plot[{-0.350781 x, -2.85078 x}, {x, -3, 3}, 
          PlotStyle -> Dashing[{0.02}] ];

Show[isoclines, saddles,  
        AxesLabel -> {"x", "y"}];

Question 3

The two equations of Example 5.4 are:
	(Subscript[x, t + 1]) = -8 - (Subscript[x, t]) + (Subscript[y, t])(Subscript[y, t + 1]) = 4 - 0.3  (Subscript[x, t]) + 0.9  (Subscript[y, t])
and although the question refers to setting up the problem on a spreadsheet, we can also use Mathematica to show a table of values. We do this by solving the system of equations for each of the initial conditions and set the results out as a table.

Clear[x, y]

For (3,10)

sol1 = 
      RSolve[{x[t + 1] == -8 - x[t] + y[t], 
          y[t + 1] == 4 - 0.3  x[t] + 0.9  y[t], x[0] == 3, 
          y[0] == 10}, {x[t], y[t]}, t]

path1 = 
      Table[{t, sol1[[1, 1, 2]], sol1[[1, 2, 2]]}, {t, 0, 20}] ; 
    TableForm[path1, 
      TableHeadings -> {{ }, {"t", "x[t]", "y[t]"}}]

Clear[x, y]

For (3,30)

sol2 = 
      RSolve[{x[t + 1] == -8 - x[t] + y[t], 
          y[t + 1] == 4 - 0.3  x[t] + 0.9  y[t], x[0] == 3, 
          y[0] == 30}, {x[t], y[t]}, t]

path2 = 
      Table[{t, sol2[[1, 1, 2]], sol2[[1, 2, 2]]}, {t, 0, 20}] ; 
    TableForm[path2, 
      TableHeadings -> {{ }, {"t", "x[t]", "y[t]"}}]

Clear[x, y]

For (10,10)

sol3 = 
      RSolve[{x[t + 1] == -8 - x[t] + y[t], 
          y[t + 1] == 4 - 0.3  x[t] + 0.9  y[t], x[0] == 10, 
          y[0] == 10}, {x[t], y[t]}, t]

path3 = 
      Table[{t, sol3[[1, 1, 2]], sol3[[1, 2, 2]]}, {t, 0, 20}] ; 
    TableForm[path3, 
      TableHeadings -> {{ }, {"t", "x[t]", "y[t]"}}]

Clear[x, y]

For (10,30)

sol4 = 
      RSolve[{x[t + 1] == -8 - x[t] + y[t], 
          y[t + 1] == 4 - 0.3  x[t] + 0.9  y[t], x[0] == 10, 
          y[0] == 30}, {x[t], y[t]}, t]

path4 = 
      Table[{t, sol4[[1, 1, 2]], sol4[[1, 2, 2]]}, {t, 0, 20}] ; 
    TableForm[path4, 
      TableHeadings -> {{ }, {"t", "x[t]", "y[t]"}}]

Question 4

The system for example 5.5 is
	\(TraditionalForm
      \`x\_t = x\_\(t - 1\) + 2  y\_\(t - 1\) + z\_\(t - 1\)\)(Subscript[y, t]) = -(Subscript[x, t - 1]) + (Subscript[y, t - 1])\(TraditionalForm
      \`z\_t = 3  x\_\(t - 1\) - 6  y\_\(t - 1\) - z\_\(t - 1\)\)

(a) We can first use Mathematica to check the solution in the text for this set of equations for the initial conditions x(0) = 3,  
      y(0) = -4 and z(0) = 3.

Clear[x, y, z, t]

\(sol = 
      RSolve[{x[t] == x[t - 1] + 2  y[t - 1] + z[t - 1], 
          y[t] == \(-x[t - 1]\) + y[t - 1], 
		z[t] == 
            3  x[t - 1] - 6  y[t - 1] - z[t - 1], x[0] == 3, y[0] == \(-4\), 
          z[0] == 3}, {x[t], y[t], z[t]}, t]\)

Note that this is a fuller answer than in the text. Here we have a qualification for the initial condition that - 5 must be added to x(t) to obtain the correct answer. Similarly for y(t) and z(t) -5 and +15 must be added respectively for the initial condition to be correct.

Clear[x, y, z]

(b)

\(path = 
      Table[{t, sol[\([1, 1, 2]\)], sol[\([1, 2, 2]\)], 
          sol[\([1, 3, 2]\)]}, {t, 0, 20}]\ ; 
    TableForm[path, 
      TableHeadings -> {{\ }, {"t", "x[t]", "y[t]", "z[t]"}}]\)

which is explosive. Plotting each against time we have,

listx = Table[sol[[1, 1, 2]], {t, 0, 20}];

listy = Table[sol[[1, 2, 2]], {t, 0, 20}];

listz = Table[sol[[1, 3, 2]], {t, 0, 20}];

linex = 
        ListPlot[listx, PlotJoined -> True,  
          AxesLabel -> {"t", "x(t)"}];

liney = 
        ListPlot[listy, PlotJoined -> True, 
          AxesLabel -> {"t", "y(t)"},  
          PlotStyle -> Dashing[{0.01}]];

linez = 
        ListPlot[listz, PlotJoined -> True, 
          AxesLabel -> {"t", "z(t)"},  
          PlotStyle -> Dashing[{0.02}]];

Show[linex, liney, linez,  
        AxesLabel -> {"t", "x,y,z"}];

(c)

pointsxy = 
        Table[{sol[[1, 1, 2]], sol[[1, 2, 2]]}, {t, 0, 
            20}];

ListPlot[pointsxy,  PlotJoined -> True, 
        Ticks -> {{100000, 300000, 500000}, Automatic}, 
        AxesLabel -> {"x", "y"}];

(d)

pointsxz = 
        Table[{sol[[1, 1, 2]], sol[[1, 3, 2]]}, {t, 0, 
            20}];

ListPlot[pointsxz,  PlotJoined -> True, 
        Ticks -> {{100000, 300000, 500000}, Automatic}, 
        AxesLabel -> {"x", "z"}];

(e)

pointsyz = 
        Table[{sol[[1, 2, 2]], sol[[1, 3, 2]]}, {t, 0, 
            20}];

ListPlot[pointsyz,  PlotJoined -> True, 
        Ticks -> {{-100000, -300000, -500000}, Automatic}, 
        AxesLabel -> {"y", "z"}];

Question 5

We already have the values for x[t],y[t] and z[t] in Question 4. These are specified as sol[[1,1,2]] for x, sol[[1,2,2]] for y and sol[[1,3,2]] for z. However, to plot such a set of points in 3-dimensional space, we need to utilise the ScatterPlot3D command which is available in the Graphics3D package. This package must first be loaded. Since the figures become large (both negatively and positively) we plot only the first 10 points.

points = 
        Table[{sol[[1, 1, 2]], sol[[1, 2, 2]], 
            sol[[1, 3, 2]]}, {t, 0, 10}];

ScatterPlot3D[points,  AxesLabel -> {x, y, z}, 
        Ticks -> {{0, 1000, 2000}, {0, -1000, -2000}, Automatic}, 
        PlotStyle -> PointSize[0.03], AspectRatio -> 1];

Question 6

Example 5.5 can be written first as a set of difference equations:
	(Subscript[Δx, t])=2  (Subscript[y, t]) + (Subscript[z, t])(Subscript[Δy, t])=-(Subscript[x, t])
(Subscript[Δz, t])=3  (Subscript[x, t]) - 6  (Subscript[y, t]) - 2  (Subscript[z, t])
We first need to load the PlotField3D command. The 3-dimensional direction field can then be obtained using the following instructions.

\(\(PlotVectorField3D[{2  y + z, \(-x\), 3  x - 6  y - 2  z}, 
	{x, 0, 
          2000}, {y, 0, \(-2000\)}, {z, 0, 6000}, PlotPoints -> 5, 
        VectorHeads -> True, ScaleFunction -> \((1 &)\), 
        AspectRatio -> 0.8];\)\)

But this does not appear to give us much additional insight.

Question 7

The matrix for each system is respectively:
	A1=(({{0, 1}, {-1, 0}}))		A2 = (({{-2, 1}, {1, -2}}))	A3 = (({{3, -4}, {1, -2}}))
If r and s are the eigenvalues of the system, then the diagonal matrix, D, is given by:
	D = (({{r, 0}, {0, s}}))

(i) System 1

matrixA1=(({{0, 1}, {-1, 0}}))

(a)

Eigenvalues[matrixA1]

(b)

Eigenvectors[matrixA1]

(c)

The diagonal matrix for this first system is
	D1 = (({{-I, 0}, {0, I}}))

(ii) System 2

matrixA2=(({{-2, 1}, {1, -2}}))

(a)

Eigenvalues[matrixA2]

(b)

Eigenvectors[matrixA2]

(c)

The diagonal matrix for this first system is
	D2 = (({{-3, 0}, {0, -1}}))

(iii) System 3

matrixA3=(({{3, -4}, {1, -2}}))

(a)

Eigenvalues[matrixA3]

(b)

Eigenvectors[matrixA3]

(c)

The diagonal matrix for this first system is
	D3 = (({{-1, 0}, {0, 2}}))