Chapter 4

Needs["Graphics`PlotField`"]

Question 1

(i) We can readily write this
	(dy/dx) = (y/(2  x))
	(dy/y) = (dx/(2  x))
Integrating both sides
	\(TraditionalForm
      \`∫\(dy\/y\) dy = 
        \(∫\(dx\/\(2  x\)\) dx
\(ln(y)\) = 
          \(1\/2\) \(ln(x)\) + c\_0\)\)
where (Subscript[c, 0]) is the constant of integration. Then
	y = c (Sqrt[ x])
If x(0) = 2 and y(0) = 3, then 3 = c (Sqrt[ 2]) and
	y = ((3 (Sqrt[ x]))/(Sqrt[2]))

(ii) To verify this result

DSolve[y'[x] == y[x]/(2  x), y[x], x]

DSolve[ {y'[x] == y[x]/(2  x), y[2] == 3}, y[x], x]

Question 2

arrows = 
        PlotVectorField[ {x - 3  y, -2 x + y},  {x, -20, 
            20},  {y, -20, 20}, Axes -> True,  
          AxesLabel -> {"x", "y"}, ScaleFunction -> (1 &), 
          AspectRatio -> 0.8];

(a)

sol1 = 
      NDSolve[ {x'[t] == 
            x[t] - 3  y[t], y'[t] == -2 x[t] + y[t],  x[0] == 4,  
          y[0] == 2}, {x, y}, {t, 0, 10}]

traj1 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol1],  {t, 0, 1},  
          DisplayFunction -> Identity];

(b)

sol2 = 
      NDSolve[ {x'[t] == 
            x[t] - 3  y[t], y'[t] == -2 x[t] + y[t],  x[0] == 4,  
          y[0] == 5}, {x, y}, {t, 0, 10}]

traj2 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol2],  {t, 0, 1},  
          DisplayFunction -> Identity];

(c)

sol3 = 
      NDSolve[ {x'[t] == 
            x[t] - 3  y[t], y'[t] == -2 x[t] + y[t],  
          x[0] == -4,  y[0] == -2}, {x, y}, {t, 0, 10}]

traj3 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol3],  {t, 0, 1},  
          DisplayFunction -> Identity];

(d)

sol4 = 
      NDSolve[ {x'[t] == 
            x[t] - 3  y[t], y'[t] == -2 x[t] + y[t],  
          x[0] == -4,  y[0] == 5}, {x, y}, {t, 0, 10}]

traj4 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol4],  {t, 0, 1},  
          DisplayFunction -> Identity];

Show[{arrows, traj1, traj2, traj3, traj4},  
        PlotRange -> {{-20, 20}, {-20, 20}}, 
        AspectRatio -> 0.8];

Question 3

Clear[x, y, sol1, sol2, sol3, sol4, traj1, traj2, traj3, traj4, 
      arrows]

(i)

arrows1 = 
        PlotVectorField[ {-3 x + y, x - 3  y},  {x, 0, 4}, {y, 0, 8},  
          Axes -> True, AxesLabel -> {"x", "y"}, 
          ScaleFunction -> (1 &), AspectRatio -> 1];

(i)

sol11 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == 4,  y[0] == 8}, {x, y}, {t, 0, 
          10}]

traj11 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol11],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

sol12 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == 4,  y[0] == 2}, {x, y}, {t, 0, 
          10}]

traj12 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol12],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

lines1 = 
        Plot[{3  x, (1/3) x}, {x, 0, 4}, 
          DisplayFunction -> Identity];

Show[{arrows1, traj11, traj12, lines1}, 
        PlotRange -> {{0, 4}, {0, 8}}];

(ii)

arrows2 = 
        PlotVectorField[ {-3 x + y, x - 3  y},  {x, 0, -4}, {y, 
            0, -8},  Axes -> True, AxesLabel -> {"x", "y"}, 
          ScaleFunction -> (1 &), AspectRatio -> 1];

sol21 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == -4,  y[0] == -8}, {x, y}, {t, 
          0, 10}]

traj21 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol21],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

sol22 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == -4,  y[0] == -2}, {x, y}, {t, 
          0, 10}]

traj22 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol22],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

lines2 = 
        Plot[{3  x, (1/3) x}, {x, 0, -4}, 
          DisplayFunction -> Identity];

Show[{arrows2, lines2, traj21, traj22}, 
        PlotRange -> {{0, -4}, {0, -8}}];

(iii)

arrows3 = 
        PlotVectorField[ {-3 x + y, x - 3  y},  {x, -3, 
            3}, {y, -10, 10},  Axes -> True, 
          AxesLabel -> {"x", "y"}, ScaleFunction -> (1 &), 
          AspectRatio -> 1];

sol31 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == 2,  y[0] == 10}, {x, y}, {t, 0, 
          10}]

traj31 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol31],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

lines3 = 
        Plot[{3  x, (1/3) x}, {x, -3, 3}, 
          DisplayFunction -> Identity];

sol32 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == -2,  y[0] == -10}, {x, y}, {t, 
          0, 10}]

traj32 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol32],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

Show[{arrows3, lines3, traj31, traj32}, 
        PlotRange -> {{-3, 3}, {-10, 10}}];

(iv)

arrows4 = 
        PlotVectorField[ {-3 x + y, x - 3  y},  {x, -1, 
            2}, {y, -5, 0},  Axes -> True, 
          AxesLabel -> {"x", "y"}, ScaleFunction -> (1 &), 
          AspectRatio -> 1];

sol4 = 
      NDSolve[ {x'[t] == -3 x[t] + y[t], y'[t] == 
            x[t] - 3  y[t],  x[0] == 2,  y[0] == -5}, {x, y}, {t, 0, 
          10}]

traj4 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol4],  {t, 0, 1}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

lines4 = 
        Plot[{3  x, (1/3) x}, {x, -1, 2}, 
          DisplayFunction -> Identity];

Show[{arrows4, lines4, traj4}, 
        PlotRange -> {{-1, 2}, {-5, 0}}];

As can be seen from the resulting figure, the trajectory passes into another quadrant before converging on the fixed point.

Question 4

Clear[x, y, sol1, sol2, sol3, sol4, sol5, traj1, traj2, traj3, traj4, 
      traj5, lines1, lines2, lines3, lines4, lines5, arrows1, arrows2, 
      arrows3, arrows4, arrows5]

Solve[{-2 x - y + 9 == 0, -y + x + 3 == 0}, {x, y}]

arrows1 = 
        PlotVectorField[{-2 x - y + 9, -y + x + 3}, {x, 0, 4}, {y, 0, 
            9},  AspectRatio -> 1, Axes -> True, 
          AxesLabel -> {"x", "y"},  
          ScaleFunction -> (1 &)];

(i)

sol1 = 
      NDSolve[ {x'[t] == -2 x[t] - y[t] + 9, y'[
              t] == -y[t] + x[t] + 3,  x[0] == 0,  y[0] == 2}, {x, 
          y}, {t, 0, 20}]

traj1 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol1],  {t, 0, 5}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

lines1 = 
        Plot[{-2 x + 9, x + 3}, {x, 0, 4}, 
          DisplayFunction -> Identity];

(ii)

sol2 = 
      NDSolve[ {x'[t] == -2 x[t] - y[t] + 9, y'[
              t] == -y[t] + x[t] + 3,  x[0] == 2,  y[0] == 8}, {x, 
          y}, {t, 0, 20}]

traj2 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol2],  {t, 0, 5}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

(iii)

sol3 = 
      NDSolve[ {x'[t] == -2 x[t] - y[t] + 9, y'[
              t] == -y[t] + x[t] + 3,  x[0] == 3,  y[0] == 1}, {x, 
          y}, {t, 0, 20}]

traj3 = 
        ParametricPlot[Evaluate[ {x[t], y[t]}  /.  sol3],  {t, 0, 5}, 
          PlotStyle -> {Thickness[0.007]},  
          DisplayFunction -> Identity];

Show[{arrows1, lines1, traj1, traj2, traj3}];

The resulting diagram shows quite clearly that all trajectories follow a counter-clockwise spiral towards the fixed point. However, it also shows that the spiral motion involves a direct convergence on the fixed point. In other words, there is no repeated over- and under-shooting.

Question 5

Clear[x, y]

(i) Characteristic roots obtained by defining the matrix of the system and using the Eigenvalues command.

a = {{2, 3}, {3, 2}} ;  MatrixForm[a]

Eigenvalues[a]

(ii) Now employ the Eigenvectors command.

Eigenvectors[a]

(iii) Solving directly we have:

DSolve[{x'[t] == 2  x[t] + 3  y[t], y'[t] == 
          3  x[t] + 2  y[t]}, {x[t], y[t]}, t]

(iv) The initial value problem is:

DSolve[{x'[t] == 2  x[t] + 3  y[t], y'[t] == 
          3  x[t] + 2  y[t],  x[0] == 1, y[0] == 0}, {x[t], y[t]}, 
      t]

Question 6

Clear[x, y, a]

(i)

a = {{1, 3}, {5, 3}};  MatrixForm[a]

Eigenvalues[a]

(ii)

Eigenvectors[a]

(iii)

DSolve[{x'[t] == x[t] + 3  y[t], y'[t] == 5  x[t] + 3  y[t],  
        x[0] == 1, y[0] == 3}, {x[t], y[t]}, t]

Question 7

Clear[x, y, a]

a = {{1, 1}, {-2, 4}};  MatrixForm[a]

Eigenvalues[a]

Eigenvectors[a]

The Wronksian matrix, here denoted w, is formed from the eigenvectors of the system. If Det[w] is nonzero then the eigenvectors are linearly independent.

w = Eigenvectors[a];  MatrixForm[w]

Det[w]

Question 8

Clear[x, y, a, w]

a = {{1, 0, 0}, {2, 3, 1}, {0, 2, 4}};  MatrixForm[a]

(i) The eigenvalues and eigenvectors are as follows:

Eigenvalues[a]

Eigenvectors[a]

(ii)The general solution is derived by solving the following differential equation system:

DSolve[ {x'[t] == x[t],  y'[t] == 
          2  x[t] + 3  y[t] + z[t],  z'[t] == 2  y[t] + 4  z[t]}, {x[t], 
        y[t], z[t]}, t]

(iii) The Wronksian is given by:

w = Eigenvectors[a];  MatrixForm[w]

Det[w]