Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Quantum Equation Addition : Step Solvers
##################################################
##### #####
##### 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment