Quantum Equation Addition : Step Solvers

qe_propstepsolvers.jl

##################################################
#####                                        #####
#####            propstepsolvers.jl          #####
#####                                        #####
##################################################


immutable QuEuler <: QuPropagatorMethod
end

immutable QuCrankNicolson <: QuPropagatorMethod
end

immutable QuKrylov <: QuPropagatorMethod
    options::Dict
end

QuKrylov() = QuKrylov(@compat Dict(:NC=>10))

function propagate(prob::QuEuler, eq::QuEquation, t, current_t, current_qustate)
    dt = t - current_t
    return (eye(operator(eq))-im*operator(eq)*dt)*vec(current_qustate)
end

function propagate(prob::QuCrankNicolson, eq::QuEquation, t, current_t, current_qustate)
    dt = t - current_t
    uni = eye(operator(eq))-im*operator(eq)*dt/2
    return \(uni', uni*vec(current_qustate))
end

function propagate(prob::QuKrylov, eq::QuEquation, t, current_t, current_qustate)
    dt = t - current_t
    basis_size = get(prob.options,:NC, length(vec(current_qustate)))
    N = min(basis_size, length(vec(current_qustate)))
    v = Array(typeof(vec(current_qustate)),0)
    @compat sizehint!(v, N+1)
    push!(v,zeros(vec(current_qustate)))
    push!(v,vec(current_qustate))
    alpha = Array(Complex{Float64},0)
    @compat sizehint!(alpha, N)
    beta = Array(Complex{Float64},0)
    @compat sizehint!(beta, N+1)
    push!(beta,0.)
    for i=2:N
        w = operator(eq)*v[i]
        push!(alpha, w'*v[i])
        w = w-alpha[i-1]*v[i]-beta[i-1]*v[i-1]
        push!(beta, norm(w))
        push!(v, w/beta[i])
    end
    w = operator(eq)*v[end]
    push!(alpha, w'*v[end])
    deleteat!(v,1)
    H_k = full(Tridiagonal(beta[2:end], alpha, beta[2:end]))
    U_k = expm(-im*dt*H_k)
    next_state = zeros(vec(current_qustate))
    for j=1:N
        next_state = next_state + v[j]*U_k[j,1]
    end
    return next_state
end