From 0143c3f3b10571dfeea834261520b9f3626df18e Mon Sep 17 00:00:00 2001 From: Qingyu Qu <2283984853@qq.com> Date: Tue, 28 Jan 2025 01:15:59 +0800 Subject: [PATCH] Change ODEIntreface to BVPInterface in extension --- Project.toml | 10 +- ... => BoundaryValueDiffEqBVPInterfaceExt.jl} | 175 +++++++++++++++++- src/BoundaryValueDiffEq.jl | 2 +- src/extension_algs.jl | 96 ++++++++-- test/wrappers/odeinterface_tests.jl | 40 ++-- 5 files changed, 282 insertions(+), 41 deletions(-) rename ext/{BoundaryValueDiffEqODEInterfaceExt.jl => BoundaryValueDiffEqBVPInterfaceExt.jl} (68%) diff --git a/Project.toml b/Project.toml index e265dcf2..649a667c 100644 --- a/Project.toml +++ b/Project.toml @@ -19,10 +19,10 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" [weakdeps] -ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c" +BVPInterface = "349b9a96-74d8-4082-8332-dd34ffd0a01f" [extensions] -BoundaryValueDiffEqODEInterfaceExt = "ODEInterface" +BoundaryValueDiffEqBVPInterfaceExt = "BVPInterface" [compat] ADTypes = "1.11" @@ -42,7 +42,7 @@ Hwloc = "3.3" InteractiveUtils = "<0.0.1, 1" JET = "0.9.12" LinearAlgebra = "1.10" -ODEInterface = "0.5" +BVPInterface = "0.1.0" OrdinaryDiffEq = "6.90.1" Pkg = "1.10.0" Random = "1.10" @@ -61,7 +61,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d" -ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c" +BVPInterface = "349b9a96-74d8-4082-8332-dd34ffd0a01f" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -71,4 +71,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "DiffEqDevTools", "Hwloc", "InteractiveUtils", "JET", "LinearSolve", "NonlinearSolveFirstOrder", "ODEInterface", "OrdinaryDiffEq", "Pkg", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"] +test = ["Aqua", "DiffEqDevTools", "Hwloc", "InteractiveUtils", "JET", "LinearSolve", "NonlinearSolveFirstOrder", "BVPInterface", "OrdinaryDiffEq", "Pkg", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"] diff --git a/ext/BoundaryValueDiffEqODEInterfaceExt.jl b/ext/BoundaryValueDiffEqBVPInterfaceExt.jl similarity index 68% rename from ext/BoundaryValueDiffEqODEInterfaceExt.jl rename to ext/BoundaryValueDiffEqBVPInterfaceExt.jl index 61c580b7..7cffe0f2 100644 --- a/ext/BoundaryValueDiffEqODEInterfaceExt.jl +++ b/ext/BoundaryValueDiffEqBVPInterfaceExt.jl @@ -1,18 +1,19 @@ -module BoundaryValueDiffEqODEInterfaceExt +module BoundaryValueDiffEqBVPInterfaceExt -using BoundaryValueDiffEq: BVPM2, BVPSOL, COLNEW +using BoundaryValueDiffEq: BVPM2, BVPSOL, COLSYS, COLNEW using BoundaryValueDiffEqCore: BoundaryValueDiffEqAlgorithm, __extract_u0, __initial_guess_length, __extract_mesh, __flatten_initial_guess, __get_bcresid_prototype, __has_initial_guess, __initial_guess using SciMLBase: SciMLBase, BVProblem, TwoPointBVProblem, ReturnCode -using ODEInterface: OptionsODE, OPT_ATOL, OPT_RTOL, OPT_METHODCHOICE, OPT_DIAGNOSTICOUTPUT, +using BVPInterface: OptionsODE, OPT_ATOL, OPT_RTOL, OPT_METHODCHOICE, OPT_DIAGNOSTICOUTPUT, OPT_ERRORCONTROL, OPT_SINGULARTERM, OPT_MAXSTEPS, OPT_BVPCLASS, OPT_SOLMETHOD, OPT_RHS_CALLMODE, OPT_COLLOCATIONPTS, OPT_ADDGRIDPOINTS, OPT_MAXSUBINTERVALS, RHS_CALL_INSITU, evalSolution -using ODEInterface: Bvpm2, bvpm2_init, bvpm2_solve, bvpm2_destroy, bvpm2_get_x -using ODEInterface: bvpsol -using ODEInterface: colnew +using BVPInterface: Bvpm2, bvpm2_init, bvpm2_solve, bvpm2_destroy, bvpm2_get_x +using BVPInterface: bvpsol +using BVPInterface: colsys +using BVPInterface: colnew using FastClosures: @closure using ForwardDiff: ForwardDiff @@ -339,6 +340,168 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, stats = destats, original = (sol, retcode, stats)) end +#------- +# COLSYS +#------- +function SciMLBase.__solve(prob::BVProblem, alg::COLSYS; maxiters = 1000, + reltol = 1e-3, dt = 0.0, verbose = true, kwargs...) + dt ≤ 0 && throw(ArgumentError("`dt` must be positive")) + + t₀, t₁ = prob.tspan + u0_ = __extract_u0(prob.u0, prob.p, t₀) + u0_size = size(u0_) + n = __initial_guess_length(prob.u0) + + u0 = __flatten_initial_guess(prob.u0) + mesh = __extract_mesh(prob.u0, t₀, t₁, ifelse(n == -1, dt, n - 1)) + if u0 === nothing + # initial_guess function was provided + u0 = mapreduce(@closure(t->vec(__initial_guess(prob.u0, prob.p, t))), hcat, mesh) + end + + no_odes = length(u0_) + + # has_initial_guess = prob.u0 isa AbstractVector{<:AbstractArray} + # dt ≤ 0 && throw(ArgumentError("dt must be positive")) + # no_odes, n, u0 = if has_initial_guess + # length(first(prob.u0)), (length(prob.u0) - 1), reduce(hcat, prob.u0) + # else + # length(prob.u0), Int(cld((prob.tspan[2] - prob.tspan[1]), dt)), prob.u0 + # end + + T = eltype(u0) + # mesh = collect(range(prob.tspan[1], stop = prob.tspan[2], length = n + 1)) + orders = ones(Int, no_odes) + _tspan = [prob.tspan[1], prob.tspan[2]] + iip = SciMLBase.isinplace(prob) + + rhs = @closure (t, u, du) -> begin + if iip + prob.f(du, u, prob.p, t) + else + (du .= prob.f(u, prob.p, t)) + end + end + + if prob.f.jac === nothing + if iip + jac = (df, u, p, t) -> begin + _du = similar(u) + prob.f(_du, u, p, t) + _f = @closure (du, u) -> prob.f(du, u, p, t) + ForwardDiff.jacobian!(df, _f, _du, u) + return + end + else + jac = (df, u, p, t) -> begin + _du = prob.f(u, p, t) + _f = @closure (du, u) -> (du .= prob.f(u, p, t)) + ForwardDiff.jacobian!(df, _f, _du, u) + return + end + end + else + jac = prob.f.jac + end + Drhs = @closure (t, u, df) -> jac(df, u, prob.p, t) + + bcresid_prototype, _ = __get_bcresid_prototype(prob.problem_type, prob, u0) + + if prob.problem_type isa TwoPointBVProblem + n_bc_a = length(first(bcresid_prototype)) + n_bc_b = length(last(bcresid_prototype)) + zeta = vcat(fill(first(prob.tspan), n_bc_a), fill(last(prob.tspan), n_bc_b)) + bc = @closure (i, z, resid) -> begin + tmpa = copy(z) + tmpb = copy(z) + tmp_resid_a = zeros(T, n_bc_a) + tmp_resid_b = zeros(T, n_bc_b) + prob.f.bc[1](tmp_resid_a, tmpa, prob.p) + prob.f.bc[2](tmp_resid_b, tmpb, prob.p) + + for j in 1:n_bc_a + if i == j + resid[1] = tmp_resid_a[j] + end + end + for j in 1:n_bc_b + if i == (j + n_bc_a) + resid[1] = tmp_resid_b[j] + end + end + end + + Dbc = @closure (i, z, dbc) -> begin + for j in 1:n_bc_a + if i == j + dbc[i] = 1.0 + end + end + for j in 1:n_bc_b + if i == (j + n_bc_a) + dbc[i] = 1.0 + end + end + end + fixed_points = nothing + else + zeta = sort(alg.zeta) + bc = alg.bc_func + Dbc = alg.dbc_func + left_index = findlast(x -> x ≈ t₀, zeta) + 1 + right_index = findfirst(x -> x ≈ t₁, zeta) - 1 + fixed_points = alg.zeta[left_index:right_index] + end + + if fixed_points === nothing + opt = OptionsODE( + OPT_BVPCLASS => alg.bvpclass, OPT_COLLOCATIONPTS => alg.collocationpts, + OPT_MAXSTEPS => maxiters, OPT_DIAGNOSTICOUTPUT => alg.diagnostic_output, + OPT_MAXSUBINTERVALS => alg.max_num_subintervals, OPT_RTOL => reltol) + else + opt = OptionsODE( + OPT_BVPCLASS => alg.bvpclass, OPT_COLLOCATIONPTS => alg.collocationpts, + OPT_MAXSTEPS => maxiters, OPT_DIAGNOSTICOUTPUT => alg.diagnostic_output, + OPT_MAXSUBINTERVALS => alg.max_num_subintervals, + OPT_RTOL => reltol, OPT_ADDGRIDPOINTS => fixed_points) + end + + # Provide initial guess could really help COLSYS to converge + guess = @closure (t, u, du) -> begin + if iip + u .= prob.u0 + prob.f(du, u, prob.p, t) + else + u .= prob.u0 + (du .= prob.f(u, prob.p, t)) + end + end + + sol, retcode, stats = colsys(_tspan, orders, zeta, rhs, Drhs, bc, Dbc, guess, opt) + + if verbose + if retcode == 0 + @warn "Collocation matrix is singular" + elseif retcode == -1 + @warn "The expected no. of subintervals exceeds storage(try to increase \ + `OPT_MAXSUBINTERVALS`)" + elseif retcode == -2 + @warn "The nonlinear iteration has not converged" + elseif retcode == -3 + @warn "There is an input data error" + end + end + + evalsol = evalSolution(sol, mesh) + destats = SciMLBase.DEStats( + stats["no_rhs_calls"], 0, 0, 0, stats["no_jac_calls"], 0, 0, 0, 0, 0, 0, 0, 0) + + return SciMLBase.build_solution( + prob, alg, mesh, collect(Vector{eltype(evalsol)}, eachrow(evalsol)); + retcode = retcode > 0 ? ReturnCode.Success : ReturnCode.Failure, + stats = destats, original = (sol, retcode, stats)) +end + export BVPM2, BVPSOL, COLNEW end diff --git a/src/BoundaryValueDiffEq.jl b/src/BoundaryValueDiffEq.jl index 830720cb..93206a03 100644 --- a/src/BoundaryValueDiffEq.jl +++ b/src/BoundaryValueDiffEq.jl @@ -31,7 +31,7 @@ export MIRKN4, MIRKN6 export Ascher1, Ascher2, Ascher3, Ascher4, Ascher5, Ascher6, Ascher7 -export BVPM2, BVPSOL, COLNEW # From ODEInterface.jl +export BVPM2, BVPSOL, COLSYS, COLNEW # From ODEInterface.jl export MIRKJacobianComputationAlgorithm, BVPJacobianAlgorithm diff --git a/src/extension_algs.jl b/src/extension_algs.jl index b6a8434a..a956de32 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -1,4 +1,4 @@ -# Algorithms from ODEInterface.jl +# Algorithms from BVPInterface.jl """ BVPM2(; max_num_subintervals = 3000, method_choice = 4, diagnostic_output = 1, error_control = 1, singular_term = nothing) @@ -33,7 +33,7 @@ Fortran code for solving two-point boundary value problems. For detailed documen !!! note - Only available if the `ODEInterface` package is loaded. + Only available if the `BVPInterface` package is loaded. """ struct BVPM2{S} <: BoundaryValueDiffEqAlgorithm max_num_subintervals::Int @@ -44,8 +44,8 @@ struct BVPM2{S} <: BoundaryValueDiffEqAlgorithm function BVPM2(max_num_subintervals::Int, method_choice::Int, diagnostic_output::Int, error_control::Int, singular_term::Union{Nothing, AbstractMatrix}) - if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqODEInterfaceExt) === nothing - error("`BVPM2` requires `ODEInterface.jl` to be loaded") + if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqBVPInterfaceExt) === nothing + error("`BVPM2` requires `BVPInterface.jl` to be loaded") end return new{typeof(singular_term)}(max_num_subintervals, method_choice, diagnostic_output, error_control, singular_term) @@ -85,7 +85,7 @@ For detailed documentation, see !!! note - Only available if the `ODEInterface` package is loaded. + Only available if the `BVPInterface` package is loaded. """ struct BVPSOL{O} <: BoundaryValueDiffEqAlgorithm bvpclass::Int @@ -93,8 +93,8 @@ struct BVPSOL{O} <: BoundaryValueDiffEqAlgorithm odesolver::O function BVPSOL(bvpclass::Int, sol_method::Int, odesolver) - if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqODEInterfaceExt) === nothing - error("`BVPSOL` requires `ODEInterface.jl` to be loaded") + if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqBVPInterfaceExt) === nothing + error("`BVPSOL` requires `BVPInterface.jl` to be loaded") end return new{typeof(odesolver)}(bvpclass, sol_method, odesolver) end @@ -104,6 +104,78 @@ function BVPSOL(; bvpclass = 2, sol_method = 0, odesolver = nothing) return BVPSOL(bvpclass, sol_method, odesolver) end +""" + COLSYS(; bvpclass = 1, collocationpts = 7, diagnostic_output = 1, + max_num_subintervals = 3000, bc_func = nothing, dbc_func = nothing, + zeta = nothing) + COLSYS(bvpclass::Int, collocationpts::Int, diagnostic_output::Int, + max_num_subintervals::Int, bc_func, dbc_func, zeta::AbstractVector) + +## Keyword Arguments: + + - `bvpclass`: Boundary value problem classification, default as nonlinear and + "extra sensitive", available choices: + + + `0`: Linear boundary value problem. + + `1`: Nonlinear and regular. + + `2`: Nonlinear and "extra sensitive" (first relax factor is rstart and the + nonlinear iteration does not rely on past convergence). + + `3`: fail-early: return immediately upon: a) two successive non-convergences. + b) after obtaining an error estimate for the first time. + + - `collocationpts`: Number of collocation points per subinterval. Require + orders[i] ≤ k ≤ 7, default as 7 + - `diagnostic_output`: Diagnostic output for COLSYS, default as no printout, available + choices: + + + `-1`: Full diagnostic printout. + + `0`: Selected printout. + + `1`: No printout. + - `max_num_subintervals`: Number of maximal subintervals, default as 3000. + - `bc_func`: Boundary condition accord with BVPInterface.jl, only used for multi-points BVP. + - `dbc_func`: Jacobian of boundary condition accord with BVPInterface.jl, only used for multi-points BVP. + - `zeta`: The points in interval where boundary conditions are specified, only used for multi-points BVP. + +A Fortran77 code solves a multi-points boundary value problems for a mixed order system of +ODEs. It incorporates a new basis representation replacing b-splines, and improvements for +the linear and nonlinear algebraic equation solvers. + +!!! warning + + Only supports two-point boundary value problems. + +!!! note + + Only available if the `BVPInterface` package is loaded. +""" +struct COLSYS <: BoundaryValueDiffEqAlgorithm + bvpclass::Int + collocationpts::Int + diagnostic_output::Int + max_num_subintervals::Int + bc_func::Union{Function, Nothing} + dbc_func::Union{Function, Nothing} + zeta::Union{AbstractVector, Nothing} + + function COLSYS(bvpclass::Int, collocationpts::Int, diagnostic_output::Int, + max_num_subintervals::Int, bc_func::Union{Function, Nothing}, + dbc_func::Union{Function, Nothing}, zeta::Union{AbstractVector, Nothing}) + if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqBVPInterfaceExt) === nothing + error("`COLSYS` requires `BVPInterface.jl` to be loaded") + end + return new(bvpclass, collocationpts, diagnostic_output, + max_num_subintervals, bc_func, dbc_func, zeta) + end +end + +function COLSYS(; bvpclass::Int = 1, collocationpts::Int = 7, diagnostic_output::Int = 1, + max_num_subintervals::Int = 4000, bc_func::Union{Function, Nothing} = nothing, + dbc_func::Union{Function, Nothing} = nothing, + zeta::Union{AbstractVector, Nothing} = nothing) + return COLSYS(bvpclass, collocationpts, diagnostic_output, + max_num_subintervals, bc_func, dbc_func, zeta) +end + """ COLNEW(; bvpclass = 1, collocationpts = 7, diagnostic_output = 1, max_num_subintervals = 3000, bc_func = nothing, dbc_func = nothing, @@ -132,8 +204,8 @@ end + `0`: Selected printout. + `1`: No printout. - `max_num_subintervals`: Number of maximal subintervals, default as 3000. - - `bc_func`: Boundary condition accord with ODEInterface.jl, only used for multi-points BVP. - - `dbc_func`: Jacobian of boundary condition accord with ODEInterface.jl, only used for multi-points BVP. + - `bc_func`: Boundary condition accord with BVPInterface.jl, only used for multi-points BVP. + - `dbc_func`: Jacobian of boundary condition accord with BVPInterface.jl, only used for multi-points BVP. - `zeta`: The points in interval where boundary conditions are specified, only used for multi-points BVP. A Fortran77 code solves a multi-points boundary value problems for a mixed order system of @@ -146,7 +218,7 @@ the linear and nonlinear algebraic equation solvers. !!! note - Only available if the `ODEInterface` package is loaded. + Only available if the `BVPInterface` package is loaded. """ struct COLNEW <: BoundaryValueDiffEqAlgorithm bvpclass::Int @@ -160,8 +232,8 @@ struct COLNEW <: BoundaryValueDiffEqAlgorithm function COLNEW(bvpclass::Int, collocationpts::Int, diagnostic_output::Int, max_num_subintervals::Int, bc_func::Union{Function, Nothing}, dbc_func::Union{Function, Nothing}, zeta::Union{AbstractVector, Nothing}) - if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqODEInterfaceExt) === nothing - error("`COLNEW` requires `ODEInterface.jl` to be loaded") + if Base.get_extension(@__MODULE__, :BoundaryValueDiffEqBVPInterfaceExt) === nothing + error("`COLNEW` requires `BVPInterface.jl` to be loaded") end return new(bvpclass, collocationpts, diagnostic_output, max_num_subintervals, bc_func, dbc_func, zeta) diff --git a/test/wrappers/odeinterface_tests.jl b/test/wrappers/odeinterface_tests.jl index de5a075d..0a6fce82 100644 --- a/test/wrappers/odeinterface_tests.jl +++ b/test/wrappers/odeinterface_tests.jl @@ -1,8 +1,8 @@ -@testsetup module ODEInterfaceWrapperTestSetup +@testsetup module BVPInterfaceWrapperTestSetup -using BoundaryValueDiffEq, LinearAlgebra, ODEInterface, Random, RecursiveArrayTools +using BoundaryValueDiffEq, LinearAlgebra, BVPInterface, Random, RecursiveArrayTools -# Adaptation of https://github.com/luchr/ODEInterface.jl/blob/958b6023d1dabf775033d0b89c5401b33100bca3/examples/BasicExamples/ex7.jl +# Adaptation of https://github.com/ErikQQY/BVPInterface.jl/blob/master/examples/colnew_example.jl function ex7_f!(du, u, p, t) ϵ = p[1] u₁, λ = u @@ -29,23 +29,23 @@ export ex7_f!, ex7_2pbc1!, ex7_2pbc2!, u0, p, tspan end -@testitem "BVPM2" setup=[ODEInterfaceWrapperTestSetup] begin - using ODEInterface, RecursiveArrayTools, LinearAlgebra +@testitem "BVPM2" setup=[BVPInterfaceWrapperTestSetup] begin + using BVPInterface, RecursiveArrayTools, LinearAlgebra tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, p; bcresid_prototype = (zeros(1), zeros(1))) - #sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) - #@test SciMLBase.successful_retcode(sol_bvpm2) - #resid_f = (Array{Float64, 1}(undef, 1), Array{Float64, 1}(undef, 1)) - #ex7_2pbc1!(resid_f[1], sol_bvpm2(tspan[1]), nothing) - #ex7_2pbc2!(resid_f[2], sol_bvpm2(tspan[2]), nothing) - #@test norm(resid_f, Inf) < 1e-6 + sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) + @test SciMLBase.successful_retcode(sol_bvpm2) + resid_f = (Array{Float64, 1}(undef, 1), Array{Float64, 1}(undef, 1)) + ex7_2pbc1!(resid_f[1], sol_bvpm2(tspan[1]), nothing) + ex7_2pbc2!(resid_f[2], sol_bvpm2(tspan[2]), nothing) + @test norm(resid_f, Inf) < 1e-6 end # Just test that it runs. BVPSOL only works with linearly separable BCs. -@testitem "BVPSOL" setup=[ODEInterfaceWrapperTestSetup] begin - using ODEInterface, OrdinaryDiffEq, RecursiveArrayTools, NonlinearSolveFirstOrder +@testitem "BVPSOL" setup=[BVPInterfaceWrapperTestSetup] begin + using BVPInterface, OrdinaryDiffEq, RecursiveArrayTools, NonlinearSolveFirstOrder tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, p; bcresid_prototype = (zeros(1), zeros(1))) @@ -86,8 +86,8 @@ end @test sol_bvpsol isa SciMLBase.ODESolution end -@testitem "COLNEW" setup=[ODEInterfaceWrapperTestSetup] begin - using ODEInterface, RecursiveArrayTools +@testitem "COLSYS/COLNEW" setup=[BVPInterfaceWrapperTestSetup] begin + using BVPInterface, RecursiveArrayTools function f!(du, u, p, t) du[1] = u[2] @@ -107,10 +107,13 @@ end prob = TwoPointBVProblem(fun, [1.0, 0.0], tspan) sol_colnew = solve(prob, COLNEW(), dt = 0.01) @test SciMLBase.successful_retcode(sol_colnew) + + sol_colsys = solve(prob, COLSYS(), dt = 0.01) + @test SciMLBase.successful_retcode(sol_colsys) end -@testitem "COLNEW for multi-points BVP" setup=[ODEInterfaceWrapperTestSetup] begin - using ODEInterface, RecursiveArrayTools +@testitem "COLSYS/COLNEW for multi-points BVP" setup=[BVPInterfaceWrapperTestSetup] begin + using BVPInterface, RecursiveArrayTools function f!(du, u, p, t) du[1] = u[2] @@ -155,4 +158,7 @@ end prob = BVProblem(f!, fakebc!, [1.0, 1.0, 0.0], tspan) sol_colnew = solve(prob, COLNEW(bc_func = bc, dbc_func = dbc, zeta = zeta), dt = 0.01) @test SciMLBase.successful_retcode(sol_colnew) + + sol_colsys = solve(prob, COLSYS(bc_func = bc, dbc_func = dbc, zeta = zeta), dt = 0.01) + @test SciMLBase.successful_retcode(sol_colsys) end