Skip to content

Commit

Permalink
Merge pull request #3238 from AayushSabharwal/as/sde-observed-noise
Browse files Browse the repository at this point in the history
fix: fix SDEs with noise dependent on observed variables
ChrisRackauckas authored Nov 25, 2024
2 parents 6ee68bb + 585de35 commit 14a7239
Showing 3 changed files with 29 additions and 0 deletions.
8 changes: 8 additions & 0 deletions src/structural_transformation/symbolics_tearing.jl
Original file line number Diff line number Diff line change
@@ -88,6 +88,14 @@ function tearing_sub(expr, dict, s)
s ? simplify(expr) : expr
end

function tearing_substitute_expr(sys::AbstractSystem, expr; simplify = false)
empty_substitutions(sys) && return expr
substitutions = get_substitutions(sys)
@unpack subs = substitutions
solved = Dict(eq.lhs => eq.rhs for eq in subs)
return tearing_sub(expr, solved, simplify)
end

"""
$(TYPEDSIGNATURES)
2 changes: 2 additions & 0 deletions src/systems/systems.jl
Original file line number Diff line number Diff line change
@@ -152,6 +152,8 @@ function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = fal
noise_eqs = sorted_g_rows
is_scalar_noise = false
end

noise_eqs = StructuralTransformations.tearing_substitute_expr(ode_sys, noise_eqs)
return SDESystem(full_equations(ode_sys), noise_eqs,
get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys);
name = nameof(ode_sys), is_scalar_noise)
19 changes: 19 additions & 0 deletions test/sdesystem.jl
Original file line number Diff line number Diff line change
@@ -780,3 +780,22 @@ end
prob = @test_nowarn SDEProblem(sys, nothing, (0.0, 1.0))
@test_nowarn solve(prob, ImplicitEM())
end

@testset "Issue#3212: Noise dependent on observed" begin
sts = @variables begin
x(t) = 1.0
input(t)
[input = true]
end
ps = @parameters a = 2
@brownian η

eqs = [D(x) ~ -a * x + (input + 1) * η
input ~ 0.0]

sys = System(eqs, t, sts, ps; name = :name)
sys = structural_simplify(sys)
@test ModelingToolkit.get_noiseeqs(sys) [1.0]
prob = SDEProblem(sys, [], (0.0, 1.0), [])
@test_nowarn solve(prob, RKMil())
end

0 comments on commit 14a7239

Please sign in to comment.