% a demonstration that considers the system
%    x' = [ -r1/4 -r1/4; r2 -r2 ] x
% which is a model for the populations of baleen whales (x2)
% and krill (x1), near their scaled equlibrium populations.
% the nonlinear system is
%    x1' = r1 x1 ( .5 - x1 - x2 )
%    x2' = r2 x2 ( 1 - x2/x1 )
%
% expected values of r1 and r2 are r1 > 0.3 and r2 = 0.4.
%
% solutions to the linear and nonlinear problem are shown, along with
% a phase portrait for the system.
%
% by Gavin LaRose, 5/16/2013
% 
% (c)2013 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.
%

function Whales_Krill_Coexist()
%-------------------------------------------------------------------------------
% configuration variables

    r1 = 0.4;
    r2 = 0.4;
    
% initial condition for solution curves
    x10 = 0.3;
    x20 = 0.1;
    
% should we show the direction field?
    show_direction_field = 1;
% number of vectors to graph in each direction
    num_vec = 15;
    
% should we put a pause between graphing the direction field and linear
% solutions and trajectories, and graphing the nonlinear trajectories?
    include_pause = 1;
    
% the distance from the equilibrium solutions to plot the trajectories
% for the linearized systems
    linear_radius = 0.15;
    
% colors
    ltraj_color = 'b';  % the color for linear trajectories    
    ntraj_color = 'm';  % and for the nonlinear trajectories

    df_color    = 'k';  % the color of the direction field, if any
    
    background_grey = [ 0.6 0.6 0.6 ];
    
% size of the phase plane and solution graph axes
    pp_size = [ 0 0.8 0 0.8 ];
    sg_size = [ 0 20 0 0.8 ];

%-------------------------------------------------------------------------------
% initialization

% the equilibrium point
    equil = [ 1/4; 1/4 ];
    
% the coefficient matrix of the linearized system at this point
    a_mat = [ -r1/4 -r1/4; r2 -r2 ];
    
% linear and nonlinear solutions 
    [ t_li, x_li ] = ode45( @(t,x) [ a_mat*x ], [ 0 20 ], ...
                            [ x10-equil(1); x20-equil(2) ] );
    [ t_nl, x_nl ] = ode45( @(t,x) [ r1*x(1)*(0.5 - x(1) - x(2)); ...
                                     r2*x(2)*(1 - x(2)/x(1)) ], ...
                            [ 0 20 ], [ x10; x20 ] );
    
% set up the figure window
    scrsz    = get(0,'ScreenSize');
    scrsz(3) = min(1500, scrsz(3));
    scrsz(4) = min(900, scrsz(4));
    
    fig1 = figure( 'OuterPosition', ...
                   [50 scrsz(4)/2 scrsz(3)/1.5 scrsz(3)/2.5], ...
                   'Name', 'Whales/Krill' );
    solnaxes = subplot( 1,2,1, 'Parent', fig1, 'Color', background_grey );
    pplnaxes = subplot( 1,2,2, 'Parent', fig1, 'Color', background_grey );

    hold(solnaxes);
    hold(pplnaxes);
    axis( pplnaxes, pp_size );
    axis( solnaxes, sg_size );
    
    xlabel( solnaxes, 't', 'FontSize', 15 );
    ylabel( solnaxes, 'x1,x2', 'FontSize', 15 );
    xlabel( pplnaxes, 'x1', 'FontSize', 15 );
    ylabel( pplnaxes, 'x2', 'FontSize', 15 );

    title( solnaxes, ['Solutions for whale/krill model'], 'FontSize', 15 );
    title( pplnaxes, ['Phase plane for whale/krill model'], 'FontSize', 15 );

%-------------------------------------------------------------------------------
% graph the trajectories and solutions for the linearized system

    if show_direction_field
	% set up x and y vectors for the direction field
        x_vec = linspace( 0, 1, num_vec );
        y_vec = linspace( 0, 1, num_vec );

	% define an x  and a y grid and the corresponding derivatives
        x_grid  = [ ];
        y_grid  = [ ];
        xp_grid = [ ];
        yp_grid = [ ];
    
        for j1=1:num_vec
            for j2=1:num_vec
                x_grid = [ x_grid x_vec(j1) ];
                y_grid = [ y_grid y_vec(j2) ];
            
                der = [ r1*x_vec(j1)*(0.5 - x_vec(j1) - y_vec(j2)); ...
                        r2*y_vec(j2)*(1 - y_vec(j2)/x_vec(j1)) ];

                xp_grid = [ xp_grid der(1) ];
                yp_grid = [ yp_grid der(2) ];
            end
        end
    
        dx_len = sqrt( xp_grid.^2 + yp_grid.^2 );
    
        quiver( pplnaxes, ...
                x_grid, y_grid, xp_grid./dx_len, yp_grid./dx_len, ...
                'Color', df_color, 'ShowArrowHead', 'on', ...
                'AutoScaleFactor', .3 );
    end
    
% plot the linear solutions
    plot( solnaxes, t_li, x_li(:,1)+equil(1), ['-' ltraj_color], ...
          t_li, x_li(:,2)+equil(2), ['--' ltraj_color], 'LineWidth', 2 );
    sleg = legend( solnaxes, 'x1 (linear)', 'x2 (linear)' );

% next plot the local trajectories near the equilibrium point
    plot_localtraj( pplnaxes, equil, a_mat, linear_radius, ...
                    [ '-' ltraj_color ], [ ], {'LineWidth', 1} );
    plot( pplnaxes, x_li(:,1)+equil(1), x_li(:,2)+equil(2), ...
          ['-' ltraj_color], 'LineWidth', 2 );
    
% add the equilibrium solution and initial condition
    plot( solnaxes, [ 0 20 ], [ 0.25 0.25 ], '-k' );
    plot( pplnaxes, [equil(1)], [equil(2)], '.k', 'MarkerSize', 20 );
    plot( pplnaxes, [x10], [x20], '.k', 'MarkerSize', 20 );
    
    if include_pause
        pause on;
        disp( 'Press any key to show nonlinear solutions' );
        pause
    end
    
% and finally plot solution curves for the nonlinear system
    plot( solnaxes, t_nl, x_nl(:,1), ['-' ntraj_color], ...
          t_nl, x_nl(:,2), ['--' ntraj_color], 'LineWidth', 2 );
    delete(sleg)
    legend( solnaxes, 'x1 (linear)', 'x2 (linear)', ...
            'x1 (nonlinear)', 'x2 (nonlinear)' );
    
% and the solution trajectory
    plot( pplnaxes, x_nl(:,1), x_nl(:,2), ['-' ntraj_color], ...
          'LineWidth', 2 );
    
    if include_pause
        disp( 'Press any key to show nonlinear trajectories' );
        pause
    end

    ic_list = [ [ 0.1; 0.1 ], [ 0.5; 0.1 ], [ 0.7; 0.25 ], ...
                [ 0.7; 0.5 ], [ 0.7; 0.7 ], [ 0.5; 0.7 ], ...
                [ 0.3; 0.7 ], [ 0.1; 0.7 ], [ 0.1; 0.5 ] ];
    
    for j1 = 1:length(ic_list)
        [t_sol, x_sol] = ode45( @(t,x) [ r1*x(1)*(0.5 - x(1) - x(2)); ...
                                         r2*x(2)*(1 - x(2)/x(1)) ], ...
                                [ 0 20 ], ic_list(:,j1) );
        plot( pplnaxes, x_sol(:,1), x_sol(:,2), ['-' ntraj_color], ...
              'LineWidth', 2 );
    end
    
end

% all done
%-------------------------------------------------------------------------------
