Revamp BVP docs for new features

docs/Project.toml

@@ -44,7 +44,7 @@ AlgebraicMultigrid = "0.5, 0.6"
 BSON = "0.3"
 BVProblemLibrary = "0.1.2"
 BenchmarkTools = "1"
-BoundaryValueDiffEq = "5.1"
+BoundaryValueDiffEq = "5.14"
 CSV = "0.10"
 ComponentArrays = "0.15"
 DAEProblemLibrary = "0.1"
@@ -70,7 +70,7 @@ OrdinaryDiffEqCore = "1"
 Plots = "1"
 RecursiveArrayTools = "3"
 SDEProblemLibrary = "0.1"
-SciMLBase = "2"
+SciMLBase = "2.71.3"
 SciMLOperators = "0.3"
 SparseConnectivityTracer = "0.6"
 StaticArrays = "1"

docs/make.jl

@@ -39,7 +39,7 @@ makedocs(
         "https://www.wolframalpha.com/input/?i=u%27%3D-sqrt%28u%29",
         "https://www.mathworks.com/help/simulink/gui/absolutetolerance.html",
         "https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html",
-        "https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml",
+        "https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml"
     ],
     doctest = false, clean = true,
     warnonly = [:missing_docs],

docs/src/basics/plot.md

@@ -75,52 +75,49 @@ the following conveniences are provided:
   - Everywhere in a tuple position where we only find an integer, this
     variable is plotted as a function of time.  For example, the list above
     is equivalent to:
-
+    
     ```julia
     idxs = [1, (1, 3), (4, 5)]
     ```
-
+    
     and
-
+    
     ```julia
     idxs = [1, 3, 4]
     ```
-
+    
     is the most concise way to plot the variables 1, 3, and 4 as a function
     of time.
 
   - It is possible to omit the list if only one plot is wanted: `(2,3)`
     and `4` are respectively equivalent to `[(2,3)]` and `[(0,4)]`.
-
   - A tuple containing one or several lists will be expanded by
     associating corresponding elements of the lists with each other:
-
+    
     ```julia
     idxs = ([1, 2, 3], [4, 5, 6])
     ```
-
+    
     is equivalent to
-
+    
     ```julia
     idxs = [(1, 4), (2, 5), (3, 6)]
     ```
-
+    
     and
-
+    
     ```julia
     idxs = (1, [2, 3, 4])
     ```
-
+    
     is equivalent to
-
+    
     ```julia
     idxs = [(1, 2), (1, 3), (1, 4)]
     ```
-
   - Instead of using integers, one can use the symbols from a `ParameterizedFunction`.
     For example, `idxs=(:x,:y)` will replace the symbols with the integer values for
     components `:x` and `:y`.
-
   - n-dimensional groupings are allowed. For example, `(1,2,3,4,5)` would be a
     5-dimensional plot between the associated variables.
 

docs/src/examples/kepler_problem.md

@@ -27,6 +27,7 @@ sol = solve(prob, KahanLi6(), dt = 1 // 10);
 ```
 
 !!! note
+
     Note that NonlinearSolve.jl is required to be imported for ManifoldProjection
 
 Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals.

docs/src/extras/timestepping.md

@@ -190,7 +190,8 @@ end
 
 function StochasticDiffEq.stepsize_controller!(integrator::StochasticDiffEq.SDEIntegrator,
         controller::CustomController, alg)
-    integrator.q11 = DiffEqBase.value(FastPower.fastpower(integrator.EEst, controller.beta1))
+    integrator.q11 = DiffEqBase.value(FastPower.fastpower(
+        integrator.EEst, controller.beta1))
     integrator.q = DiffEqBase.value(integrator.q11 /
                                     FastPower.fastpower(integrator.qold, controller.beta2))
     integrator.q = DiffEqBase.value(max(inv(integrator.opts.qmax),

docs/src/features/io.md

@@ -81,11 +81,11 @@ work. If none of these are put into scope, the solution type
 will still load and hold all of the values (so `sol.u` and `sol.t`
 will work), but none of the interface will be available.
 
-If you want a copy of the solution that contains no function information 
+If you want a copy of the solution that contains no function information
 you can use the function `SciMLBase.strip_solution(sol)`.
-This will return a copy of the solution that doesn't have any functions, 
+This will return a copy of the solution that doesn't have any functions,
 which you can serialize and deserialize without having any of the problems
-that typically come with serializing functions. 
+that typically come with serializing functions.
 
 ## JLD
 

docs/src/solvers/bvp_solve.md

@@ -1,8 +1,13 @@
-# BVP Solvers
+# [BVP Solvers](@id bvp_solve)
 
-`solve(prob::BVProblem,alg,dt=0.0;kwargs)`
+```julia
+solve(prob::BVProblem, alg, dt; kwargs)
+solve(prob::TwoPointBVProblem, alg, dt; kwargs)
+solve(prob::SecondOrderBVProblem, alg, dt; kwargs)
+solve(prob::SecondOrderTwoPointBVProblem, alg, dt; kwargs)
+```
 
-Solves the BVP defined by `prob` using the algorithm `alg`. All algorithms except `Shooting` methods should specify a `dt` which is the step size for the discretized mesh.
+Solves the BVP defined by `prob` using the algorithm `alg`. All algorithms except `Shooting` and `MultipleShooting` methods should specify a `dt` which is the step size for the discretized mesh.
 
 ## Recommended Methods
 
@@ -52,6 +57,30 @@ Similar to `MIRK` methods, fully implicit Runge-Kutta methods construct nonlinea
   - `RadauIIa5` - A 5th stage Radau collocation method.
   - `RadauIIa7` - A 7th stage Radau collocation method.
 
+#### Gauss Legendre collocation methods
+
+The `Ascher` collocation methods are similar with `MIRK` and `FIRK` methods but have extension for BVDAE prblem solving, the error control is based on instead of defect control adaptivity.
+
+  - `Ascher1` - A 1st stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+  - `Ascher2` - A 2nd stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+  - `Ascher3` - A 3rd stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+  - `Ascher4` - A 4th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+  - `Ascher5` - A 5th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+  - `Ascher6` - A 6th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+  - `Ascher7` - A 7th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
+
+#### MIRKN(Monotonic Implicit Runge-Kutta-Nystöm) methods
+
+  - `MIRKN4` - A 4th order collocation method using an implicit Runge-Kutta-Nyström tableau without defect control adaptivity.
+  - `MIRKN6` - A 6th order collocation method using an implicit Runge-Kutta-Nyström tableau without defect control adaptivity.
+
+### SimpleBoundaryValueDiffEq.jl
+
+  - `SimpleMIRK4` - A simplified 4th order collocation method using an implicit Runge-Kutta tableau.
+  - `SimpleMIRK5` - A simplified 5th order collocation method using an implicit Runge-Kutta tableau.
+  - `SimpleMIRK6` - A simplified 6th order collocation method using an implicit Runge-Kutta tableau.
+  - `SimpleShooting` - A simplified single Shooting method.
+
 ### ODEInterface.jl
 
 ODEInterface.jl can be used seamlessly with BoundaryValueDiffEq.jl, after we define our model using `BVProblem` or `TwoPointBVProblem`, we can directly call the solvers from ODEInterface.jl.

docs/src/solvers/ode_solve.md

@@ -783,9 +783,9 @@ Examples:
 
 ```julia
 sol = solve(prob, Rosenbrock23()) # Standard, uses autodiff
-sol = solve(prob, Rosenbrock23(autodiff = AutoForwardDiff(chunksize = 10)) # Autodiff with chunksize of 10
-sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff()) # Numerical differentiation with central differencing
-sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff(fdtype = Val{:forward})) # Numerical differentiation with forward differencing
+sol = solve(prob, Rosenbrock23(autodiff = AutoForwardDiff(chunksize = 10))) # Autodiff with chunksize of 10
+sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff())) # Numerical differentiation with central differencing
+sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff(fdtype = Val{:forward}))) # Numerical differentiation with forward differencing
 ```
 
 #### Tableau Method
@@ -959,20 +959,25 @@ using IRKGaussLegendre
 
 `IRKGL16(;kwargs...)` has the following arguments:
 
-- second_order_ode (boolean):
-  - =false (default): for a ODEProblem type 
-  - =true: for a second-order differential equation 
-- simd (boolean):
-  - =true: SIMD-vectorized implementation only available for Float32 or Float64 computations
-  - =false (default):  generic implementation that can use with arbitrary Julia-defined number systems 
-- mstep: output saved at every 'mstep' steps. Default 1.
-- initial_extrapolation: initialization method for stages.
-  - =false: simplest initialization
-  - =true (default): extrapolation from the stage values of previous step
-- maxtrials: maximum number of attempts to accept adaptive step size
-- threading
-  - =false (default): sequential execution of the numerical integration
-  - =true: computations using threads (shared memory multi-threading) for stage-wise parallelization
+  - second_order_ode (boolean):
+
+      + =false (default): for a ODEProblem type
+      + =true: for a second-order differential equation
+
+  - simd (boolean):
+
+      + =true: SIMD-vectorized implementation only available for Float32 or Float64 computations
+      + =false (default):  generic implementation that can use with arbitrary Julia-defined number systems
+  - mstep: output saved at every 'mstep' steps. Default 1.
+  - initial_extrapolation: initialization method for stages.
+
+      + =false: simplest initialization
+      + =true (default): extrapolation from the stage values of previous step
+  - maxtrials: maximum number of attempts to accept adaptive step size
+  - threading
+
+      + =false (default): sequential execution of the numerical integration
+      + =true: computations using threads (shared memory multi-threading) for stage-wise parallelization
 
 ### SimpleDiffEq.jl
 

docs/src/tutorials/bvp_example.md

@@ -31,32 +31,26 @@ function simplependulum!(du, u, p, t)
 end
 ```
 
-### Boundary Condition
-
 There are two problem types available:
 
-  - A problem type for general boundary conditions `BVProblem` (including conditions that may be anywhere/ everywhere on the integration interval).
-  - A problem type for boundaries that are specified at the beginning and the end of the integration interval `TwoPointBVProblem`
-
-#### `BVProblem`
+  - A problem type for general boundary conditions `BVProblem` (including conditions that may be anywhere/ everywhere on the integration interval, aka multi-points BVP).
+  - A problem type for boundaries that are specified at the beginning and the end of the integration interval `TwoPointBVProblem`(aka two-points BVP)
 
 The boundary conditions are specified by a function that calculates the residual in-place from the problem solution, such that the residual is $\vec{0}$ when the boundary condition is satisfied.
 
+There are collocation and shooting methods for addressing boundary value problems in DifferentialEquations.jl. We need to use appropriate [available BVP solvers](@ref bvp_solve) to solve `BVProblem`. In this example, we use `MIRK4` to solve the simple pendulum example.
+
 ```@example bvp
 function bc1!(residual, u, p, t)
-    residual[1] = u[end ÷ 2][1] + pi / 2 # the solution at the middle of the time span should be -pi/2
-    residual[2] = u[end][1] - pi / 2 # the solution at the end of the time span should be pi/2
+    residual[1] = u(pi / 4)[1] + pi / 2 # the solution at the middle of the time span should be -pi/2
+    residual[2] = u(pi / 2)[1] - pi / 2 # the solution at the end of the time span should be pi/2
 end
 bvp1 = BVProblem(simplependulum!, bc1!, [pi / 2, pi / 2], tspan)
 sol1 = solve(bvp1, MIRK4(), dt = 0.05)
 plot(sol1)
 ```
 
-The third argument of `BVProblem`  is the initial guess of the solution, which is constant in this example.
-
-<!-- add examples of more general initial conditions -->
-We need to use `MIRK4` or `Shooting` methods to solve `BVProblem`. `MIRK4` is a collocation method, whereas `Shooting` treats the problem as an IVP and varies the initial conditions until the boundary conditions are met.
-If you can have a good initial guess, `Shooting` method works very well.
+The third argument of `BVProblem` or `TwoPointBVProblem` is the initial guess of the solution, which can be specified as a `Vector`, a `Function` of `t` or solution object from previous solving, in this example the initial guess is set as a `Vector`.
 
 ```@example bvp
 using OrdinaryDiffEq
@@ -70,14 +64,12 @@ sol3 = solve(bvp3, Shooting(Vern7()))
 ```
 
 The initial guess can also be supplied via a function of `t` or a previous solution type, this is especially handy for parameter analysis.
-We changed `u` to `sol` to emphasize the fact that in this case, the boundary condition can be written on the solution object. Thus, all the features on the solution type such as interpolations are available when using the `Shooting` method. (i.e. you can have a boundary condition saying that the maximum over the interval is `1` using an optimization function on the continuous output). Note that user has to import the IVP solver before it can be used. Any common interface ODE solver is acceptable.
+We changed `u` to `sol` to emphasize the fact that in this case, the boundary condition can be written on the solution object. Thus, all the features on the solution type such as interpolations are available when using both collocation and shooting method. (i.e. you can have a boundary condition saying that the maximum over the interval is `1` using an optimization function on the continuous output).
 
 ```@example bvp
 plot(sol3)
 ```
 
-#### `TwoPointBVProblem`
-
 `TwoPointBVProblem` is operationally the same as `BVProblem` but allows for the solver
 to specialize on the common form of being a two-point BVP, i.e. a BVP which only has
 boundary conditions at the start and the finish of the time interval.
@@ -99,5 +91,89 @@ plot(sol2)
 Note here that `bc2a!` is a boundary condition for the first time point, and `bc2b!` is a boundary condition
 for the final time point. `bcresid_prototype` is a prototype array which is passed in order to know the size of
 `resid_a` and `resid_b`. In this case, we have one residual term for the start and one for the final time point,
-and thus we have `bcresid_prototype = (zeros(1), zeros(1))`. If we wanted to only have boundary conditions at the
-final time, we could instead have done `bcresid_prototype = (zeros(0), zeros(2))`.
+and thus we have `bcresid_prototype = (zeros(1), zeros(1))`.
+
+## Example 2: Directly Solving with Second Order BVP
+
+Suppose we want to solve the second order BVP system which can be formulated as
+
+```math
+\begin{cases}
+u_1''(x)= u_2(x),\\
+\epsilon u_2''(x)=-u_1(x)u_2'(x)- u_3(x)u_3'(x),\\
+\epsilon u_3''(x)=u_1'(x)u_3(x)- u_1(x) u_3 '(x)
+\end{cases}
+```
+
+with initial conditions:
+
+```math
+u_1(0) = u_1'(0)= u_1(1)=u_1'(1)=0,u_3(0)=
+-1, u_3(1)=1
+```
+
+The common way of solving the second order BVP is to define intermediate variables and transform the second order system into first order one, however, DifferentialEquations.jl allows the direct solving of second order BVP system to achieve more efficiency and higher continuity of the numerical solution.
+
+```@example bvp
+function f!(ddu, du, u, p, t)
+    ϵ = 0.1
+    ddu[1] = u[2]
+    ddu[2] = (-u[1] * du[2] - u[3] * du[3]) / ϵ
+    ddu[3] = (du[1] * u[3] - u[1] * du[3]) / ϵ
+end
+function bc!(res, du, u, p, t)
+    res[1] = u(0.0)[1]
+    res[2] = u(1.0)[1]
+    res[3] = u(0.0)[3] + 1
+    res[4] = u(1.0)[3] - 1
+    res[5] = du(0.0)[1]
+    res[6] = du(1.0)[1]
+end
+u0 = [1.0, 1.0, 1.0]
+tspan = (0.0, 1.0)
+prob = SecondOrderBVProblem(f!, bc!, u0, tspan)
+sol = solve(prob, MIRKN4(; jac_alg = BVPJacobianAlgorithm(AutoForwardDiff())), dt = 0.01)
+```
+
+## Example 3: Semi-Explicit Boundary Value Differential-Algebraic Equations
+
+Consider a semi-explicit boundary value differential-algebraic equation formulated as
+
+```math
+\begin{cases}
+x_1'=(\epsilon+x_2-p_2(t))y+p_1'(t) \\
+x_2'=p_2'(t) \\
+x_3'=y \\
+0=(x_1-p_1(t))(y-e^t)
+\end{cases}
+```
+
+with boundary conditions
+
+```math
+x_1(0)=0,x_3(0)=1,x_2(1)=\sin(1)
+```
+
+We need to choose the Ascher methods for solving BVDAEs.
+
+```@example bvp
+function f!(du, u, p, t)
+    e = 2.7
+    du[1] = (1 + u[2] - sin(t)) * u[4] + cos(t)
+    du[2] = cos(t)
+    du[3] = u[4]
+    du[4] = (u[1] - sin(t)) * (u[4] - e^t)
+end
+function bc!(res, u, p, t)
+    res[1] = u[1]
+    res[2] = u[3] - 1
+    res[3] = u[2] - sin(1.0)
+end
+u0 = [0.0, 0.0, 0.0, 0.0]
+tspan = (0.0, 1.0)
+fun = BVPFunction(f!, bc!, mass_matrix = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 0])
+prob = BVProblem(fun, u0, tspan)
+sol = solve(prob,
+    Ascher4(zeta = [0.0, 0.0, 1.0], jac_alg = BVPJacobianAlgorithm(AutoForwardDiff())),
+    dt = 0.01)
+```

docs/src/types/bvp_types.md

@@ -2,6 +2,7 @@
 
 ```@docs
 SciMLBase.BVProblem
+SciMLBase.SecondOrderBVProblem
 ```
 
 ## Solution Type