function x = S(A,i,j,t,index)
% S = S(A,i,j,t,index)
%
% This function calculates the quantity S(t), 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.
% S(t) can be defined as the contribution to the near-term sensitivity which requires no
% additional (compared to an asymptotic sensitivity analysis) data collection.
% S(t) describes the near-term sensitivity when a population is initially at a stable
% stage distribution, and it can be calculated from the Leslie matrix alone.
%
% 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 an initial perturbation for which the sensitivity is
% being calculated
%
% Optional argument:
% index defines whether a component of S is being computed. S can be
% decomposed into components from each eigenvector.
% If index=0 then calculate S(t) in full
% If index=i (i>0) then calculate ith term corresponding to the
% contribution from the ith eigenvector
%
% Output:
% The value for S(t), which is a vector with the same dimensions as the
% population structure vector.
% 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
% Largest eigenvalue is lambda(N), smallest is lambda(1)
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 S
x = 0;
for k=list,
if k==N,
x = x + lambda(N)^t * (conj(v(i,k))*w(j,k) * t*w(:,N)/lambda(N) + ...
sum(w(:,1:N-1).*(ones(N,1)*(conj(v(i,1:N-1))*w(j,N)./(lambda(N)-lambda(1:N-1).'))),2));
else
x = x - lambda(k)^t * conj(v(i,k))*w(j,N)' * w(:,k)/(lambda(N)-lambda(k));
end
end