Discrete and Continous Simulation

<< Graphics`PlotField`

The System

Here, we find any fixed points:

f1[x_, y_] = y^2 - 1  -  x;
f2[x_, y_] = x^2 - 1  -  y;
Solve[{f1[x, y] == 0, f2[x, y] == 0}, {x, 
          y}];
N[%]

And plot the vector field, from which it looks like the (0, -1) is a stable point:

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

The Discrete Model

We can iterate to get a sequence of points in the phase space:

Clear[x1, x2]
x1[0] :=  α1
x2[0] :=  α2
x1[n_] :=  x1[n]  =  
        f1[x1[n - 1], x2[n - 1]] + x1[n - 1]
x2[n_] :=  x2[n]  =  
        f2[x1[n - 1], x2[n - 1]] + x2[n - 1]
num  =  10;
data  =  Table[{n, x1[n], x2[n]}, {n, 0, num}];
data  =  data  /.  α1 -> 
            0.5  /. α2 -> -1.0
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}];
ListPlot[phase, PlotJoined -> True, Frame -> True, 
      FrameLabel -> {"x1", "x2"}]
ListPlot[datax1, PlotJoined -> True, Frame -> True, 
      FrameLabel -> {"n", "x1"}]
ListPlot[datax2, PlotJoined -> True, Frame -> True, 
      FrameLabel -> {"n", "x2"}]

The Continuous Model

Clear[x1, x2]
x1[0] :=  α1
x2[0] :=  α2
num = 10;
to  =  0;
tf  =  10;
h  =  N[(tf - to)/num]
x1[n_] :=  x1[n]  =  
        h f1[x1[n - 1], x2[n - 1]] + x1[n - 1]
x2[n_] :=  x2[n]  =  
        h f2[x1[n - 1], x2[n - 1]] + x2[n - 1]
data  =  Table[{h*n, x1[n], x2[n]}, {n, 0, num}];
data  =  data  /.  α1 -> 
            0.5  /. α2 -> -1.0
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}];
discretephase = 
      ListPlot[phase, PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"x1", "x2"}]
ListPlot[datax1, PlotJoined -> True, Frame -> True, 
      FrameLabel -> {"n", "x1"}]
ListPlot[datax2, PlotJoined -> True, Frame -> True, 
      FrameLabel -> {"n", "x2"}]

Clear[x1, x2]
x1[0] :=  0.5
x2[0] :=  -1.0
num = 500;
to  =  0;
tf  =  10;
h  =  N[(tf - to)/num]
x1[n_] :=  x1[n]  =  
        h f1[x1[n - 1], x2[n - 1]] + x1[n - 1]
x2[n_] :=  x2[n]  =  
        h f2[x1[n - 1], x2[n - 1]] + x2[n - 1]
data  =  Table[{h*n, x1[n], x2[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}];
contphase = 
      ListPlot[phase, PlotRange -> All, PlotJoined -> True, 
        Frame -> True, 
        FrameLabel -> {"x1", "x2"}]
ListPlot[datax1, PlotRange -> All, PlotJoined -> True, 
      Frame -> True, 
      FrameLabel -> {"Time, t", "x1"}]
ListPlot[datax2, PlotRange -> All, PlotJoined -> True, 
      Frame -> True, 
      FrameLabel -> {"Time, t", "x2"}]

Animation of the Discrete to the Continuous

I have set things up so that num=10 corresponds to the discrete case. Higher values of num will approach the continuous case.

\(Do[{
Clear[x1, x2], 
x1[0] := \ 
          0.5, 
x2[
            0] := \ \(-1.0\), 

x1[
            n_] := \ \(x1[n]\  = \ h\ f1[x1[n - 1], x2[n - 1]] + x1[n - 1]\), 
        x2[n_] := \ \(x2[n]\  = \ 
            h\ f2[x1[n - 1], x2[n - 1]] + 
              x2[n - 1]\), 

num = 
          i, 
to\  = \ 0, 
tf\  = \ 
          10, 
h\  = \ 
          N[\((tf - to)\)/num], 
data\  = \ 
          Table[{h*n, x1[n], x2[n]}, {n, 0, 
              num}], 
dataphase\  = \ 
          Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
              num}], 
ListPlot[dataphase, 
          PlotRange -> {{\(-0.5\), 0.5}, {\(-1\), \(-0.75\)}}, 
          PlotJoined -> True, Frame -> True, 
          FrameLabel -> {"Time, t", "x1"}]}, {i, 10, 15, 
        1}]\)

\(Do[{
Clear[x1, x2], 
x1[0] := \ 
          0.5, 
x2[
            0] := \ \(-1.0\), 

x1[
            n_] := \ \(x1[n]\  = \ h\ f1[x1[n - 1], x2[n - 1]] + x1[n - 1]\), 
        x2[n_] := \ \(x2[n]\  = \ 
            h\ f2[x1[n - 1], x2[n - 1]] + 
              x2[n - 1]\), 

num = 
          i, 
to\  = \ 0, 
tf\  = \ 
          10, 
h\  = \ 
          N[\((tf - to)\)/num], 
data\  = \ 
          Table[{h*n, x1[n], x2[n]}, {n, 0, 
              num}], 
dataphase\  = \ 
          Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
              num}], 
ListPlot[dataphase, 
          PlotRange -> {{\(-0.5\), 0.5}, {\(-1\), \(-0.75\)}}, 
          PlotJoined -> True, Frame -> True, 
          FrameLabel -> {"Time, t", "x1"}]}, {i, 10, 100, 
        5}]\)

Show[contphase, discretephase, field, AspectRatio -> 1]

An Animated Vector Field

Let's examine how the solutoin changes when we vary the parameter a:

f1[x_, y_] = y^2 - 1  -  x;
f2[x_, y_, a_] = a*x^2 - 1  -  y;
Solve[{f1[x, y] == 0, f2[x, y] == 0}, {x, y}]

The fixed points are still the same, and do not depend on a.

fp1  =  
        Graphics[Disk[{-1, 0}, 0.05]];
fp2  =  
        Graphics[Disk[{0, -1}, 0.05]];
fp3  =  
        Graphics[
          Disk[{(1/2) - ((Sqrt[5])/2), (1/2) - ((Sqrt[5])/2)}, 
            0.05]];
fp4  =  
        Graphics[Disk[{(1/2) + ((Sqrt[5])/2), (1/2) + ((Sqrt[5])/2)}, 0.05]];

step  =  0.2
Clear[field]
Table[
        field[i_]  =  
          Show[PlotVectorField[{f1[x, y], f2[x, y, i]}, {x, -3, 
                3}, {y, -3, 3}, 
              PlotRange -> {{-3, 3}, {-3, 3}}, 
              Frame -> True, FrameLabel -> {"x1", "x2"}, 
              AspectRatio -> 1], fp1, fp2, fp3, fp4], {i, 0, 4, 
          step}];

Problem 6.5.1

We have added weather conditions to the battle simulation, with weather only affecting the effectiveness of direct fire on both sides. w=0 corresponds to horrible weather, w=1 corresponds to good weather. In good weather, Red will win after 120 hours with 3.00 divisions remaining. Will bad weather improve the performance of Blue?

λ  =  3
w  =  1
f1[x_, 
          y_] = -w λ (0.05) y  - λ (0.005) x y;
f2[x_, y_] = -w 0.05  x  -  
          0.005  x y;

Clear[x1, x2]
x1[0] :=  5;
x2[0] :=  2;
num = 20;
x1[n_] :=  x1[n]  =   f1[x1[n - 1], x2[n - 1]] + 
          x1[n - 1]
x2[n_] :=  x2[n]  =   f2[x1[n - 1], x2[n - 1]] + 
          x2[n - 1]
data  =  Table[{12*n, x1[n], x2[n]}, {n, 0, num}];
\(TableForm[data, 
      TableHeadings -> {None, {"Time (hours)", "Red Divisions", "Blue Divisions"}}]\)
phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
contphase = 
      ListPlot[phase, PlotRange -> {{0, 5}, {0, 2}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Divisions", "Blue Divisions"}]

In good weather, Red wins the battle after 120 hours, with 3.00 divisions remaining.

λ  =  1
w  =  1
f1[x_, 
          y_] = -w λ (0.05) y  - λ (0.005) x y;
f2[x_, y_] = -w 0.05  x  -  
          0.005  x y;

Clear[x1, x2]
x1[0] :=  5;
x2[0] :=  2;
num = 25;
x1[n_] :=  x1[n]  =   f1[x1[n - 1], x2[n - 1]] + 
          x1[n - 1]
x2[n_] :=  x2[n]  =   f2[x1[n - 1], x2[n - 1]] + 
          x2[n - 1]
data  =  Table[{n, x1[n], x2[n]}, {n, 0, num}];
\(TableForm[data, 
      TableHeadings -> {None, {"Time (hours)", "Red Divisions", "Blue Divisions"}}]\)
phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
contphase = 
      ListPlot[phase, PlotRange -> {{0, 5}, {0, 2}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Divisions", "Blue Divisions"}]

\(Do[
{Clear[x1, x2, f1, 
          f2], 
f1[x_, 
            y_] = \(-w\)\ λ \((0.05)\)\ y\  - λ\ \((0.005)\)\ x\ y, 
f2[x_, y_] = \(-w\)\ 0.05  x\  - \ 
            0.005  x\ y, 
x1[0] := \ 
          5, 
x2[0] := \ 2, 
num = 
          70, 
x1[
            n_] := \ \(x1[n]\  = \ \ f1[x1[n - 1], x2[n - 1]] + x1[n - 1]\), 
        x2[n_] := \ \(x2[n]\  = \ \ f2[x1[n - 1], x2[n - 1]] + x2[n - 1]\), 
        data\  = \ 
          Table[{n, x1[n], x2[n]}, {n, 0, 
              num}], 

Do[{If[
              Part[data, n + 1, 2] < 0 || \ 
                Part[data, n + 1, 3] < 0, {Print["Battle ends after ", 
                  Part[data, n + 1, 1], " hours. w = ", 
                  w, " lambda = ", λ], Break[]}]}, {n, 0, 
            num}], 

phase\  = \ 
          Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, num}], 
        contphase = 
          ListPlot[phase, PlotRange -> {{0, 5}, {0, 2}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> {{0, 1, 2, 3, 4, 5}, {0, 1, 
                  2}, {}, {}}]}, {λ, 1, 5}, {w, 0.2, 1, 
        0.2}]\)

Continuous version of the battle problem

\(Do[
{Clear[x1, x2, f1, 
          f2], 
f1[x_, 
            y_] = \(-w\)\ λ \((0.05)\)\ y\  - λ\ \((0.005)\)\ x\ y, 
f2[x_, y_] = \(-w\)\ 0.05  x\  - \ 
            0.005  x\ y, 
x1[0] := \ 
          5, 
x2[0] := \ 2, 
num = 
          7000, 
h\  = \ 
          0.01, 
x1[
            n_] := \ \(x1[n]\  = \ \ h\ f1[x1[n - 1], x2[n - 1]] + 
              x1[n - 1]\), 
        x2[n_] := \ \(x2[n]\  = \ \ h\ f2[x1[n - 1], x2[n - 1]] + 
              x2[n - 1]\), 
        data\  = \ 
          Table[{h*n, x1[n], x2[n]}, {n, 0, 
              num}], 

Do[{If[
              Part[data, n + 1, 2] < 0 || \ 
                Part[data, n + 1, 3] < 0, {Print["Battle ends after ", 
                  Part[data, n + 1, 1], " hours. w = ", 
                  w, " lambda = ", λ], Break[]}]}, {n, 0, 
            num}], 

phase\  = \ 
          Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, num}], 
        contphase = 
          ListPlot[phase, PlotRange -> {{0, 5}, {0, 2}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> {{0, 1, 2, 3, 4, 5}, {0, 1, 
                  2}, {}, {}}]}, {λ, 1, 5}, {w, 0.2, 1, 
        0.2}]\)