Add maxsol and minsol

docs/pages.jl

@@ -2,7 +2,8 @@
 
 pages = ["index.md",
     "Getting Started with BVP solving in Julia" => "tutorials/getting_started.md",
-    "Tutorials" => Any["tutorials/continuation.md", "tutorials/solve_nlls_bvp.md"],
+    "Tutorials" => Any["tutorials/continuation.md",
+        "tutorials/solve_nlls_bvp.md", "tutorials/inequality.md"],
     "Basics" => Any["basics/bvp_problem.md", "basics/bvp_functions.md",
         "basics/solve.md", "basics/autodiff.md", "basics/error_control.md"],
     "Solver Summaries and Recommendations" => Any[

docs/src/tutorials/inequality.md

@@ -0,0 +1,29 @@
+# Solve BVP with Inequality Boundary Conditions
+
+When dealing with scenarios that boundary conditions are stronger than just explicit values at specific points—such as inequality boundary conditions, where the solution must satisfy constraints like staying within upper and lower bounds—we need a more flexible approach. In such cases, we can use the `maxsol` and `minsol` functions provided by BoundaryValueDiffEq.jl to specify such inequality conditions.
+
+Let's walk through this functionlity with an intuitive example. We still revisit the simple pendulum example here, but this time, we’ll impose upper and lower bound constraints on the solution, specified as:
+
+```math
+lb \leq u \leq ub
+```
+
+where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must be bigger than `lb` but smaller than `ub`. To solve such problems, we can simply use the `minsol` and `maxsol` functions when defining the boundary value problem in BoundaryValueDiffEq.jl.
+
+```@example inequality
+using BoundaryValueDiffEq, Plots
+tspan = (0.0, pi / 2)
+function simplependulum!(du, u, p, t)
+    θ = u[1]
+    dθ = u[2]
+    du[1] = dθ
+    du[2] = -9.81 * sin(θ)
+end
+function bc!(residual, u, p, t)
+    residual[1] = maxsol(u, (0.0, pi / 2)) - 5.0496477654230745
+    residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
+end
+prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
+sol = solve(prob, MIRK4(), dt = 0.05)
+plot(sol)
+```

lib/BoundaryValueDiffEqCore/src/solution_utils.jl

@@ -29,3 +29,6 @@ Base.eltype(e::EvalSol) = eltype(e.u)
 function Base.show(io::IO, m::MIME"text/plain", e::EvalSol)
     (print(io, "t: "); show(io, m, e.t); println(io); print(io, "u: "); show(io, m, e.u))
 end
+
+Base.maximum(sol::EvalSol) = maximum(Iterators.flatten(sol.u))
+Base.minimum(sol::EvalSol) = minimum(Iterators.flatten(sol.u))

lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl

@@ -160,5 +160,6 @@ end
 
 export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6, MIRK6I
 export BVPJacobianAlgorithm
+export maxsol, minsol
 
 end

lib/BoundaryValueDiffEqMIRK/src/interpolation.jl

@@ -158,6 +158,89 @@ function (s::EvalSol{C})(tval::Number) where {C <: MIRKCache}
     return z
 end
 
+# Interpolate intermidiate solution at multiple points
+function (s::EvalSol{C})(tvals::AbstractArray{<:Number}) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, stage, k_discrete, mesh_dt) = cache
+    # Quick handle for the case where tval is at the boundary
+    zvals = [zero(last(u)) for _ in tvals]
+    for (i, tval) in enumerate(tvals)
+        (tval == t[1]) && return first(u)
+        (tval == t[end]) && return last(u)
+        ii = interval(t, tval)
+        dt = mesh_dt[ii]
+        τ = (tval - t[ii]) / dt
+        w, _ = evalsol_interp_weights(τ, alg)
+        K = __needs_diffcache(alg.jac_alg) ? @view(k_discrete[ii].du[:, 1:stage]) :
+            @view(k_discrete[ii][:, 1:stage])
+        __maybe_matmul!(zvals[i], K, @view(w[1:stage]))
+        zvals[i] .= zvals[i] .* dt .+ u[ii]
+    end
+    return zvals
+end
+
+# Intermidiate derivative solution for evaluating boundry conditions
+function (s::EvalSol{C})(tval::Number, ::Type{Val{1}}) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, stage, k_discrete, mesh_dt) = cache
+    z′ = zeros(typeof(tval), 2)
+    ii = interval(t, tval)
+    dt = mesh_dt[ii]
+    τ = (tval - t[ii]) / dt
+    _, w′ = interp_weights(τ, alg)
+    __maybe_matmul!(z′, @view(k_discrete[ii].du[:, 1:stage]), @view(w′[1:stage]))
+    return z′
+end
+
+"""
+n root-finding problems to find the critical points with continuous derivative polynomials
+"""
+function __construct_then_solve_root_problem(
+        sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
+    (; alg) = sol.cache
+    n = first(size(sol))
+    nlprobs = Vector{NonlinearProblem}(undef, n)
+    for i in 1:n
+        f = @closure (t, p) -> sol(t, Val{1})[i]
+        nlprob = NonlinearProblem(f, sol.cache.prob.u0[i], tspan)
+        nlprobs[i] = nlprob
+    end
+    nlsols = Vector{SciMLBase.NonlinearSolution}(undef, length(nlprobs))
+    nlsolve_alg = __concrete_nonlinearsolve_algorithm(first(nlprobs), alg.nlsolve)
+    for (i, nlprob) in enumerate(nlprobs)
+        nlsols[i] = solve(nlprob, nlsolve_alg)
+    end
+    return nlsols
+end
+
+# It turns out the critical points can't cover all possible maximum/minimum values
+# especially when the solution are monotonic, we still need to compare the extremes with
+# value at critical points to find the maximum/minimum
+
+"""
+    maxsol(sol::EvalSol, tspan::Tuple)
+
+Find the maximum of the solution over the time span `tspan`.
+"""
+function maxsol(sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
+    nlsols = __construct_then_solve_root_problem(sol, tspan)
+    tvals = map(nlsol -> (SciMLBase.successful_retcode(nlsol); return nlsol.u), nlsols)
+    u = sol(tvals)
+    return max(maximum(sol), maximum(Iterators.flatten(u)))
+end
+
+"""
+    minsol(sol::EvalSol, tspan::Tuple)
+
+Find the minimum of the solution over the time span `tspan`.
+"""
+function minsol(sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
+    nlsols = __construct_then_solve_root_problem(sol, tspan)
+    tvals = map(nlsol -> (SciMLBase.successful_retcode(nlsol); return nlsol.u), nlsols)
+    u = sol(tvals)
+    return min(minimum(sol), minimum(Iterators.flatten(u)))
+end
+
 @inline function evalsol_interp_weights(τ::T, ::MIRK2) where {T}
     w = [0, τ * (1 - τ / 2), τ^2 / 2]
 

lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl

@@ -401,3 +401,22 @@ end
         @test SciMLBase.successful_retcode(sol)
     end
 end
+
+@testitem "Test maxsol and minsol" setup=[MIRKConvergenceTests] begin
+    tspan = (0.0, pi / 2)
+    function simplependulum!(du, u, p, t)
+        θ = u[1]
+        dθ = u[2]
+        du[1] = dθ
+        du[2] = -9.81 * sin(θ)
+    end
+    function bc!(residual, u, p, t)
+        residual[1] = maxsol(u, (0.0, pi / 2)) - 5.0496477654230745
+        residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
+    end
+    prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
+    for order in (4, 6)
+        sol = solve(prob, mirk_solver(Val(order)), dt = 0.05)
+        @test SciMLBase.successful_retcode(sol)
+    end
+end

src/BoundaryValueDiffEq.jl

@@ -32,4 +32,6 @@ export BVPM2, BVPSOL, COLNEW # From ODEInterface.jl
 
 export BVPJacobianAlgorithm
 
+export maxsol, minsol
+
 end