Actually extremum boundary conditions

Author: ErikQQY

docs/pages.jl

@@ -2,8 +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/inequality.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)
+```

docs/src/tutorials/inequality.md

@@ -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/default_nlsolve.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

lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl

@@ -193,21 +193,18 @@ function (s::EvalSol{C})(tval::Number, ::Type{Val{1}}) where {C <: MIRKCache}
 end
 
 """
-n root-finding problems to find the critical points with continuous derivative polynomials
+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{NonlinearProblem}(undef, n)
+    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)
-        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