Math 448/548 - Example 6.1 - War Problem

This Matlab code illustrates the dependence on paramenter analysis for a discrete dymanical system via repeated simulations.

clear all, close all, clc

endtime = []; winner = [];
values = [1 1.5 2 3 5 6]; % choose a set of values for lambda

Plot the direction field (no rescaling!)

for lambda = values
    a2=0.05; b2=0.005;
    a1=lambda*a2; b1 =lambda*b2;

    figure
    M=10; % number of samples points
    x1min=0; x1max=5; % domain specification
    x2min=0; x2max=2;
    [X1,X2]=meshgrid(x1min:(x1max-x1min)/M:x1max,x2min:(x2max-x2min)/M:x2max);
    dX1=-a1*X2-b1*X1.*X2; % x1-component
    dX2=-a2*X1-b2*X1.*X2;         % x2-component
    quiver(X1,X2,dX1,dX2,0); % matlab routine
    axis([x1min x1max x2min x2max]);
    hold on
    xlabel('Red Army'); ylabel('Blue Army');

    x=[5;2];
    n=1;
    sprintf(strcat('Battle evolution for lambda =',num2str(lambda),'\n'))
    fprintf('n    Red Army    Blue Army  \n')
    fprintf('%2.0f      %5.2f        %5.2f\n', 0, x(1), x(2))

    while x(1)>0 & x(2)>0
        n=n+1;
        rhs = [ -a1*x(2)-b1*x(1)*x(2);
                -a2*x(1)-b2*x(1)*x(2)];
        xnew = x + rhs;
        plot([x(1),xnew(1)], [x(2),xnew(2)],'--ro',...
        'MarkerFaceColor','k','MarkerSize',4)
        x = xnew;
        fprintf('%2.0f      %5.2f        %5.2f\n', n, x(1), x(2))
    end
    if x(1)>0 winner = [winner,1];
        else winner = [winner,2];
    end
    title(strcat('Battle outcome for \lambda =',num2str(lambda)))
    xlabel('Red Army'); ylabel('Blue Army');
    endtime = [endtime n];
end

%endtime=endtime(2:end);
ans =

Battle evolution for lambda =1


n    Red Army    Blue Army  
 0       5.00         2.00
 2       4.85         1.70
 3       4.72         1.42
 4       4.62         1.15
 5       4.54         0.89
 6       4.47         0.64
 7       4.42         0.40
 8       4.40         0.17
 9       4.38        -0.05

ans =

Battle evolution for lambda =1.5


n    Red Army    Blue Army  
 0       5.00         2.00
 2       4.78         1.70
 3       4.59         1.42
 4       4.43         1.16
 5       4.31         0.91
 6       4.21         0.68
 7       4.14         0.45
 8       4.09         0.24
 9       4.06         0.03
10       4.06        -0.18

ans =

Battle evolution for lambda =2


n    Red Army    Blue Army  
 0       5.00         2.00
 2       4.70         1.70
 3       4.45         1.43
 4       4.24         1.17
 5       4.08         0.93
 6       3.95         0.71
 7       3.85         0.50
 8       3.78         0.30
 9       3.74         0.10
10       3.72        -0.09

ans =

Battle evolution for lambda =3


n    Red Army    Blue Army  
 0       5.00         2.00
 2       4.55         1.70
 3       4.18         1.43
 4       3.87         1.19
 5       3.63         0.98
 6       3.43         0.78
 7       3.27         0.59
 8       3.15         0.42
 9       3.07         0.26
10       3.02         0.10
11       3.00        -0.05

ans =

Battle evolution for lambda =5


n    Red Army    Blue Army  
 0       5.00         2.00
 2       4.25         1.70
 3       3.64         1.45
 4       3.15         1.24
 5       2.74         1.07
 6       2.40         0.91
 7       2.12         0.78
 8       1.88         0.67
 9       1.68         0.57
10       1.52         0.48
11       1.38         0.40
12       1.26         0.33
13       1.17         0.26
14       1.10         0.20
15       1.04         0.15
16       1.00         0.09
17       0.98         0.04
18       0.96        -0.01

ans =

Battle evolution for lambda =6


n    Red Army    Blue Army  
 0       5.00         2.00
 2       4.10         1.70
 3       3.38         1.46
 4       2.79         1.27
 5       2.31         1.11
 6       1.90         0.98
 7       1.55         0.88
 8       1.25         0.79
 9       0.98         0.73
10       0.74         0.67
11       0.52         0.63
12       0.32         0.61
13       0.14         0.59
14      -0.04         0.58
disp('lambda    battle length    winner')
disp([values' endtime' winner'])
lambda    battle length    winner
    1.0000    9.0000    1.0000
    1.5000   10.0000    1.0000
    2.0000   10.0000    1.0000
    3.0000   11.0000    1.0000
    5.0000   18.0000    1.0000
    6.0000   14.0000    2.0000