Merge pull request #332 from SciML/qqy/max_min_sol

Add maxsol and minsol

Author: ChrisRackauckas

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/extremum.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/extremum.md

@@ -0,0 +1,30 @@
+# Solve BVP with Extremum Boundary Conditions
+
+In many physical systems, boundary conditions are not always defined at fixed points such as intitial or terminal ends of the domain. Instead, we may encounter scenarios where constraints are imposed on the maximum or minimum values that the solution must attain somewhere within the domain. In such cases, we can use the `maxsol` and `minsol` functions provided by BoundaryValueDiffEq.jl to specify such extremum boundary conditions.
+
+Let's walk through this functionality with an intuitive example. We still revisit the simple pendulum example here, but this time, suppose we need to impose the maximum and minimum value to our boundary conditions, specified as:
+
+```math
+\max{u}=ub\\
+\min{u}=lb
+```
+
+where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must conform that the maximum value of the state should be `lb` while the minimum value of the state should be `ub`. To solve such problems, we can simply use the `maxsol` and `minsol` 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/default_nlsolve.jl

@@ -36,6 +36,18 @@ function __FastShortcutBVPCompatibleNonlinearPolyalg(
     return NonlinearSolvePolyAlgorithm(algs)
 end
 
+function __FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothing,
+        linsolve = nothing, autodiff = nothing) where {T}
+    if T <: Complex
+        algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),)
+    else
+        algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),
+            NewtonRaphson(; concrete_jac, linsolve, linesearch = BackTracking(), autodiff),
+            TrustRegion(; concrete_jac, linsolve, autodiff))
+    end
+    return NonlinearSolvePolyAlgorithm(algs)
+end
+
 @inline __concrete_nonlinearsolve_algorithm(prob, alg) = alg
 @inline function __concrete_nonlinearsolve_algorithm(prob, ::Nothing)
     if prob isa NonlinearLeastSquaresProblem

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

@@ -21,7 +21,8 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                DefectControl, GlobalErrorControl, SequentialErrorControl,
                                HybridErrorControl, HOErrorControl, __use_both_error_control,
                                __default_coloring_algorithm, DiffCacheNeeded,
-                               NoDiffCacheNeeded, __split_kwargs
+                               NoDiffCacheNeeded, __split_kwargs,
+                               __FastShortcutNonlinearPolyalg
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
@@ -160,5 +161,6 @@ end
 
 export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6, MIRK6I
 export BVPJacobianAlgorithm
+export maxsol, minsol
 
 end

lib/BoundaryValueDiffEqMIRK/src/interpolation.jl

@@ -158,6 +158,86 @@ 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
+
+"""
+Construct n root-finding problems and solve them 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{SciMLBase.NonlinearProblem}(undef, n)
+    nlsols = Vector{SciMLBase.NonlinearSolution}(undef, length(nlprobs))
+    nlsolve_alg = __FastShortcutNonlinearPolyalg(eltype(sol.cache))
+    for i in 1:n
+        f = @closure (t, p) -> sol(t, Val{1})[i]
+        nlprob = NonlinearProblem(f, sol.cache.prob.u0[i], tspan)
+        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