Instantly share code, notes, and snippets.

# amitjamadagni/qe_propstepsolvers.jl Last active Aug 29, 2015

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