Relative entropy as applied to an evolutionary competition between two species

evolution.m

function [T Y] = relent(Tspan)

    %% Configuration

    leg = {'Species 1','Species 2'};

    opts = odeset();
    opts.AbsTol = 1e-8;
    opts.RelTol = 1e-4;

    %% Parameters and initial condition

    a = 1;
    b = 1;
    c = 1;
    d = 1;
    e = 0.1;

    %% Define the evolutionarily stable state

    C0 = [c/d + e*a/(b*d), a/b]
    c0 = C0 / sum(C0)

    %% Initial conditions

    Y0 = [1 2];

    %% Solve the equation

    [T,Y] = ode45(@f,[0 Tspan],Y0,opts);

    %% Make plots

    subplot(2,2,1)
    plot(T,Y)
    ylabel('Absolute number')
    legend(leg)

    hline(C0)

    P = sum(Y,2);
    y = bsxfun(@rdivide,Y,P);

    subplot(2,2,2)
    plot(T,y)
    ylabel('Proportion')
    legend(leg)

    hline(c0)

    S = entropy(y);

    subplot(2,2,3)
    plot(T,S)
    ylabel('Entropy')

    c0 = repmat(c0,[size(Y,1) 1]);
    I = relative_entropy(c0,y);

    subplot(2,2,4)
    plot(T,I)
    ylabel('Relative entropy')

    %% Function definitions

    function f = f(t,y)
        f = zeros(2,1);
        f(1) = a * y(1) - b * y(1) * y(2);
        f(2) = -(c * y(2) - d * y(1) * y(2) + e* y(2)^2);
    end

    function S = entropy(p)
        S = -sum(p .* log(p),2);
    end

    function I = relative_entropy(q,p)
        I = sum(q .* log(q./p),2);
    end

end