function delta_S = delta_S(A,i,j,t,index)
% delta_S = delta_S(A,i,j,t,index)
%
% This function calculates the quantity Delta S, defined in the section
% "decomposing near-term sensitivity analysis" of the paper by J. Yearsley,
% "Transient population dynamics and short-term sensitivity analysis of matrix
% population models", Ecological Modelling.
% Delta S can be defined as the contribution to the near-term sensitivity
% which requires additional (compared to an asymptotic sensitivity analysis)
% data collection.
% Delta S describes the contribution to the near-term sensitivity due to
% deviations from a stable stage distribution, and it can be calculated from the Leslie matrix.
%
% Input arguments:
% A is the Leslie matrix
% i,j are the indices describing the element of the Leslie matrix for
% which a sensitivity is being calculated.
% t is the the time after the perturbation for which the sensitivity is
% being calculated
%
% Optional argument:
% index defines which a component of Delta S is being computed. Delta S can be
% decomposed into components from each eigenvector.
% If index=0 then calculate Delta S in full
% If index=i (i>0) then calculate ith term corresponding to the
% contribution from the ith eigenvector
%
% Output:
% delta_S, which is a matrix with the same dimensions as the Leslie
% matrix
% Jon Yearsley 22/7/2003
% j.yearsley@macaulay.ac.uk
if nargin<5,
index=0;
end
% Right eigenvectors and eigenvalues of Leslie matrix A
[w, lambda] = eig(A);
[~, ind]=sort(real(diag(lambda)));
lambda = diag(lambda);
lambda = lambda(ind);
% Sort the right eigenvectors into order
w = w(:,ind);
% Calculate left eigenvectors
v = inv(w)';
N = length(lambda);
if index==0 || index>N,
% Add up all components
list = 1:N;
else
% Look at only one component
list = index;
end
% Calculate delta S
delta_S = 0;
for k=list,
list2 = 1:N;
list2(k) = [];
delta_S = delta_S + t * lambda(k)^(t-1) * conj(v(i,k))*w(j,k).' * w(:,k)*v(:,k)';
for m=list2,
delta_S = delta_S + lambda(k)^t * ...
(conj(v(i,k))*w(j,m)* w(:,k)*v(:,m)' + conj(v(i,m))*w(j,k) * w(:,m)*v(:,k)')./(lambda(k)-lambda(m));
end
end