up

Author: vyudu

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

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]