Audit docs for DifferentialEquations.jl v8 solver provenance

docs/src/examples/classical_physics.md

@@ -70,7 +70,9 @@ we can use `SecondOrderODEProblem` as follows.
 
 ```@example physics
 # Simple Harmonic Oscillator Problem
-import OrdinaryDiffEq as ODE, Plots
+import OrdinaryDiffEq as ODE
+import OrdinaryDiffEqRKN as ODERKN # DPRKN6
+import Plots
 
 #Parameters
 ω = 1
@@ -90,7 +92,7 @@ end
 
 #Pass to solvers
 prob = ODE.SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
-sol = ODE.solve(prob, ODE.DPRKN6())
+sol = ODE.solve(prob, ODERKN.DPRKN6())
 
 #Plot
 Plots.plot(sol, idxs = [2, 1], linewidth = 2, title = "Simple Harmonic Oscillator",
@@ -382,6 +384,7 @@ Plots.plot(sol.t, energy .- energy[1], title = "Change in Energy over Time",
 To prevent energy drift, we can instead use a symplectic integrator. We can directly define and solve the `SecondOrderODEProblem`:
 
 ```@example physics
+import OrdinaryDiffEqSymplecticRK as ODESymp # KahanLi8
 function HH_acceleration!(dv, v, u, p, t)
     x, y = u
     dx, dy = dv
@@ -391,7 +394,7 @@ end
 initial_positions = [0.0, 0.1]
 initial_velocities = [0.5, 0.0]
 prob = ODE.SecondOrderODEProblem(HH_acceleration!, initial_velocities, initial_positions, tspan)
-sol2 = ODE.solve(prob, ODE.KahanLi8(), dt = 1 / 10);
+sol2 = ODE.solve(prob, ODESymp.KahanLi8(), dt = 1 / 10);
 ```
 
 Notice that we get the same results:
@@ -424,7 +427,7 @@ Plots.plot(sol2.t, energy .- energy[1], title = "Change in Energy over Time",
 And let's try to use a Runge-Kutta-Nyström solver to solve this. Note that Runge-Kutta-Nyström isn't symplectic.
 
 ```@example physics
-sol3 = ODE.solve(prob, ODE.DPRKN6());
+sol3 = ODE.solve(prob, ODERKN.DPRKN6());
 energy = map(x -> E(x[3], x[4], x[1], x[2]), sol3.u)
 @show ΔE = energy[1] - energy[end]
 Plots.gr()

docs/src/examples/kepler_problem.md

@@ -19,7 +19,11 @@ Also, we know that
 ```
 
 ```@example kepler
-import OrdinaryDiffEq as ODE, ForwardDiff, Plots
+import OrdinaryDiffEq as ODE
+import OrdinaryDiffEqSymplecticRK as ODESymp # KahanLi6
+import OrdinaryDiffEqRKN as ODERKN           # DPRKN6, ERKN4
+import OrdinaryDiffEqLowOrderRK as ODELow    # RK4
+import ForwardDiff, Plots
 import LinearAlgebra: norm
 H(q, p) = norm(p)^2 / 2 - inv(norm(q))
 L(q, p) = q[1] * p[2] - p[1] * q[2]
@@ -33,7 +37,7 @@ initial_cond = (initial_position, initial_velocity)
 initial_first_integrals = (H(initial_cond...), L(initial_cond...))
 tspan = (0, 20.0)
 prob = ODE.DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)
-sol = ODE.solve(prob, ODE.KahanLi6(), dt = 1 // 10);
+sol = ODE.solve(prob, ODESymp.KahanLi6(), dt = 1 // 10);
 ```
 
 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.
@@ -59,7 +63,7 @@ analysis_plot(sol, H, L)
 Let's try to use a Runge-Kutta-Nyström solver to solve this problem and check the first integrals' variation.
 
 ```@example kepler
-sol2 = ODE.solve(prob, ODE.DPRKN6())  # dt is not necessary, because unlike symplectic
+sol2 = ODE.solve(prob, ODERKN.DPRKN6())  # dt is not necessary, because unlike symplectic
 # integrators DPRKN6 is adaptive
 @show sol2.u |> length
 analysis_plot(sol2, H, L)
@@ -68,7 +72,7 @@ analysis_plot(sol2, H, L)
 Let's then try to solve the same problem by the `ERKN4` solver, which is specialized for sinusoid-like periodic function
 
 ```@example kepler
-sol3 = ODE.solve(prob, ODE.ERKN4()) # dt is not necessary, because unlike symplectic
+sol3 = ODE.solve(prob, ODERKN.ERKN4()) # dt is not necessary, because unlike symplectic
 # integrators ERKN4 is adaptive
 @show sol3.u |> length
 analysis_plot(sol3, H, L)
@@ -127,7 +131,7 @@ function hamiltonian(du, u, params, t)
 end
 
 prob2 = ODE.ODEProblem(hamiltonian, [initial_position; initial_velocity], tspan)
-sol_ = ODE.solve(prob2, ODE.RK4(), dt = 1 // 5, adaptive = false)
+sol_ = ODE.solve(prob2, ODELow.RK4(), dt = 1 // 5, adaptive = false)
 analysis_plot2(sol_, H, L)
 ```
 
@@ -140,7 +144,7 @@ function first_integrals_manifold(residual, u, p, t)
 end
 
 cb = CB.ManifoldProjection(first_integrals_manifold, autodiff = NLS.AutoForwardDiff())
-sol5 = ODE.solve(prob2, ODE.RK4(), dt = 1 // 5, adaptive = false, callback = cb)
+sol5 = ODE.solve(prob2, ODELow.RK4(), dt = 1 // 5, adaptive = false, callback = cb)
 analysis_plot2(sol5, H, L)
 ```
 
@@ -152,7 +156,7 @@ function energy_manifold(residual, u, p, t)
     residual[3:4] .= 0
 end
 energy_cb = CB.ManifoldProjection(energy_manifold, autodiff = NLS.AutoForwardDiff())
-sol6 = ODE.solve(prob2, ODE.RK4(), dt = 1 // 5, adaptive = false, callback = energy_cb)
+sol6 = ODE.solve(prob2, ODELow.RK4(), dt = 1 // 5, adaptive = false, callback = energy_cb)
 analysis_plot2(sol6, H, L)
 ```
 
@@ -164,7 +168,7 @@ function angular_manifold(residual, u, p, t)
     residual[3:4] .= 0
 end
 angular_cb = CB.ManifoldProjection(angular_manifold, autodiff = NLS.AutoForwardDiff())
-sol7 = ODE.solve(prob2, ODE.RK4(), dt = 1 // 5, adaptive = false, callback = angular_cb)
+sol7 = ODE.solve(prob2, ODELow.RK4(), dt = 1 // 5, adaptive = false, callback = angular_cb)
 analysis_plot2(sol7, H, L)
 ```
 

docs/src/extras/timestepping.md

@@ -232,6 +232,6 @@ end
 and used via
 
 ```julia
-# Assuming: import DifferentialEquations as DE
-sol = DE.solve(prob, DE.EM(), dt = dt, adaptive = true, controller = CustomController())
+import StochasticDiffEq as SDE # EM
+sol = SDE.solve(prob, SDE.EM(), dt = dt, adaptive = true, controller = CustomController())
 ```

docs/src/features/callback_functions.md

@@ -298,6 +298,7 @@ that we implemented. Now
 
 ```@example callback3
 import DifferentialEquations as DE
+import OrdinaryDiffEqLowOrderRK as ODELow # BS3
 cb = AutoAbstol(true; init_curmax = 1e-6)
 ```
 
@@ -313,7 +314,7 @@ end
 u0 = 10.0
 const V = 1
 prob = DE.ODEProblem(g, u0, (0.0, 10.0))
-integrator = DE.init(prob, DE.BS3(), callback = cb)
+integrator = DE.init(prob, ODELow.BS3(), callback = cb)
 at1 = integrator.opts.abstol
 DE.step!(integrator)
 at2 = integrator.opts.abstol

docs/src/features/dae_initialization.md

@@ -68,7 +68,9 @@ DiffEqBase.ShampineCollocationInit
 ### Example 1: Simple Pendulum DAE
 
 ```julia
-using DifferentialEquations
+using OrdinaryDiffEqBDF       # DFBDF
+using SciMLBase               # DAEProblem, solve
+using DiffEqBase              # BrownFullBasicInit, CheckInit, NoInit
 
 function pendulum!(res, du, u, p, t)
     x, y, T = u
@@ -114,7 +116,9 @@ sol = solve(prob2, DFBDF(), initializealg = CheckInit())
 When using ModelingToolkit, initialization information is often included automatically:
 
 ```julia
-using ModelingToolkit, DifferentialEquations
+using ModelingToolkit
+using OrdinaryDiffEqBDF       # DFBDF
+using SciMLBase               # DAEProblem, solve
 
 @variables t x(t) y(t) T(t)
 @parameters g L
@@ -143,9 +147,9 @@ Both OrdinaryDiffEq and Sundials support the same initialization algorithms thro
 ### OrdinaryDiffEq and Sundials
 
 ```julia
-using OrdinaryDiffEq
+using OrdinaryDiffEqBDF   # DFBDF
 # or
-using Sundials
+using Sundials            # IDA
 
 # Use Brown's algorithm to fix inconsistent conditions
 sol = solve(prob, DFBDF(), initializealg = BrownFullBasicInit())  # OrdinaryDiffEq

docs/src/features/ensemble.md

@@ -390,8 +390,9 @@ Now we build the SDE with these functions:
 
 ```@example ensemble2
 import DifferentialEquations as DE
+import StochasticDiffEq as SDE # SDEProblem, SRIW1
 p = [1.5, 1.0, 0.1, 0.1]
-prob = DE.SDEProblem(f, g, [1.0, 1.0], (0.0, 10.0), p)
+prob = SDE.SDEProblem(f, g, [1.0, 1.0], (0.0, 10.0), p)
 ```
 
 This is the base problem for our study. What would like to do with this experiment
@@ -416,8 +417,8 @@ end
 Now we solve the problem 10 times and plot all of the trajectories in phase space:
 
 ```@example ensemble2
-ensemble_prob = DE.EnsembleProblem(prob, prob_func = prob_func)
-sim = DE.solve(ensemble_prob, DE.SRIW1(), trajectories = 10);
+ensemble_prob = SDE.EnsembleProblem(prob, prob_func = prob_func)
+sim = SDE.solve(ensemble_prob, SDE.SRIW1(), trajectories = 10);
 nothing # hide
 ```
 

docs/src/features/io.md

@@ -12,11 +12,13 @@ type as the data source to convert to other tabular data formats. For example,
 let's solve a 4x2 system of ODEs and get the DataFrame:
 
 ```@example IO
-import OrdinaryDiffEq as ODE, DataFrames
+import OrdinaryDiffEq as ODE
+import OrdinaryDiffEqLowOrderRK as ODELow # Euler
+import DataFrames
 f_2dlinear = (du, u, p, t) -> du .= 1.01u;
 tspan = (0.0, 1.0)
 prob = ODE.ODEProblem(f_2dlinear, rand(2, 2), tspan);
-sol = ODE.solve(prob, ODE.Euler(); dt = 1 // 2^(4));
+sol = ODE.solve(prob, ODELow.Euler(); dt = 1 // 2^(4));
 df = DataFrames.DataFrame(sol)
 ```
 
@@ -42,7 +44,7 @@ JLD2.jl and BSON.jl will work with the full solution type if you bring the requi
 back into scope before loading. For example, if we save the solution:
 
 ```@example IO
-sol = ODE.solve(prob, ODE.Euler(); dt = 1 // 2^(4))
+sol = ODE.solve(prob, ODELow.Euler(); dt = 1 // 2^(4))
 import JLD2
 JLD2.@save "out.jld2" sol
 ```
@@ -60,7 +62,7 @@ JLD2.@load "out.jld2" sol
 The example with BSON.jl is:
 
 ```@example IO
-sol = ODE.solve(prob, ODE.Euler(); dt = 1 // 2^(4))
+sol = ODE.solve(prob, ODELow.Euler(); dt = 1 // 2^(4))
 import BSON
 BSON.bson("test.bson", Dict(:sol => sol))
 ```

docs/src/features/linear_nonlinear.md

@@ -24,9 +24,13 @@ For linear solvers, DifferentialEquations.jl uses
 [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl). Any
 [LinearSolve.jl algorithm](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/)
 can be used as the linear solver simply by passing the algorithm choice to
-linsolve. For example, the following tells `TRBDF2` to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl)
+linsolve. For example, the following tells `TRBDF2` (from
+`OrdinaryDiffEqSDIRK`) to use [KLU.jl](https://github.com/JuliaSparse/KLU.jl)
+(via `KLUFactorization` from `LinearSolve`):
 
 ```julia
+using OrdinaryDiffEqSDIRK: TRBDF2
+using LinearSolve: KLUFactorization
 TRBDF2(linsolve = KLUFactorization())
 ```
 
@@ -41,7 +45,8 @@ can be used as a left or right preconditioner. Preconditioners are configured on
 object itself, through LinearSolve's `Pl` / `Pr` interface. For example:
 
 ```julia
-using LinearSolve
+using LinearSolve              # KrylovJL_GMRES
+using OrdinaryDiffEqSDIRK: TRBDF2
 alg = TRBDF2(linsolve = KrylovJL_GMRES(precs = mypreconditioner))
 ```
 

docs/src/features/low_dep.md

@@ -102,12 +102,15 @@ so there is no issue here!
 
 For the add-on packages, you will normally need SciMLBase, the solver package
 you choose, and the add-on package. So for example, for predefined callbacks you
-would likely want SciMLBase+OrdinaryDiffEq+DiffEqCallbacks. If you aren't sure
-which package a specific command is from, then use `@which`. For example, from
-the callback docs we have:
+would likely want SciMLBase+OrdinaryDiffEq+DiffEqCallbacks. For example, the
+callback `ProbIntsUncertainty` lives in `DiffEqCallbacks` and the explicit Euler
+solver lives in `OrdinaryDiffEqLowOrderRK`, so this fully-explicit form replaces
+the old `using DifferentialEquations` shape:
 
 ```@example low_dep_1
-import DifferentialEquations as DE
+import OrdinaryDiffEq as ODE
+import OrdinaryDiffEqLowOrderRK as ODELow # Euler
+import DiffEqCallbacks as CB              # ProbIntsUncertainty
 function fitz(du, u, p, t)
     V, R = u
     a, b, c = p
@@ -117,41 +120,21 @@ end
 u0 = [-1.0; 1.0]
 tspan = (0.0, 20.0)
 p = (0.2, 0.2, 3.0)
-prob = DE.ODEProblem(fitz, u0, tspan, p)
-cb = DE.ProbIntsUncertainty(0.2, 1)
-ensemble_prob = DE.EnsembleProblem(prob)
-sim = DE.solve(ensemble_prob, DE.Euler(), trajectories = 100, callback = cb, dt = 1 / 10)
+prob = ODE.ODEProblem(fitz, u0, tspan, p)
+cb = CB.ProbIntsUncertainty(0.2, 1)
+ensemble_prob = ODE.EnsembleProblem(prob)
+sim = ODE.solve(ensemble_prob, ODELow.Euler(), trajectories = 100, callback = cb, dt = 1 / 10)
 ```
 
-If we wanted to know where `ProbIntsUncertainty(0.2,1)` came from, we can do:
+If you encounter an unfamiliar name and want to discover which package owns it,
+use `@which`:
 
 ```@example low_dep_1
 import InteractiveUtils # hide
-InteractiveUtils.@which DE.ProbIntsUncertainty(0.2, 1)
-```
-
-This says it's in the DiffEqCallbacks.jl package. Thus in this case, we could have
-done
-
-```@example low_dep_2
-import OrdinaryDiffEq as ODE, DiffEqCallbacks as CB
-function fitz(du, u, p, t)
-    V, R = u
-    a, b, c = p
-    du[1] = c * (V - V^3 / 3 + R)
-    du[2] = -(1 / c) * (V - a - b * R)
-end
-u0 = [-1.0; 1.0]
-tspan = (0.0, 20.0)
-p = (0.2, 0.2, 3.0)
-prob = ODE.ODEProblem(fitz, u0, tspan, p)
-cb = CB.ProbIntsUncertainty(0.2, 1)
-ensemble_prob = ODE.EnsembleProblem(prob)
-sim = ODE.solve(ensemble_prob, ODE.Euler(), trajectories = 100, callback = cb, dt = 1 / 10)
+InteractiveUtils.@which CB.ProbIntsUncertainty(0.2, 1)
 ```
 
-instead of the full `using DifferentialEquations`. Note that due to the way
-Julia dependencies work, any internal function in the package will work. The only
-dependencies you need to explicitly `using` are the functions you are specifically
-calling. Thus, this method can be used to determine all of the DiffEq packages
-you are using.
+Note that due to the way Julia dependencies work, any internal function in the
+package will work. The only dependencies you need to explicitly `using` are the
+functions you are specifically calling. Thus, this method can be used to
+determine all of the DiffEq packages you are using.

docs/src/features/verbosity.md

@@ -33,6 +33,7 @@ When using auto-switching algorithms, see when and why switches occur:
 
 ```julia
 using ODEProblemLibrary: prob_ode_vanderpol_stiff
+using OrdinaryDiffEqRosenbrock: Rodas5
 
 verbose = DEVerbosity(alg_switch = SciMLLogging.InfoLevel())
 sol = solve(prob_ode_vanderpol_stiff, AutoTsit5(Rodas5()), verbose = verbose)
@@ -110,6 +111,7 @@ The simplest approach is to pass a SciMLLogging preset, which is automatically c
 
 ```julia
 using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg
+using OrdinaryDiffEqSDIRK: ImplicitEuler
 
 # Enable detailed nonlinear solver diagnostics
 verbose = DEVerbosity(nonlinear_verbosity = SciMLLogging.Detailed())