More docstrings

Author: ErikQQY

docs/src/assets/Project.toml

@@ -0,0 +1,40 @@
+[deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
+BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e"
+BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
+BoundaryValueDiffEqFIRK = "85d9eb09-370e-4000-bb32-543851f73618"
+BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
+BoundaryValueDiffEqMIRKN = "9255f1d6-53bf-473e-b6bd-23f1ff009da4"
+BoundaryValueDiffEqShooting = "ed55bfe0-3725-4db6-871e-a1dc9f42a757"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SimpleBoundaryValueDiffEq = "be0294bd-f90f-4760-ac4e-3421ce2b2da0"
+
+[compat]
+ADTypes = "1.9"
+BoundaryValueDiffEq = "5.13.0"
+BoundaryValueDiffEqAscher = "1.2.0"
+BoundaryValueDiffEqCore = "1.3.0"
+BoundaryValueDiffEqFIRK = "1.3.0"
+BoundaryValueDiffEqMIRK = "1.3.0"
+BoundaryValueDiffEqShooting = "1.3.0"
+DiffEqBase = "6.158.3"
+Documenter = "1"
+DocumenterCitations = "1"
+DocumenterInterLinks = "1.0.0"
+LineSearch = "0.1.4"
+LinearAlgebra = "1.10"
+LinearSolve = "2.36.2, 3"
+OrdinaryDiffEq = "6.90.1"
+Plots = "1"
+SciMLBase = "2.60.0"
+SimpleBoundaryValueDiffEq = "1.1.0"

docs/src/basics/error_control.md

@@ -4,13 +4,13 @@ Adaptivity helps ensure the quality of the our numerical solution, and when our
 
 When comes to solving ill-conditioned BVP, for example the singular pertubation problem where the small parameters become extremally small leading to the layers phonemona, the error control adaptivity becomes even more critical, because the minor pertubations can lead to large deviation in the solution. In such cases, adaptivity autimatically figure out where to use refined mesh and where to use coarse mesh to achieve the balance of computational efficiency and accuracy.
 
-BoundaryValuDiffEq.jl support error control adaptivity, and the adaptivity is default as on when using adaptive methods:
+BoundaryValuDiffEq.jl support error control adaptivity for collocation methods, and the adaptivity is default as defect control adaptivity when using adaptive collocation solvers:
 
 ```julia
 sol = solve(prob, MIRK4(), dt = 0.01, adaptive = true)
 ```
 
-Actually, BoundaryValueDiffEq.jl supports both defect and global error control adaptivity(while the defect control is the default), to specify different error control metods, we simply need to specify the `controller` keyword in `solve`:
+Actually, BoundaryValueDiffEq.jl supports both defect and global error control adaptivity(while the defect control is the default controller) [boisvert2013runge](@Citet), to specify different error control metods, we simply need to specify the `controller` keyword in `solve`:
 
 ```julia
 sol = solve(prob, MIRK4(), dt = 0.01, controller = GlobalErrorControl()) # Use global error control

docs/src/refs.bib

@@ -27,3 +27,14 @@ @book{ascher1995numerical
   URL = {https://epubs.siam.org/doi/abs/10.1137/1.9781611971231},
   eprint = {https://epubs.siam.org/doi/pdf/10.1137/1.9781611971231}
 }
+
+@article{boisvert2013runge,
+  title={A Runge-Kutta BVODE solver with global error and defect control},
+  author={Boisvert, Jason J and Muir, Paul H and Spiteri, Raymond J},
+  journal={ACM Transactions on Mathematical Software (TOMS)},
+  volume={39},
+  number={2},
+  pages={1--22},
+  year={2013},
+  publisher={ACM New York, NY, USA}
+}

lib/BoundaryValueDiffEqCore/src/calc_errors.jl

@@ -21,7 +21,7 @@ abstract type GlobalErrorControlMethod end
 Defect estimation method with defect as
 
 ```math
-DE = max\frac{|S'(x) - f(x,S(x))|}{1 + |f(x,S(x))}
+defect = \\max\\frac{S'(x) - f(x,S(x))}{1 + |f(x,S(x))|}
 ```
 
 Defect controller, with the maximum `defect_threshold` as 0.1, when the estimating defect is greater than the `defect_threshold`, the mesh will be refined.
@@ -87,7 +87,7 @@ Higher order global error estimation method
 Uses a solution from order+2 method on the original mesh and calculate the error with
 
 ```math
-GE = max\frac{|u_p - u_{p+2}|}{1 + |u_p|}
+error = \\max\\frac{u_p - u_{p+2}}{1 + |u_p|}
 ```
 """
 struct HOErrorControl <: GlobalErrorControlMethod end
@@ -99,7 +99,9 @@ Richardson extrapolation global error estimation method
 
 Use Richardson extrapolation to calculate the error on the doubled mesh with
 
-GE = \frac{2^p}{2^p-1} * max\frac{|u_h - u_{h/2}|}{1 + |u_h|}
+```math
+error = \\frac{2^p}{2^p-1} * \\max\\frac{u_h - u_{h/2}}{1 + |u_h|}
+```
 """
 struct REErrorControl <: GlobalErrorControlMethod end
 

lib/BoundaryValueDiffEqCore/src/utils.jl

@@ -144,8 +144,17 @@ function eval_bc_residual!(
     bcb!(resid[2], dub, ub, p)
 end
 
-# NOTE: We use `last` since the `first` might not conform to the same structure. For eg,
-#       in the case of residuals
+"""
+    __resize!(x, n, M)
+
+Resizes the input `x` to length `n` and returns the resized array. If `n` is less than the
+length of `x`, it truncates the array. If `n` is greater than the length of `x`, it appends
+zeros to the array.
+
+!!! note
+
+    We use `last` since the `first` might not conform to the same structure. For example, in the case of residuals
+"""
 function __resize!(x::AbstractVector{<:AbstractArray}, n, M)
     N = n - length(x)
     N == 0 && return x

lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl

@@ -24,7 +24,8 @@ end
 """
     mesh_selector!(cache::MIRKCache, controller::DefectControl)
     mesh_selector!(cache::MIRKCache, controller::GlobalErrorControl)
-    mesh_selector!(cache::MIRKCache, controller::Union{SequentialErrorControl, HybridErrorControl})
+    mesh_selector!(cache::MIRKCache, controller::SequentialErrorControl)
+    mesh_selector!(cache::MIRKCache, controller::HybridErrorControl)
 
 Generate new mesh based on the defect or the global error.
 """
@@ -35,7 +36,7 @@ Generate new mesh based on the defect or the global error.
     N = length(mesh)
 
     safety_factor = T(1.3)
-    ρ = T(1.0) # Set rho=1 means mesh distribution will take place everytime.
+    ρ = T(1.0)
     Nsub_star = 0
     Nsub_star_ub = 4 * (N - 1)
     Nsub_star_lb = N ÷ 2
@@ -95,7 +96,7 @@ end
 
     info = ReturnCode.Success
 
-    ŝ = [maximum(abs, d) for d in errors]  # Broadcasting breaks GPU Compilation
+    ŝ = [maximum(abs, d) for d in errors]
     ŝ .= (ŝ ./ abstol) .^ (T(1) / order)
     r₁ = maximum(ŝ)
     r₂ = sum(ŝ)
@@ -108,7 +109,8 @@ end
 
     if r₁ ≤ ρ * r₃
         Nsub_star = 2 * (N - 1)
-        if Nsub_star > cache.alg.max_num_subintervals # Need to determine the too large threshold
+        # Need to determine the too large threshold
+        if Nsub_star > cache.alg.max_num_subintervals
             info = ReturnCode.Failure
             meshₒ = mesh
             mesh_dt₀ = mesh_dt
@@ -148,7 +150,7 @@ end
 
     info = ReturnCode.Success
 
-    ŝ = [maximum(abs, d) for d in errors]  # Broadcasting breaks GPU Compilation
+    ŝ = [maximum(abs, d) for d in errors]
     ŝ .= (ŝ ./ abstol) .^ (T(1) / (order + 1))
     r₁ = maximum(ŝ)
     r₂ = sum(ŝ)
@@ -161,7 +163,8 @@ end
 
     if r₁ ≤ ρ * r₃
         Nsub_star = 2 * (N - 1)
-        if Nsub_star > cache.alg.max_num_subintervals # Need to determine the too large threshold
+        # Need to determine the too large threshold
+        if Nsub_star > cache.alg.max_num_subintervals
             info = ReturnCode.Failure
             meshₒ = mesh
             mesh_dt₀ = mesh_dt
@@ -216,7 +219,8 @@ end
 
     if r₁ ≤ ρ * r₃
         Nsub_star = 2 * (N - 1)
-        if Nsub_star > cache.alg.max_num_subintervals # Need to determine the too large threshold
+        # Need to determine the too large threshold
+        if Nsub_star > cache.alg.max_num_subintervals
             info = ReturnCode.Failure
             meshₒ = mesh
             mesh_dt₀ = mesh_dt
@@ -305,8 +309,7 @@ end
 """
     halve_sol(sol)
 
-The input sol has length of `n + 1`. Divide the original solution into two equal length
-solution.
+The input sol has length of `n + 1`. Divide the original mesh and u from original solution into `2n + 1` one.
 """
 function halve_sol(sol::AbstractVectorOfArray{T}, mesh) where {T}
     new_sol = copy(sol)
@@ -319,17 +322,17 @@ function halve_sol(sol::AbstractVectorOfArray{T}, mesh) where {T}
     @simd for i in (2n):-2:2
         new_sol[i] = (new_sol[i + 1] + new_sol[i - 1]) ./ T(2)
     end
-    mes = deepcopy(mesh)
-    resize!(mes, 2 * n + 1)
-    mes[1] = mesh[1]
-    mes[end] = mesh[end]
+    new_mesh = deepcopy(mesh)
+    resize!(new_mesh, 2 * n + 1)
+    new_mesh[1] = mesh[1]
+    new_mesh[end] = mesh[end]
     for i in (2n - 1):-2:1
-        mes[i] = mes[(i + 1) ÷ 2]
+        new_mesh[i] = new_mesh[(i + 1) ÷ 2]
     end
     for i in (2n):-2:2
-        mes[i] = (mes[i + 1] + mes[i - 1]) / 2
+        new_mesh[i] = (new_mesh[i + 1] + new_mesh[i - 1]) / 2
     end
-    return DiffEqArray(new_sol.u, mes)
+    return DiffEqArray(new_sol.u, new_mesh)
 end
 
 """
@@ -338,16 +341,21 @@ end
     error_estimate!(cache::MIRKCache, controller::SequentialErrorControl)
     error_estimate!(cache::MIRKCache, controller::HybridErrorControl)
 
+## Defect Control
+
 error_estimate for the defect uses the discrete solution approximation Y, plus stages of
 the RK method in 'k_discrete', plus some new stages in 'k_interp' to construct
 an interpolant.
 
+## Global Error Control
 error_estimate for the global error use the higher order or doubled mesh to estiamte the
 global error according to err = max(abs(Y_high - Y_low)) / (1 + abs(Y_low))
 
+## Sequential Error Control
 error_estimate for the sequential error first uses the defect controller, if the defect is
 satisfying, then use the global error controller.
 
+## Hybrid Error Control
 error_estimate for the hybrid error control uses the linear combination of defect and global
 error to estimate the error norm.
 """