Deterministic Battle Simulations

f1[x_, y_, w_] = -0.015 y ;
f2[x_, y_, w_] = -0.010 
          x ;

Clear[x1, x2]
x1[0] :=  3000;
x2[0] :=  3000;
\(\(num = 120;\)
\)
x1[n_] :=  x1[n] = 
        Round[f1[x1[n - 1], x2[n - 1], w[n - 1]] + x1[n - 1]]
x2[n_] :=  x2[n] = 
        Round[f2[x1[n - 1], x2[n - 1], w[n - 1]] + 
            x2[n - 1]]

data  =  Table[{0.1*n, x1[n], x2[n]}, {n, 0, num}];
TableForm[data, 
      TableHeadings -> {None, {"Time (hours)", "Red Soldiers", "\
Blue Soldiers"}}]
phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
contphase = 
      ListPlot[phase, PlotRange -> {{0, 3000}, {0, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]

Adding weather: Stochastic Effect

f1[x_, y_, w_] = -w 0.015 y ;
f2[x_, y_, w_] = - 0.010 x ;

    
Clear[x1, x2, w]
w[0]  := 1.0;
x1[0] :=  3000;
x2[0] :=  3000;
\(\(\(num = 200;\)\(
\)
    \)\)
w[n_]  :=  w[n]  =  
        Module[{}, {next = w[n - 1]  +  0.1*Random[Integer, {-1, 1}], 
            If[next ≥ 1,  next = 1.0], 
            If[next  ≤  0,  next = 0.0]}; 
          next]
x1[n_] :=  x1[n] = 
        Round[f1[x1[n - 1], x2[n - 1], w[n - 1]] + x1[n - 1]]
\(\(\(x2[n_] := \ \(x2[n] = 
        Round[f2[x1[n - 1], x2[n - 1], w[n - 1]] + x2[n - 1]]\)\)\(
\)
    \)\)
data  =  
        Table[{0.1*n, x1[n], x2[n], w[n]}, {n, 0, num}];
TableForm[data, 
      TableHeadings -> {None, {"Time (hours)", "Red Soldiers", "\
Blue Soldiers", "Weather"}}]
phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
\(\(\(contphase = 
      ListPlot[phase, PlotRange -> {{0, 3000}, {0, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]\)\(
\)
    \)\)
weather  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 4]}, {n, 0, 
            num}];
weatherplot = 
      ListPlot[weather, PlotRange -> All, PlotJoined -> True, 
        Frame -> True, 
        FrameLabel -> {"Time", "Weather"}]

Simulation run for 20 times.

f1[x_, y_, w_] = -w 0.015 y ;
f2[x_, y_, w_] = - 0.010 x ;

    
\(Do[
{Clear[x1, x2, w], 
w[0]\  := 
          1.0, 
x1[0] := \ 3000, 
x2[0] := \ 3000, 

        num = 200, 


        w[n_]\  := \ \(w[n]\  = \ 
            Module[{}, {next = 
                  w[n - 1]\  + \ 0.1*Random[Integer, {\(-1\), 1}], 
                If[next ≥ 1, \ next = 1.0], 
                If[next\  ≤ \ 0, \ next = 0.0]}; 
              next]\), 
x1[n_] := \ \(x1[n] = 
            Round[f1[x1[n - 1], x2[n - 1], w[n - 1]] + x1[n - 1]]\), 

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


        data\  = \ Table[{0.1*n, x1[n], x2[n], w[n]}, {n, 0, num}], 

        TableForm[data, 
          TableHeadings -> {None, {"Time (hours)", "Red Soldiers", "Blue Soldiers", "Weather"}}]; 
phase\  = \ Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, num}], 

        weather\  = \ 
          Table[{Part[data, n + 1, 1], Part[data, n + 1, 4]}, {n, 0, 
              num}], 

p1\  = \ 
          ListPlot[phase, PlotRange -> {{0, 3000}, {0, 3000}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> False, 
            FrameLabel -> {"Red", "Blue"}, 
            DisplayFunction -> Identity]; 
p2\  = \ 
          ListPlot[weather, PlotRange -> {{0, 20}, {0, 1}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> False, 
            FrameLabel -> {"Time", "Weather"}, 
            DisplayFunction -> 
              Identity]; 

Show[
          GraphicsArray[{p1, 
              p2}]]

}, {i, 1, 
        20}]\)

A more realistic model of weather

f1[x_, y_, w_] = -w 0.015 y ;
f2[x_, y_, w_] = - 0.010 x ;

    
Clear[x1, x2, w]
w[0]  := 0.8;
x1[0] :=  3000;
x2[0] :=  3000;
\(\(num = 200;\)
\)
w[n_]  :=  w[n]  =  
        Module[{}, {next = w[n - 1]  +  0.03*Random[Integer, {-1, 1}], 
            If[next ≥ 1,  next = 1.0], 
            If[next  ≤  0,  next = 0.0]}; 
          next]
x1[n_] :=  x1[n] = 
        Round[f1[x1[n - 1], x2[n - 1], w[n - 1]] + x1[n - 1]]
x2[n_] :=  x2[n] = 
        Round[f2[x1[n - 1], x2[n - 1], w[n - 1]] + 
            x2[n - 1]]

data  =  
        Table[{0.1*n, x1[n], x2[n], w[n]}, {n, 0, num}];
TableForm[data, 
      TableHeadings -> {None, {"Time (hours)", "Red Soldiers", "\
Blue Soldiers", "Weather"}}]
phase  =  
        Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, 
            num}];
\(contphase = 
      ListPlot[phase, PlotRange -> {{0, 3000}, {0, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]
\)
weather  =  
        Table[{Part[data, n + 1, 1], Part[data, n + 1, 4]}, {n, 0, 
            num}];
weatherplot = 
      ListPlot[weather, PlotRange -> {{0, 20}, {0, 1}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Time", "Weather"}]

f1[x_, y_, w_] = -w 0.015 y ;
f2[x_, y_, w_] = - 0.010 x ;

    
\(Do[
{Clear[x1, x2, w], 
w[0]\  := 
          0.8, 
x1[0] := \ 3000, 
x2[0] := \ 3000, 

        num = 300, 


        w[n_]\  := \ \(w[n]\  = \ 
            Module[{}, {next = 
                  w[n - 1]\  + \ 0.03*Random[Integer, {\(-1\), 1}], 
                If[next ≥ 1, \ next = 1.0], 
                If[next\  ≤ \ 0, \ next = 0.0]}; 
              next]\), 
x1[n_] := \ \(x1[n] = 
            Round[f1[x1[n - 1], x2[n - 1], w[n - 1]] + x1[n - 1]]\), 

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


        data\  = \ Table[{0.1*n, x1[n], x2[n], w[n]}, {n, 0, num}], 

        TableForm[data, 
          TableHeadings -> {None, {"Time (hours)", "Red Soldiers", "Blue Soldiers", "Weather"}}]; 
phase\  = \ Table[{Part[data, n + 1, 2], Part[data, n + 1, 3]}, {n, 0, num}], 

        weather\  = \ 
          Table[{Part[data, n + 1, 1], Part[data, n + 1, 4]}, {n, 0, 
              num}], 

p1\  = \ 
          ListPlot[phase, PlotRange -> {{0, 3000}, {0, 3000}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> False, 
            FrameLabel -> {"Red", "Blue"}, 
            DisplayFunction -> Identity]; 
p2\  = \ 
          ListPlot[weather, PlotRange -> {{0, 20}, {0, 1}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> False, 
            FrameLabel -> {"Time", "Weather"}, 
            DisplayFunction -> 
              Identity]; 

Show[
          GraphicsArray[{p1, 
              p2}]]

}, {i, 1, 
        20}]\)

A Stochastic Model

Clear[x, y, data]
Probability[x_, y_]  =  
        0.15 y/(0.15 y  +  
                0.1 x);

x[0] := 3000
y[0] := 3000
data  =  {};

n = 0;
While[
      x[n] > 0 && y[n] > 0, 
n = 
        n + 1; 
death  =  
        Random[Real, {0, 1}]; 
If[
        death  < Probability[x[n - 1], y[n - 1]] , {y[n]  =  y[n - 1], 
          x[n]  =  x[n - 1] - 1}, {y[n]  =  y[n - 1] - 1, 
          x[n]  =  x[n - 1]}]; 
data  =  
        Append[data, {x[n], 
            y[n]}]
]

Last[data]

\(p1\  = \ 
      ListPlot[data, PlotRange -> {{0, 3000}, {0, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]
\)
p2  =  
      ListPlot[data, PlotRange -> {{2950, 3000}, {2950, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]

Here is 100 runs of the stochastic model, to see how many Blue soldiers are left 
after the battle. We need to load the Graphics package to do histograms:

<< Graphics`Graphics`

Clear[data2]
data2 = {};
Do[

{Clear[x, y, 
          data], 
Probability[x_, y_]  =  
          0.15 y/(0.15 y  +  
                  0.1 x), 

x[0] := 
          3000, 
y[0] := 
          3000, 
data  =  {}, 

n = 0, 
While[
          x[n] > 0 && y[n] > 0, 
n = 
            n + 1; 
death  =  
            Random[Real, {0, 1}]; 
If[
            death  < 
              Probability[x[n - 1], y[n - 1]] , {y[n]  =  y[n - 1], 
              x[n]  =  x[n - 1] - 1}, {y[n]  =  y[n - 1] - 1, 
              x[n]  =  x[n - 1]}]; 
data  =  
            Append[data, {x[n], 
                y[n]}]
], 

final  =  
          Part[Last[data], 2], 
data2 = 
          Append[data2, final]}, 
{i, 1, 100}]

data2

histy1 = 
      Histogram[data2, HistogramCategories -> 5, Frame -> True, 
        FrameLabel -> {"Final Blue Soldiers", "Frequency"}]

\(histy2 = 
      Histogram[data2, HistogramCategories -> 5, Frame -> True, 
        HistogramScale -> 1, 
        FrameLabel -> {"Final Blue Soldiers", "Probability"}]\)

s  =  50
fit2  =  
      Plot[1/Sqrt[2  Pi] /s*Exp[-(x - 1725)^2/2/s^2], {x, 
          1550, 1900}]
Show[histy2, fit2]

Include time in the stochastic model

Clear[x, y, t, data]
Probability[x_, y_]  =  
        0.15 y/(0.15 y  +  0.1 x);
Probabilitytime[x_, y_]  =  
        1.0/(0.15 y  +  
              0.1 x);

x[0] := 3000
y[0] := 3000
t[0]  :=  0
data  =  {};

n = 0;
While[
      x[n] > 0 && y[n] > 0, 
n = 
        n + 1; 
death  =  
        Random[Real, {0, 1}]; 
t[n]  =  
        t[n - 1]  - 
          Log[death]*
            Probabilitytime[x[n - 1], 
              y[n - 1]]; 
nfinal  =  
        n; 
If[
        death  < Probability[x[n - 1], y[n - 1]] , {y[n]  =  y[n - 1], 
          x[n]  =  x[n - 1] - 1}, {y[n]  =  y[n - 1] - 1, 
          x[n]  =  x[n - 1]}]; 
data  =  
        Append[data, {x[n], 
            y[n]}]
]

t[nfinal]

\(p1\  = \ 
      ListPlot[data, PlotRange -> {{0, 3000}, {0, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]
\)
p2  =  
      ListPlot[data, PlotRange -> {{2950, 3000}, {2950, 3000}}, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel -> {"Red Soldiers", "Blue Soldiers"}]

Include time,  and do many runs and get a histogram for the duration of the battle.

Clear[data2]
data2 = {};
Do[

{Clear[x, y, t, 
          data], 
Probability[x_, y_]  =  
          0.15 y/(0.15 y  +  
                  0.1 x), 
Probabilitytime[x_, y_]  =  
          1.0/(0.15 y  +  
                0.1 x), 

x[0] := 
          3000, 
y[0] := 3000, 
t[0] := 
          0, 
data  =  {}, 

n = 0, 
While[
          x[n] > 0 && y[n] > 0, 
n = 
            n + 1; 
nfinal  =  
            n; 
death  =  
            Random[Real, {0, 1}]; 
t[n]  =  
            t[n - 1]  - 
              Log[death]*
                Probabilitytime[x[n - 1], y[n - 1]]; 
If[
            death  < 
              Probability[x[n - 1], y[n - 1]] , {y[n]  =  y[n - 1], 
              x[n]  =  x[n - 1] - 1}, {y[n]  =  y[n - 1] - 1, 
              x[n]  =  x[n - 1]}]; 
data  =  
            Append[data, {x[n], 
                y[n]}]
], 

data2 = Append[data2, t[nfinal]]}, 
{i, 
        1, 100}]

data2

What is the mean duration of the battle:

Sum[Part[data2, n], {n, 1, 100}]/100

histy1 = 
      Histogram[data2, HistogramCategories -> 8, Frame -> True, 
        FrameLabel -> {"Duration of Battle", "Frequency"}]

A Stochastic Model with weather

Clear[data2]
data2 = {};
\(Do[

{Clear[x, y, w, t, nfinal, 
          data, weather], 
Probability[x_, y_, w_]\  = \ 
          w*0.15\ y/\((w*0.15\ y\  + \ 
                  0.1\ x)\), 
Probabilitytime[x_, y_, 
            w_]\  = \ 
          1.0/\((w*0.15\ y\  + \ 
                0.1\ x)\), 

x[0] := 
          3000, 
y[0] := 3000, 
t[0] := 
          0, 
w[0] := 
          0.8, 
data\  = \ {}, 
weather = {}, 

n = 0, 
While[
          x[n] > 0 && y[n] > 0, 
n = 
            n + 1; 
nfinal\  = \ 
            n; 
death\  = \ 
            Random[Real, {0, 1}]; 
t[n]\  = \ 
            t[n - 1]\  - 
              Log[death]*
                Probabilitytime[x[n - 1], y[n - 1], 
                  w[n - 1]]; 
w[n]\  = \ 
            Module[{}, {next = 
                  w[n - 1]\  + \ 0.005*Random[Integer, {\(-1\), 1}], 
                If[next ≥ 1, \ next = 1.0], 
                If[next\  ≤ \ 0, \ next = 0.0]}; 
              next]; 
If[
            death\  < 
              Probability[x[n - 1], y[n - 1], w[n - 1]]\ , {y[n]\  = \ 
                y[n - 1], x[n]\  = \ x[n - 1] - 1}, {y[n]\  = \ y[n - 1] - 1, 
              x[n]\  = \ x[n - 1]}]; 
data\  = \ 
            Append[data, {x[n], y[n]}]; 
weather = 
            Append[weather, {t[n], 
                w[n]}]
], 

p1\  = \ 
          ListPlot[data, PlotRange -> {{0, 3000}, {0, 3000}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> False, 
            FrameLabel -> {"Red", "Blue"}, 
            DisplayFunction -> Identity], 
p2\  = \ 
          ListPlot[weather, PlotRange -> {{0, t[nfinal]}, {0, 1}}, 
            PlotJoined -> True, Frame -> True, 
            FrameTicks -> False, 
            FrameLabel -> {"Time", "Weather"}, 
            DisplayFunction -> Identity], 
Show[
          GraphicsArray[{p1, 
              p2}]], 

Print["Final Population: Red ", Part[Last[data], 1], ", ", " Blue ", 
          Part[Last[data], 2], ".", " Duration ", 
          Part[Last[weather], 
            1], "."], 
}, 
{i, 1, 
        20}]\)