diff --git a/docs/Project.toml b/docs/Project.toml index e47d5e719..c13ec2261 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index 5603a57cb..17f0b7dbb 100644 --- a/docs/make.jl +++ b/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], diff --git a/docs/src/basics/plot.md b/docs/src/basics/plot.md index 2e27fea2e..8c8cff542 100644 --- a/docs/src/basics/plot.md +++ b/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. diff --git a/docs/src/examples/kepler_problem.md b/docs/src/examples/kepler_problem.md index 019d8167c..c8a22c066 100644 --- a/docs/src/examples/kepler_problem.md +++ b/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. diff --git a/docs/src/extras/timestepping.md b/docs/src/extras/timestepping.md index be3de8870..9dfe3bfba 100644 --- a/docs/src/extras/timestepping.md +++ b/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), diff --git a/docs/src/features/io.md b/docs/src/features/io.md index 93bf2c4ce..bbc82316e 100644 --- a/docs/src/features/io.md +++ b/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 diff --git a/docs/src/solvers/bvp_solve.md b/docs/src/solvers/bvp_solve.md index 826240a5f..d6e294512 100644 --- a/docs/src/solvers/bvp_solve.md +++ b/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. diff --git a/docs/src/solvers/ode_solve.md b/docs/src/solvers/ode_solve.md index 1dbbf401f..6c5c155d8 100644 --- a/docs/src/solvers/ode_solve.md +++ b/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 diff --git a/docs/src/tutorials/bvp_example.md b/docs/src/tutorials/bvp_example.md index b55f1ab42..8d74e0240 100644 --- a/docs/src/tutorials/bvp_example.md +++ b/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. - - -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) +``` diff --git a/docs/src/types/bvp_types.md b/docs/src/types/bvp_types.md index c809a9074..cf51e710f 100644 --- a/docs/src/types/bvp_types.md +++ b/docs/src/types/bvp_types.md @@ -2,6 +2,7 @@ ```@docs SciMLBase.BVProblem +SciMLBase.SecondOrderBVProblem ``` ## Solution Type