Stiff methods throw `Cannot create a dual over scalar type ComplexF64` in `SDEProblem`

Specific Domains Modelling & Simulations
1

I am solving the following von Neumann master equation with a stochastic term

\dot{\rho} = -i[H(t), \rho]

where [., .] is the commutator operator and H(t) is given by

H(t) = s(t)\hat{\sigma_x} + \hat{\sigma_z} W_t

s(t) is a control signal and is \sim10^8.
Here is the snippet of the code to solve this using SDEProblem

noise_factor = 1e3 # a prefactor that scales the Wiener process
sigma_x = Matrix{ComplexF64}([[0., 1.] [1., 0.]]);
sigma_z = Matrix{ComplexF64}([[1., 0.] [0., -1.]]);
ham_z = noise_factor * sigma_z
    
# von Neumann master equation
function von_neumann_eq!(dρ, ρ, p, t)
    ham(t) = p(t) * sigma_x
    commutator = ham(t) * ρ .- ρ * ham(t) 
    dρ .=  -im * commutator
end

function stochastic_von_neumann!(dρ, ρ, p, t)
    commutator = ham_z * ρ .- ρ * ham_z
    dρ .= -im * commutator
end

ρ0 = Matrix{ComplexF64}([[1., 0] [0., 0.]])
tspan = (0.0, T) 

prob = SDEProblem(von_neumann_eq!, stochastic_von_neumann!, ρ0, tspan, st)
sol = solve(prob, ImplicitRKMil(), save_everystep=false);

When I specifically choose one of the stiff methods to solve the problem I get the following error

ERROR: LoadError: ArgumentError: Cannot create a dual over scalar type ComplexF64. If the type behaves
 as a scalar, define ForwardDiff.can_dual(::Type{ComplexF64}) = true.

However, when I use alg_hints=[:stiff] without explicitly providing the method, I can see that the algorithm in sol.alg is ImplicitRKMil.

Any suggestions on how to choose explicitly a stiff method in this case?

2

The suggestion I’ve seen is to re-write your equation as two coupled equations in the real and imaginary parts of your solution, so that you are giving the code real variables and autodiff works better.

1 Like
3

You should be able to use ImplicitRKMil(autodiff=false) to disable autodiff and use finite differencing instead.

3 Likes
4

Thanks! I will consider that.

5

This seems nice for the sake of just simulating the problem. My ultimate goal is to optimise the signal s(t). But it was nice to know that there is an option to disable differentiability for a stiff solver. Could you point me to the documentation?