Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Complex array-valued Wiener process only produces real noise #101

Closed
danielwe opened this issue May 17, 2022 · 4 comments
Closed

Complex array-valued Wiener process only produces real noise #101

danielwe opened this issue May 17, 2022 · 4 comments

Comments

@danielwe
Copy link
Contributor

WienerProcess(0.0, [0.0im]) only produces real-valued increments. MWE:

julia> using DiffEqNoiseProcess

julia> using DiffEqNoiseProcess: NoiseProblem, solve

julia> prob = NoiseProblem(WienerProcess(0.0, [0.0im]), (0.0, 1.0));

julia> sol = solve(prob; dt=0.1);

julia> sol.W[end]
1-element Vector{ComplexF64}:
 2.3820176552198165 + 0.0im

The scalar-valued equivalent produces complex noise as expected:

julia> prob = NoiseProblem(WienerProcess(0.0, 0.0im), (0.0, 1.0));

julia> sol = solve(prob; dt=0.1);

julia> sol.W[end]
-1.1828100925799103 + 1.4213457901396127im

Note that this doesn't seem to affect the default noise process for complex-valued SDEs in StochasticDiffEq. I'm just not able to manually instantiate a noise process that reproduces the default behavior:

julia> using StochasticDiffEq

julia> f!(du, u, p, t) = (du .= 0)
f! (generic function with 1 method)

julia> g!(du, u, p, t) = (du .= 1)
g! (generic function with 1 method)

julia> prob = SDEProblem(f!, g!, [0.0im], (0.0, 1.0));

julia> sol = solve(prob, EM(); dt=0.1);

julia> sol.u[end]
1-element Vector{ComplexF64}:
 -0.9054409033854848 - 0.4148766416920265im
@danielwe
Copy link
Contributor Author

Oh, I see, the in-place WienerProcess! works. I'm confused: what's the difference between the two when initialized with an array?

@danielwe
Copy link
Contributor Author

Looks like WienerProcess without the bang will eventually hit the following method, which fails to make use of the eltype T in the call to randn. That explains it.

@inline function wiener_randn(rng::AbstractRNG,proto::Array{T}) where T
randn(rng,size(proto))
end

@ChrisRackauckas
Copy link
Member

Thanks for the investigation. Yeah this shows we should probably upgrade it to the new Base stuff.

@ChrisRackauckas
Copy link
Member

Fixed by #157

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants