gusucode.com > qit_matlab_0.10.0工具箱源码程序 > qit/@state/propagate.m

    function out = propagate(s, H, t, varargin)
% PROPAGATE  Propagate the state continuously in time.
%
%  out = propagate(s, H, t [, out_func, odeopts])
%  out = propagate(s, L, t [, out_func, odeopts])
%  out = propagate(s, {H, {A_i}}, t [, out_func, odeopts])
%
%  Propagates the state s using the given generator(s) for the time t,
%  returns the resulting state.
%
%  The generator can either be a Hamiltonian matrix H or, for time-dependent
%  Hamiltonians, a function handle H(t) which takes a time instance t
%  as input and return the corresponding H matrix.
%
%  Alternatively, the generator can also be a Liouvillian superoperator, or
%  a list consisting of a Hamiltonian and a list of Lindblad operators.
%
%  If t is a vector of increasing time instances, returns a cell array
%  containing the propagated state for each time given in t.
%
%  Optional parameters (can be given in any order):
%    out_func: If given, for each time instance propagate returns out_func(s(t), H(t)).
%    odeopts:  Options struct for MATLAB ODE solvers from the odeset function.

%  out == expm(-i*H*t)*|s>
%  out == inv_vec(expm(L*t)*vec(\rho_s))

% Ville Bergholm 2008-2010
% James Whitfield 2009


if (nargin < 3)
  error('Needs a state, a generator and a time.');
end

out_func = @(x,h) x; % if no out_func is given, use a NOP

odeopts = odeset('RelTol', 1e-4, 'AbsTol', 1e-6, 'Vectorized', 'on');

n = length(t); % number of time instances we are interested in
out = cell(1, n);
dim = size(s.data); % system dimension

if (isa(H, 'function_handle'))
  % time dependent
  t_dependent = true;
  F = H;
  H = F(0);
else
  % time independent
  t_dependent = false;
end

if (isnumeric(H))
  % matrix
  dim_H = size(H, 2);

  if (dim_H == dim(1))
    gen = 'H';
  elseif (dim_H == dim(1)^2)
    gen = 'L';
    s = to_op(s);
  else
    error('Dimension of the generator does not match the dimension of the state.');
  end
  
elseif (iscell(H))
  % list of Lindblad operators
  dim_H = size(H{1}, 2);
  if (dim_H == dim(1))
    gen = 'A';
    s = to_op(s);

    % HACK, in this case we use ode45 anyway
    if (~t_dependent)
      t_dependent = true; 
      F = @(t) H; % ops stay constant
    end
  else
    error('Dimension of the Lindblad ops does not match the dimension of the state.');
  end

else
  error(['The second parameter has to be either a matrix, a cell array, '...
         'or a function handle that returns a matrix or a cell array.'])
end

dim = size(s.data); % may have been switched to operator representation


% process optional arguments
for k=1:nargin-3
  switch class(varargin{k})
    case 'function_handle'
      out_func = varargin{k};

    case 'struct'
      odeopts = odeset(odeopts, varargin{k});

    otherwise
      error('Unknown optional parameter type.');
  end
end


% derivative functions for the solver

function d = lindblad_fun(t, y, F)
  X = F(t);
  A = X{2};
  A = A(:);

  d = zeros(size(y));
  % lame vectorization
  for loc1_k=1:size(y, 2)
    x = reshape(y(:,loc1_k), dim); % into a matrix
  
    % Hamiltonian
    temp = -1i * (X{1} * x  -x * X{1});
    % Lindblad operators
    for j=1:length(A)
      ac = A{j}'*A{j};
      temp = temp +A{j}*x*A{j}' -0.5*(ac*x + x*ac);
    end
    d(:,loc1_k) = temp(:); % back into a vector
  end
end

function d = mixed_fun(t, y, F)
  H = F(t);
  
  d = zeros(size(y));
  % vectorization
  for loc2_k=1:size(y, 2)
    x = reshape(y(:,loc2_k), dim); % into a matrix
    temp = -1i * (H * x  -x * H);
    d(:,loc2_k) = temp(:); % back into a vector
  end
end


if (t_dependent)
  % time dependent case, use ODE solver

  switch (gen)
    case 'H'
      % Hamiltonian
      if (dim(2) == 1)
        % pure state
        odefun = @(t, y, F) -1i * F(t) * y; % derivative function for the solver
      else
        % mixed state
        odefun = @mixed_fun;
      end

    case 'L'
      % Liouvillian
      odefun = @(t, y, F) F(t) * y;
      %odeopts = odeset(odeopts, 'Jacobian', F);

    case 'A'
      % Hamiltonian and Lindblad operators in a list
      odefun = @lindblad_fun;
  end

  skip = 0;
  if (t(1) ~= 0)
    t = [0, t]; % ODE solver needs to be told that t0 = 0
    skip = 1;
  end

  if (length(t) < 3)
    t = [t, t(end)+1e5*eps]; % add another time point to get reasonable output from solver
  end

  %odeopts = odeset(odeopts, 'OutputFcn', @(t,y,flag) odeout(t, y, flag, H));

  [t_out, s_out] = ode45(odefun, t, s.data, odeopts, F);
  % s_out has states in columns, row i corresponds to t(i)

  % apply out_func
  for k=1:n
    % this works because ode45 automatically expands input data into a col vector
    s.data = inv_vec(s_out(k+skip,:), dim);
    out{k} = out_func(s, F(t_out(k+skip)));
  end

else
  % time independent case

  switch (gen)
    case 'H'
      if (length(H) < 500)
        % eigendecomposition
        [v, d] = eig(full(H)); % TODO eigs?
        d = diag(d);
        for k=1:n
          U = v * diag(exp(-1i * t(k) * d)) * v';
          out{k} = out_func(u_propagate(s, U), H);
          %out{k} = out_func(u_propagate(s, expm(-i*H*t(k))), H);
        end
      else
        % Krylov subspace method
        [w, err] = expv(-1i*t, H, s.data);
        for k=1:n
          s.data = w(:,k);
          out{k} = out_func(s, H);
        end
      end
      
    case 'L'
      % Krylov subspace method
      [w, err] = expv(t, H, vec(s.data));
      for k=1:n
        s.data = inv_vec(w(:,k));
        out{k} = out_func(s, H);
      end
  end
end

if (n == 1)
  % single output, don't bother with a list
  out = out{1};
end
end


%function status = odeout(t, y, flag, H)
%if isempty(flag)
%  sdfsd
%end
%status = 0;
%end