From add7059c8a52df7401aee9af3f7e6d8833b59db4 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 28 Oct 2024 17:14:30 -0400 Subject: [PATCH 01/18] make Observable definition comply with new record_from_solution definition --- ext/MTKBifurcationKitExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl index a0bd3bfc02..e5b57efc5a 100644 --- a/ext/MTKBifurcationKitExt.jl +++ b/ext/MTKBifurcationKitExt.jl @@ -59,7 +59,7 @@ struct ObservableRecordFromSolution{S, T} end end # Functor function that computes the value. -function (orfs::ObservableRecordFromSolution)(x, p) +function (orfs::ObservableRecordFromSolution)(x, p; k...) # Updates the state values (in subs_vals). for state_idx in 1:(orfs.state_end_idxs) orfs.subs_vals[state_idx] = orfs.subs_vals[state_idx][1] => x[state_idx] From eb681c10403031e3a1ee28c96f9d4b2defc3e389 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 28 Oct 2024 17:26:25 -0400 Subject: [PATCH 02/18] add test --- test/extensions/bifurcationkit.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl index b8c95dca99..629edf46a6 100644 --- a/test/extensions/bifurcationkit.jl +++ b/test/extensions/bifurcationkit.jl @@ -162,4 +162,12 @@ let bf = bifurcationdiagram(bp, PALC(), 2, opts_br) @test bf.γ.specialpoint[1].param≈0.1 atol=1e-4 rtol=1e-4 + + # Test with plot variable as observable + pvar = ModelingToolkit.get_var_to_name(fol)[:RHS] + bp = BifurcationProblem(fol, u0, par, bif_par; plot_var = pvar) + opts_br = ContinuationPar(p_min = -1.0, + p_max = 1.0) + bf = bifurcationdiagram(bp, PALC(), 2, opts_br) + @test bf.γ.specialpoint[1].param≈0.1 atol=1e-4 rtol=1e-4 end From 9733460668a662b41fd4cf87102d08521c1e4f92 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 13:25:30 -0500 Subject: [PATCH 03/18] init --- src/systems/diffeqs/abstractodesystem.jl | 110 +++++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 8d9e0b5381..2579ebafcf 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -466,6 +466,116 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end +""" +```julia +SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` + +Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and +`ps` are used to set the order of the dependent variable and parameter vectors, +respectively. +""" +function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem(sys::AbstractODESystem, + u0map::StaticArray, + args...; + kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) +end + +function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], + tspan = get_tspan(sys), + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + sparsity = true, + callback = nothing, + check_length = true, + warn_initialize_determined = true, + eval_expression = false, + eval_module = @__MODULE__, + kwargs...) where {iip, specialize} + + if !iscomplete(sys) + error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") + end + + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, + check_length, warn_initialize_determined, eval_expression, eval_module, jac, sparse, sparsity, kwargs...) + + # if jac + # jac_gen = generate_jacobian(sys, dvs, ps; + # simplify = simplify, sparse = sparse, + # expression = Val{true}, + # expression_module = eval_module, + # checkbounds = checkbounds, kwargs...) + # jac_oop, jac_iip = eval_or_rgf.(jac_gen; eval_expression, eval_module) + + # _jac(u, p, t) = jac_oop(u, p, t) + # _jac(J, u, p, t) = jac_iip(J, u, p, t) + # _jac(u, p::Tuple{Vararg{Number}}, t) = jac_oop(u, p, t) + # _jac(J, u, p::Tuple{Vararg{Number}}, t) = jac_iip(J, u, p, t) + # _jac(u, p::Tuple, t) = jac_oop(u, p..., t) + # _jac(J, u, p::Tuple, t) = jac_iip(J, u, p..., t) + # _jac(u, p::MTKParameters, t) = jac_oop(u, p..., t) + # _jac(J, u, p::MTKParameters, t) = jac_iip(J, u, p..., t) + # else + # _jac = nothing + # end + + # jac_prototype = if sparse + # uElType = u0 === nothing ? Float64 : eltype(u0) + # if jac + # similar(calculate_jacobian(sys, sparse = sparse), uElType) + # else + # similar(jacobian_sparsity(sys), uElType) + # end + # else + # nothing + # end + + # f.jac = _jac + # f.jac_prototype = jac_prototype + # f.sparsity = jac ? jacobian_sparsity(sys) : nothing + + cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) + + kwargs = filter_kwargs(kwargs) + + kwargs1 = (;) + if cbs !== nothing + kwargs1 = merge(kwargs1, (callback = cbs,)) + end + + # Define the boundary conditions + bc = if iip + (residual, u, p, t) -> (residual = u[1] - u0) + else + (u, p, t) -> (u[1] - u0) + end + + return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) +end + +get_callback(prob::BVProblem) = prob.kwargs[:callback] + """ ```julia DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys), From b3da8137a9288b9c3eeff209a541dd41a5f97524 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 16:18:18 -0500 Subject: [PATCH 04/18] up --- src/systems/diffeqs/abstractodesystem.jl | 64 +++++++++--------------- 1 file changed, 23 insertions(+), 41 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 7a8bf56562..bb4cbbb6c3 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -469,6 +469,17 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end +""" +```julia +SciMLBase.BVPFunction{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` +""" + """ ```julia SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, @@ -481,7 +492,7 @@ SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and `ps` are used to set the order of the dependent variable and parameter vectors, -respectively. +respectively. `u0` should be either the initial condition, a vector of values `u(t_i)` for collocation methods, or a function returning one or the other. """ function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) BVProblem{true}(sys, args...; kwargs...) @@ -502,12 +513,13 @@ function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end +# figure out what's going on when we try to set `sparse`? + function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); version = nothing, tgrad = false, jac = true, sparse = true, - sparsity = true, callback = nothing, check_length = true, warn_initialize_determined = true, @@ -521,57 +533,27 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, jac, sparse, sparsity, kwargs...) - - # if jac - # jac_gen = generate_jacobian(sys, dvs, ps; - # simplify = simplify, sparse = sparse, - # expression = Val{true}, - # expression_module = eval_module, - # checkbounds = checkbounds, kwargs...) - # jac_oop, jac_iip = eval_or_rgf.(jac_gen; eval_expression, eval_module) - - # _jac(u, p, t) = jac_oop(u, p, t) - # _jac(J, u, p, t) = jac_iip(J, u, p, t) - # _jac(u, p::Tuple{Vararg{Number}}, t) = jac_oop(u, p, t) - # _jac(J, u, p::Tuple{Vararg{Number}}, t) = jac_iip(J, u, p, t) - # _jac(u, p::Tuple, t) = jac_oop(u, p..., t) - # _jac(J, u, p::Tuple, t) = jac_iip(J, u, p..., t) - # _jac(u, p::MTKParameters, t) = jac_oop(u, p..., t) - # _jac(J, u, p::MTKParameters, t) = jac_iip(J, u, p..., t) - # else - # _jac = nothing - # end - - # jac_prototype = if sparse - # uElType = u0 === nothing ? Float64 : eltype(u0) - # if jac - # similar(calculate_jacobian(sys, sparse = sparse), uElType) - # else - # similar(jacobian_sparsity(sys), uElType) - # end - # else - # nothing - # end - - # f.jac = _jac - # f.jac_prototype = jac_prototype - # f.sparsity = jac ? jacobian_sparsity(sys) : nothing + check_length, warn_initialize_determined, eval_expression, eval_module, jac, kwargs...) cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) - kwargs = filter_kwargs(kwargs) kwargs1 = (;) if cbs !== nothing kwargs1 = merge(kwargs1, (callback = cbs,)) end + + # Construct initial conditions + _u0 = prepare_initial_state(u0) + __u0 = if _u0 isa Function + _u0(t_i) + end # Define the boundary conditions bc = if iip - (residual, u, p, t) -> (residual = u[1] - u0) + (residual, u, p, t) -> (residual = u[1] - __u0) else - (u, p, t) -> (u[1] - u0) + (u, p, t) -> (u[1] - __u0) end return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) From 4affeac4b340a06c770373b498ec6b0e94c25a06 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 17:35:05 -0500 Subject: [PATCH 05/18] up --- src/systems/diffeqs/abstractodesystem.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 263a75d153..33bddece30 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -545,9 +545,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] # Construct initial conditions _u0 = prepare_initial_state(u0) - __u0 = if _u0 isa Function - _u0(t_i) - end + __u0 = _u0 isa Function ? _u0(tspan[1]) : _u0 # Define the boundary conditions bc = if iip From a3429ea2b7d9e67898023eac1599e71c3d1b7bec Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 17:42:39 -0500 Subject: [PATCH 06/18] up --- src/systems/diffeqs/abstractodesystem.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 33bddece30..84f40a01f7 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -513,8 +513,6 @@ function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end -# figure out what's going on when we try to set `sparse`? - function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); @@ -544,14 +542,13 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] end # Construct initial conditions - _u0 = prepare_initial_state(u0) - __u0 = _u0 isa Function ? _u0(tspan[1]) : _u0 + _u0 = u0 isa Function ? u0(tspan[1]) : u0 # Define the boundary conditions bc = if iip - (residual, u, p, t) -> (residual = u[1] - __u0) + (residual, u, p, t) -> (residual = u[1] - _u0) else - (u, p, t) -> (u[1] - __u0) + (u, p, t) -> (u[1] - _u0) end return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) From f751fbb51c49e0b35ee84a12c54ec2de919660a5 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 21:22:08 -0500 Subject: [PATCH 07/18] up --- src/systems/diffeqs/abstractodesystem.jl | 18 ++------ test/bvproblem.jl | 56 ++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 15 deletions(-) create mode 100644 test/bvproblem.jl diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 84f40a01f7..112419dea8 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -469,17 +469,6 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end -""" -```julia -SciMLBase.BVPFunction{iip}(sys::AbstractODESystem, u0map, tspan, - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - jac = true, sparse = true, - simplify = false, - kwargs...) where {iip} -``` -""" - """ ```julia SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, @@ -492,7 +481,7 @@ SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and `ps` are used to set the order of the dependent variable and parameter vectors, -respectively. `u0` should be either the initial condition, a vector of values `u(t_i)` for collocation methods, or a function returning one or the other. +respectively. `u0map` should be used to specify the initial condition, or be a function returning an initial condition. """ function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) BVProblem{true}(sys, args...; kwargs...) @@ -517,7 +506,6 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); version = nothing, tgrad = false, - jac = true, sparse = true, callback = nothing, check_length = true, warn_initialize_determined = true, @@ -531,7 +519,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, jac, kwargs...) + check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) kwargs = filter_kwargs(kwargs) @@ -546,7 +534,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] # Define the boundary conditions bc = if iip - (residual, u, p, t) -> (residual = u[1] - _u0) + (residual, u, p, t) -> residual .= u[1] - _u0 else (u, p, t) -> (u[1] - _u0) end diff --git a/test/bvproblem.jl b/test/bvproblem.jl new file mode 100644 index 0000000000..90d41a96ad --- /dev/null +++ b/test/bvproblem.jl @@ -0,0 +1,56 @@ +using BoundaryValueDiffEq, OrdinaryDiffEq +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D + +@parameters σ = 10. ρ = 28 β = 8/3 +@variables x(t) = 1 y(t) = 0 z(t) = 0 + +eqs = [D(x) ~ σ*(y-x), + D(y) ~ x*(ρ-z)-y, + D(z) ~ x*y - β*z] + +u0map = [:x => 1., :y => 0., :z => 0.] +parammap = [:ρ => 28., :β => 8/3, :σ => 10.] +tspan = (0., 10.) + +@mtkbuild lorenz = ODESystem(eqs, t) + +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) +sol = solve(bvp, MIRK4(), dt = 0.1); + +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) +sol2 = solve(bvp, MIRK4(), dt = 0.1); + +op = ODEProblem(lorenz, u0map, tspan, parammap) +osol = solve(op) + +@test sol.u[end] ≈ osol.u[end] +@test sol2.u[end] ≈ osol.u[end] +@test sol.u[1] == [1., 0., 0.] +@test sol2.u[1] == [1., 0., 0.] + +### Testing on pendulum + +@parameters g = 9.81 L = 1. +@variables θ(t) = π/2 + +eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] + +@mtkbuild pend = ODESystem(eqs, t) + +u0map = [θ => π/2, D(θ) => π/2] +parammap = [:L => 2.] +tspan = (0., 10.) + +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) +sol = solve(bvp, MIRK4(), dt = 0.05); + +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) +sol2 = solve(bvp2, MIRK4(), dt = 0.05); + +osol = solve(pend) + +@test sol.u[end] ≈ osol.u[end] +@test sol.u[1] == [π/2, π/2] +@test sol2.u[end] ≈ osol.u[end] +@test sol2.u[1] == [π/2, π/2] From a9fdfd6a3a114c4df28e4781c6b667990c01bafc Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 3 Dec 2024 05:04:45 -0500 Subject: [PATCH 08/18] up --- src/systems/diffeqs/abstractodesystem.jl | 6 ++-- test/bvproblem.jl | 44 ++++++++++++------------ 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 112419dea8..b795bb981d 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -529,12 +529,12 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] kwargs1 = merge(kwargs1, (callback = cbs,)) end - # Construct initial conditions + # Construct initial conditions. _u0 = u0 isa Function ? u0(tspan[1]) : u0 - # Define the boundary conditions + # Define the boundary conditions. bc = if iip - (residual, u, p, t) -> residual .= u[1] - _u0 + (residual, u, p, t) -> (residual .= u[1] - _u0) else (u, p, t) -> (u[1] - _u0) end diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 90d41a96ad..4865c867b3 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -2,32 +2,31 @@ using BoundaryValueDiffEq, OrdinaryDiffEq using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D -@parameters σ = 10. ρ = 28 β = 8/3 -@variables x(t) = 1 y(t) = 0 z(t) = 0 +@parameters α = 7.5 β = 4. γ = 8. δ = 5. +@variables x(t) = 1. y(t) = 2. -eqs = [D(x) ~ σ*(y-x), - D(y) ~ x*(ρ-z)-y, - D(z) ~ x*y - β*z] +eqs = [D(x) ~ α*x - β*x*y, + D(y) ~ -γ*y + δ*x*y] -u0map = [:x => 1., :y => 0., :z => 0.] -parammap = [:ρ => 28., :β => 8/3, :σ => 10.] +u0map = [:x => 1., :y => 2.] +parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.] tspan = (0., 10.) -@mtkbuild lorenz = ODESystem(eqs, t) +@mtkbuild lotkavolterra = ODESystem(eqs, t) -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.1); +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) +sol = solve(bvp, MIRK4(), dt = 0.01); -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) -sol2 = solve(bvp, MIRK4(), dt = 0.1); +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) +sol2 = solve(bvp, MIRK4(), dt = 0.01); -op = ODEProblem(lorenz, u0map, tspan, parammap) -osol = solve(op) +op = ODEProblem(lotkavolterra, u0map, tspan, parammap) +osol = solve(op, Vern9()) -@test sol.u[end] ≈ osol.u[end] -@test sol2.u[end] ≈ osol.u[end] -@test sol.u[1] == [1., 0., 0.] -@test sol2.u[1] == [1., 0., 0.] +@test isapprox(sol.u[end],osol.u[end]; atol = 0.001) +@test isapprox(sol2.u[end],osol.u[end]; atol = 0.001) +@test sol.u[1] == [1., 2.] +@test sol2.u[1] == [1., 2.] ### Testing on pendulum @@ -39,16 +38,17 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) u0map = [θ => π/2, D(θ) => π/2] -parammap = [:L => 2.] +parammap = [:L => 2., :g => 9.81] tspan = (0., 10.) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.05); +sol = solve(bvp, MIRK4(), dt = 0.01); bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -sol2 = solve(bvp2, MIRK4(), dt = 0.05); +sol2 = solve(bvp2, MIRK4(), dt = 0.01); -osol = solve(pend) +op = ODEProblem(pend, u0map, tspan, parammap) +osol = solve(op, Vern9()) @test sol.u[end] ≈ osol.u[end] @test sol.u[1] == [π/2, π/2] From a9f210691c1f38ef0f575eb89a939763b81c9c1b Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 3 Dec 2024 19:11:35 -0500 Subject: [PATCH 09/18] up --- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- test/bvproblem.jl | 8 ++++---- test/runtests.jl | 1 + 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index b795bb981d..8dac19f296 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -534,7 +534,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] # Define the boundary conditions. bc = if iip - (residual, u, p, t) -> (residual .= u[1] - _u0) + (residual, u, p, t) -> (residual .= u[1] .- _u0) else (u, p, t) -> (u[1] - _u0) end @@ -542,7 +542,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) end -get_callback(prob::BVProblem) = prob.kwargs[:callback] +get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") """ ```julia diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 4865c867b3..7235638cf0 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -23,8 +23,8 @@ sol2 = solve(bvp, MIRK4(), dt = 0.01); op = ODEProblem(lotkavolterra, u0map, tspan, parammap) osol = solve(op, Vern9()) -@test isapprox(sol.u[end],osol.u[end]; atol = 0.001) -@test isapprox(sol2.u[end],osol.u[end]; atol = 0.001) +@test isapprox(sol.u[end],osol.u[end]; atol = 0.01) +@test isapprox(sol2.u[end],osol.u[end]; atol = 0.01) @test sol.u[1] == [1., 2.] @test sol2.u[1] == [1., 2.] @@ -50,7 +50,7 @@ sol2 = solve(bvp2, MIRK4(), dt = 0.01); op = ODEProblem(pend, u0map, tspan, parammap) osol = solve(op, Vern9()) -@test sol.u[end] ≈ osol.u[end] +@test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [π/2, π/2] -@test sol2.u[end] ≈ osol.u[end] +@test isapprox(sol2.u[end], osol.u[end]; atol = 0.01) @test sol2.u[1] == [π/2, π/2] diff --git a/test/runtests.jl b/test/runtests.jl index 44846eed57..eaa87e6407 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -80,6 +80,7 @@ end @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") + @safetestset "BVProblem Test" include("bvproblem.jl") @safetestset "print_tree" include("print_tree.jl") @safetestset "Constraints Test" include("constraints.jl") end From 18fdd5f79702627e30ce0acc971677a5cd1d303f Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 13 Dec 2024 10:04:59 +0900 Subject: [PATCH 10/18] up --- Project.toml | 6 ++++-- test/bvproblem.jl | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 093c632c42..19063a0342 100644 --- a/Project.toml +++ b/Project.toml @@ -72,14 +72,15 @@ MTKBifurcationKitExt = "BifurcationKit" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKHomotopyContinuationExt = "HomotopyContinuation" -MTKLabelledArraysExt = "LabelledArrays" MTKInfiniteOptExt = "InfiniteOpt" +MTKLabelledArraysExt = "LabelledArrays" [compat] AbstractTrees = "0.3, 0.4" ArrayInterface = "6, 7" BifurcationKit = "0.4" BlockArrays = "1.1" +BoundaryValueDiffEq = "5.12.0" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" @@ -145,6 +146,7 @@ julia = "1.9" [extras] AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" @@ -174,4 +176,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"] diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7235638cf0..7787bd9c3e 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -44,7 +44,7 @@ tspan = (0., 10.) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) sol = solve(bvp, MIRK4(), dt = 0.01); -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) +bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) sol2 = solve(bvp2, MIRK4(), dt = 0.01); op = ODEProblem(pend, u0map, tspan, parammap) From 9d65a3345ade530654a64ace2069ff24e8c67875 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 16 Dec 2024 23:43:35 +0800 Subject: [PATCH 11/18] fixing create_array --- Project.toml | 3 ++ src/systems/diffeqs/abstractodesystem.jl | 10 ++++- test/bvproblem.jl | 53 +++++++++++++++--------- 3 files changed, 45 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 19063a0342..ab854b80bd 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,9 @@ version = "9.54.0" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" +BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a" +BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 8dac19f296..3a502d3183 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -539,11 +539,19 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] (u, p, t) -> (u[1] - _u0) end - return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) + return BVProblem{iip}(f, bc, _u0, tspan, p; kwargs1..., kwargs...) end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") +@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where dims + [elems...] +end + +@inline function create_array(::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where dims + T[elems...] +end + """ ```julia DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys), diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7787bd9c3e..7cdb7e1837 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -2,6 +2,8 @@ using BoundaryValueDiffEq, OrdinaryDiffEq using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D +solvers = [MIRK4, RadauIIa5, LobattoIIIa3] + @parameters α = 7.5 β = 4. γ = 8. δ = 5. @variables x(t) = 1. y(t) = 2. @@ -13,20 +15,26 @@ parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.] tspan = (0., 10.) @mtkbuild lotkavolterra = ODESystem(eqs, t) +op = ODEProblem(lotkavolterra, u0map, tspan, parammap) +osol = solve(op, Vern9()) -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.01); +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) -sol2 = solve(bvp, MIRK4(), dt = 0.01); +for solver in solvers + println("$solver") + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1., 2.] +end -op = ODEProblem(lotkavolterra, u0map, tspan, parammap) -osol = solve(op, Vern9()) +# Test out of place +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) -@test isapprox(sol.u[end],osol.u[end]; atol = 0.01) -@test isapprox(sol2.u[end],osol.u[end]; atol = 0.01) -@test sol.u[1] == [1., 2.] -@test sol2.u[1] == [1., 2.] +for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + @test sol.u[1] == [1., 2.] +end ### Testing on pendulum @@ -38,19 +46,24 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) u0map = [θ => π/2, D(θ) => π/2] -parammap = [:L => 2., :g => 9.81] +parammap = [:L => 1., :g => 9.81] tspan = (0., 10.) +op = ODEProblem(pend, u0map, tspan, parammap) +osol = solve(op, Vern9()) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.01); +for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + @test sol.u[1] == [π/2, π/2] +end +# Test out-of-place bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) -sol2 = solve(bvp2, MIRK4(), dt = 0.01); - -op = ODEProblem(pend, u0map, tspan, parammap) -osol = solve(op, Vern9()) -@test isapprox(sol.u[end], osol.u[end]; atol = 0.01) -@test sol.u[1] == [π/2, π/2] -@test isapprox(sol2.u[end], osol.u[end]; atol = 0.01) -@test sol2.u[1] == [π/2, π/2] +for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + @test sol.u[1] == [π/2, π/2] +end From 999ec308d58e32c51941743a6ba5a8f096c56ac2 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 16 Dec 2024 23:44:24 +0800 Subject: [PATCH 12/18] revert Project.toml --- Project.toml | 3 --- 1 file changed, 3 deletions(-) diff --git a/Project.toml b/Project.toml index ab854b80bd..19063a0342 100644 --- a/Project.toml +++ b/Project.toml @@ -7,9 +7,6 @@ version = "9.54.0" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" -BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a" -BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" From 9226ad687ff06cbc1dbbdcaa16ba3df85ef32d4e Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 16 Dec 2024 23:56:04 +0800 Subject: [PATCH 13/18] Up --- test/bvproblem.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7cdb7e1837..86e3722eec 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -21,7 +21,6 @@ osol = solve(op, Vern9()) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) for solver in solvers - println("$solver") sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [1., 2.] @@ -47,15 +46,15 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] u0map = [θ => π/2, D(θ) => π/2] parammap = [:L => 1., :g => 9.81] -tspan = (0., 10.) +tspan = (0., 6.) op = ODEProblem(pend, u0map, tspan, parammap) osol = solve(op, Vern9()) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) for solver in solvers - sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [π/2, π/2] end From 67d8164c591b74a9b985ffbcaf1ef248fb0efaaa Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 00:09:52 +0800 Subject: [PATCH 14/18] formatting --- src/systems/diffeqs/abstractodesystem.jl | 11 ++++--- test/bvproblem.jl | 42 +++++++++++++----------- 2 files changed, 28 insertions(+), 25 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 87ee6e823d..06c83073bf 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -512,7 +512,6 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip, specialize} - if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") end @@ -528,12 +527,12 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] if cbs !== nothing kwargs1 = merge(kwargs1, (callback = cbs,)) end - + # Construct initial conditions. _u0 = u0 isa Function ? u0(tspan[1]) : u0 # Define the boundary conditions. - bc = if iip + bc = if iip (residual, u, p, t) -> (residual .= u[1] .- _u0) else (u, p, t) -> (u[1] - _u0) @@ -544,11 +543,13 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where dims +@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, + ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function create_array(::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where dims +@inline function create_array( + ::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 86e3722eec..1072874917 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -4,49 +4,51 @@ using ModelingToolkit: t_nounits as t, D_nounits as D solvers = [MIRK4, RadauIIa5, LobattoIIIa3] -@parameters α = 7.5 β = 4. γ = 8. δ = 5. -@variables x(t) = 1. y(t) = 2. +@parameters α=7.5 β=4.0 γ=8.0 δ=5.0 +@variables x(t)=1.0 y(t)=2.0 -eqs = [D(x) ~ α*x - β*x*y, - D(y) ~ -γ*y + δ*x*y] +eqs = [D(x) ~ α * x - β * x * y, + D(y) ~ -γ * y + δ * x * y] -u0map = [:x => 1., :y => 2.] -parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.] -tspan = (0., 10.) +u0map = [:x => 1.0, :y => 2.0] +parammap = [:α => 7.5, :β => 4, :γ => 8.0, :δ => 5.0] +tspan = (0.0, 10.0) @mtkbuild lotkavolterra = ODESystem(eqs, t) op = ODEProblem(lotkavolterra, u0map, tspan, parammap) osol = solve(op, Vern9()) -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) for solver in solvers sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [1., 2.] + @test sol.u[1] == [1.0, 2.0] end # Test out of place -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) for solver in solvers sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) - @test sol.u[1] == [1., 2.] + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1.0, 2.0] end ### Testing on pendulum -@parameters g = 9.81 L = 1. -@variables θ(t) = π/2 +@parameters g=9.81 L=1.0 +@variables θ(t) = π / 2 eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) -u0map = [θ => π/2, D(θ) => π/2] -parammap = [:L => 1., :g => 9.81] -tspan = (0., 6.) +u0map = [θ => π / 2, D(θ) => π / 2] +parammap = [:L => 1.0, :g => 9.81] +tspan = (0.0, 6.0) op = ODEProblem(pend, u0map, tspan, parammap) osol = solve(op, Vern9()) @@ -55,7 +57,7 @@ bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, pa for solver in solvers sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [π/2, π/2] + @test sol.u[1] == [π / 2, π / 2] end # Test out-of-place @@ -63,6 +65,6 @@ bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, for solver in solvers sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) - @test sol.u[1] == [π/2, π/2] + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [π / 2, π / 2] end From 25988f3bc6b6d66b008b6a40341f75e4547e574a Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 10:52:37 +0800 Subject: [PATCH 15/18] up --- src/systems/diffeqs/abstractodesystem.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 06c83073bf..c84f5ff5be 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -543,13 +543,13 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, +@inline function SciML.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function create_array( - ::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} +@inline function SciML.Code.create_array( + ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end From bb28d4fe2a0d753221104186fe42808c7657f5d6 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 10:53:22 +0800 Subject: [PATCH 16/18] up --- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index c84f5ff5be..eac2df16aa 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -543,12 +543,12 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function SciML.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, +@inline function SciMLBase.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function SciML.Code.create_array( +@inline function SciMLBase.Code.create_array( ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end From b2bf7c05532bdf68d93e566d07c5e7444473dd7a Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 11:00:30 +0800 Subject: [PATCH 17/18] fix --- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index eac2df16aa..3d143d49fb 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -543,12 +543,12 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function SciMLBase.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, +@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function SciMLBase.Code.create_array( +@inline function SymbolicUtils.Code.create_array( ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end From 3751c2a92c282593ef52b67b727c7db635ea677e Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 20 Dec 2024 23:58:28 +0900 Subject: [PATCH 18/18] up --- src/systems/diffeqs/abstractodesystem.jl | 165 +++++++++++------------ test/bvproblem.jl | 4 +- 2 files changed, 83 insertions(+), 86 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 3d143d49fb..d8ff71c324 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -469,90 +469,6 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end -""" -```julia -SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - jac = true, sparse = true, - simplify = false, - kwargs...) where {iip} -``` - -Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and -`ps` are used to set the order of the dependent variable and parameter vectors, -respectively. `u0map` should be used to specify the initial condition, or be a function returning an initial condition. -""" -function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) - BVProblem{true}(sys, args...; kwargs...) -end - -function SciMLBase.BVProblem(sys::AbstractODESystem, - u0map::StaticArray, - args...; - kwargs...) - BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) -end - -function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) - BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) - BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], - tspan = get_tspan(sys), - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - callback = nothing, - check_length = true, - warn_initialize_determined = true, - eval_expression = false, - eval_module = @__MODULE__, - kwargs...) where {iip, specialize} - if !iscomplete(sys) - error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") - end - - f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; - t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) - - cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) - kwargs = filter_kwargs(kwargs) - - kwargs1 = (;) - if cbs !== nothing - kwargs1 = merge(kwargs1, (callback = cbs,)) - end - - # Construct initial conditions. - _u0 = u0 isa Function ? u0(tspan[1]) : u0 - - # Define the boundary conditions. - bc = if iip - (residual, u, p, t) -> (residual .= u[1] .- _u0) - else - (u, p, t) -> (u[1] - _u0) - end - - return BVProblem{iip}(f, bc, _u0, tspan, p; kwargs1..., kwargs...) -end - -get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") - -@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, - ::Val{1}, ::Val{dims}, elems...) where {dims} - [elems...] -end - -@inline function SymbolicUtils.Code.create_array( - ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} - T[elems...] -end - """ ```julia DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys), @@ -943,6 +859,87 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = end get_callback(prob::ODEProblem) = prob.kwargs[:callback] +""" +```julia +SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` + +Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and +`ps` are used to set the order of the dependent variable and parameter vectors, +respectively. `u0map` should be used to specify the initial condition. +""" +function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem(sys::AbstractODESystem, + u0map::StaticArray, + args...; + kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) +end + +function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], + tspan = get_tspan(sys), + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + callback = nothing, + check_length = true, + warn_initialize_determined = true, + eval_expression = false, + eval_module = @__MODULE__, + kwargs...) where {iip, specialize} + if !iscomplete(sys) + error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") + end + + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, + check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) + + cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) + kwargs = filter_kwargs(kwargs) + + kwargs1 = (;) + if cbs !== nothing + kwargs1 = merge(kwargs1, (callback = cbs,)) + end + + # Define the boundary conditions. + bc = if iip + (residual, u, p, t) -> (residual .= u[1] .- u0) + else + (u, p, t) -> (u[1] - u0) + end + + return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) +end + +get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") + +@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, + ::Val{1}, ::Val{dims}, elems...) where {dims} + [elems...] +end + +@inline function SymbolicUtils.Code.create_array( + ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} + T[elems...] +end + """ ```julia DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 1072874917..c5a302147d 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -10,8 +10,8 @@ solvers = [MIRK4, RadauIIa5, LobattoIIIa3] eqs = [D(x) ~ α * x - β * x * y, D(y) ~ -γ * y + δ * x * y] -u0map = [:x => 1.0, :y => 2.0] -parammap = [:α => 7.5, :β => 4, :γ => 8.0, :δ => 5.0] +u0map = [x => 1.0, y => 2.0] +parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] tspan = (0.0, 10.0) @mtkbuild lotkavolterra = ODESystem(eqs, t)