Skip to content

Commit 10ea8e6

Browse files
Merge pull request #769 from ErikQQY/qqy/revamp_bvp
Revamp BVP docs for new features
2 parents 3d98a4c + 4ee8b46 commit 10ea8e6

File tree

10 files changed

+168
-58
lines changed

10 files changed

+168
-58
lines changed

docs/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ AlgebraicMultigrid = "0.5, 0.6"
4444
BSON = "0.3"
4545
BVProblemLibrary = "0.1.2"
4646
BenchmarkTools = "1"
47-
BoundaryValueDiffEq = "5.1"
47+
BoundaryValueDiffEq = "5.14"
4848
CSV = "0.10"
4949
ComponentArrays = "0.15"
5050
DAEProblemLibrary = "0.1"

docs/make.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@ makedocs(
3939
"https://www.wolframalpha.com/input/?i=u%27%3D-sqrt%28u%29",
4040
"https://www.mathworks.com/help/simulink/gui/absolutetolerance.html",
4141
"https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html",
42-
"https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml",
42+
"https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml"
4343
],
4444
doctest = false, clean = true,
4545
warnonly = [:missing_docs],

docs/src/basics/plot.md

Lines changed: 11 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -75,52 +75,49 @@ the following conveniences are provided:
7575
- Everywhere in a tuple position where we only find an integer, this
7676
variable is plotted as a function of time. For example, the list above
7777
is equivalent to:
78-
78+
7979
```julia
8080
idxs = [1, (1, 3), (4, 5)]
8181
```
82-
82+
8383
and
84-
84+
8585
```julia
8686
idxs = [1, 3, 4]
8787
```
88-
88+
8989
is the most concise way to plot the variables 1, 3, and 4 as a function
9090
of time.
9191

9292
- It is possible to omit the list if only one plot is wanted: `(2,3)`
9393
and `4` are respectively equivalent to `[(2,3)]` and `[(0,4)]`.
94-
9594
- A tuple containing one or several lists will be expanded by
9695
associating corresponding elements of the lists with each other:
97-
96+
9897
```julia
9998
idxs = ([1, 2, 3], [4, 5, 6])
10099
```
101-
100+
102101
is equivalent to
103-
102+
104103
```julia
105104
idxs = [(1, 4), (2, 5), (3, 6)]
106105
```
107-
106+
108107
and
109-
108+
110109
```julia
111110
idxs = (1, [2, 3, 4])
112111
```
113-
112+
114113
is equivalent to
115-
114+
116115
```julia
117116
idxs = [(1, 2), (1, 3), (1, 4)]
118117
```
119-
120118
- Instead of using integers, one can use the symbols from a `ParameterizedFunction`.
121119
For example, `idxs=(:x,:y)` will replace the symbols with the integer values for
122120
components `:x` and `:y`.
123-
124121
- n-dimensional groupings are allowed. For example, `(1,2,3,4,5)` would be a
125122
5-dimensional plot between the associated variables.
126123

docs/src/examples/kepler_problem.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ sol = solve(prob, KahanLi6(), dt = 1 // 10);
2727
```
2828

2929
!!! note
30+
3031
Note that NonlinearSolve.jl is required to be imported for ManifoldProjection
3132

3233
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

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -190,7 +190,8 @@ end
190190

191191
function StochasticDiffEq.stepsize_controller!(integrator::StochasticDiffEq.SDEIntegrator,
192192
controller::CustomController, alg)
193-
integrator.q11 = DiffEqBase.value(FastPower.fastpower(integrator.EEst, controller.beta1))
193+
integrator.q11 = DiffEqBase.value(FastPower.fastpower(
194+
integrator.EEst, controller.beta1))
194195
integrator.q = DiffEqBase.value(integrator.q11 /
195196
FastPower.fastpower(integrator.qold, controller.beta2))
196197
integrator.q = DiffEqBase.value(max(inv(integrator.opts.qmax),

docs/src/features/io.md

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -81,11 +81,11 @@ work. If none of these are put into scope, the solution type
8181
will still load and hold all of the values (so `sol.u` and `sol.t`
8282
will work), but none of the interface will be available.
8383

84-
If you want a copy of the solution that contains no function information
84+
If you want a copy of the solution that contains no function information
8585
you can use the function `SciMLBase.strip_solution(sol)`.
86-
This will return a copy of the solution that doesn't have any functions,
86+
This will return a copy of the solution that doesn't have any functions,
8787
which you can serialize and deserialize without having any of the problems
88-
that typically come with serializing functions.
88+
that typically come with serializing functions.
8989

9090
## JLD
9191

docs/src/solvers/bvp_solve.md

Lines changed: 32 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,13 @@
1-
# BVP Solvers
1+
# [BVP Solvers](@id bvp_solve)
22

3-
`solve(prob::BVProblem,alg,dt=0.0;kwargs)`
3+
```julia
4+
solve(prob::BVProblem, alg, dt; kwargs)
5+
solve(prob::TwoPointBVProblem, alg, dt; kwargs)
6+
solve(prob::SecondOrderBVProblem, alg, dt; kwargs)
7+
solve(prob::SecondOrderTwoPointBVProblem, alg, dt; kwargs)
8+
```
49

5-
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.
10+
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.
611

712
## Recommended Methods
813

@@ -52,6 +57,30 @@ Similar to `MIRK` methods, fully implicit Runge-Kutta methods construct nonlinea
5257
- `RadauIIa5` - A 5th stage Radau collocation method.
5358
- `RadauIIa7` - A 7th stage Radau collocation method.
5459

60+
#### Gauss Legendre collocation methods
61+
62+
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.
63+
64+
- `Ascher1` - A 1st stage Gauss Legendre collocation method with Ascher's error control adaptivity.
65+
- `Ascher2` - A 2nd stage Gauss Legendre collocation method with Ascher's error control adaptivity.
66+
- `Ascher3` - A 3rd stage Gauss Legendre collocation method with Ascher's error control adaptivity.
67+
- `Ascher4` - A 4th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
68+
- `Ascher5` - A 5th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
69+
- `Ascher6` - A 6th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
70+
- `Ascher7` - A 7th stage Gauss Legendre collocation method with Ascher's error control adaptivity.
71+
72+
#### MIRKN(Monotonic Implicit Runge-Kutta-Nystöm) methods
73+
74+
- `MIRKN4` - A 4th order collocation method using an implicit Runge-Kutta-Nyström tableau without defect control adaptivity.
75+
- `MIRKN6` - A 6th order collocation method using an implicit Runge-Kutta-Nyström tableau without defect control adaptivity.
76+
77+
### SimpleBoundaryValueDiffEq.jl
78+
79+
- `SimpleMIRK4` - A simplified 4th order collocation method using an implicit Runge-Kutta tableau.
80+
- `SimpleMIRK5` - A simplified 5th order collocation method using an implicit Runge-Kutta tableau.
81+
- `SimpleMIRK6` - A simplified 6th order collocation method using an implicit Runge-Kutta tableau.
82+
- `SimpleShooting` - A simplified single Shooting method.
83+
5584
### ODEInterface.jl
5685

5786
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

Lines changed: 22 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -783,9 +783,9 @@ Examples:
783783

784784
```julia
785785
sol = solve(prob, Rosenbrock23()) # Standard, uses autodiff
786-
sol = solve(prob, Rosenbrock23(autodiff = AutoForwardDiff(chunksize = 10)) # Autodiff with chunksize of 10
787-
sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff()) # Numerical differentiation with central differencing
788-
sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff(fdtype = Val{:forward})) # Numerical differentiation with forward differencing
786+
sol = solve(prob, Rosenbrock23(autodiff = AutoForwardDiff(chunksize = 10))) # Autodiff with chunksize of 10
787+
sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff())) # Numerical differentiation with central differencing
788+
sol = solve(prob, Rosenbrock23(autodiff = AutoFiniteDiff(fdtype = Val{:forward}))) # Numerical differentiation with forward differencing
789789
```
790790

791791
#### Tableau Method
@@ -959,20 +959,25 @@ using IRKGaussLegendre
959959

960960
`IRKGL16(;kwargs...)` has the following arguments:
961961

962-
- second_order_ode (boolean):
963-
- =false (default): for a ODEProblem type
964-
- =true: for a second-order differential equation
965-
- simd (boolean):
966-
- =true: SIMD-vectorized implementation only available for Float32 or Float64 computations
967-
- =false (default): generic implementation that can use with arbitrary Julia-defined number systems
968-
- mstep: output saved at every 'mstep' steps. Default 1.
969-
- initial_extrapolation: initialization method for stages.
970-
- =false: simplest initialization
971-
- =true (default): extrapolation from the stage values of previous step
972-
- maxtrials: maximum number of attempts to accept adaptive step size
973-
- threading
974-
- =false (default): sequential execution of the numerical integration
975-
- =true: computations using threads (shared memory multi-threading) for stage-wise parallelization
962+
- second_order_ode (boolean):
963+
964+
+ =false (default): for a ODEProblem type
965+
+ =true: for a second-order differential equation
966+
967+
- simd (boolean):
968+
969+
+ =true: SIMD-vectorized implementation only available for Float32 or Float64 computations
970+
+ =false (default): generic implementation that can use with arbitrary Julia-defined number systems
971+
- mstep: output saved at every 'mstep' steps. Default 1.
972+
- initial_extrapolation: initialization method for stages.
973+
974+
+ =false: simplest initialization
975+
+ =true (default): extrapolation from the stage values of previous step
976+
- maxtrials: maximum number of attempts to accept adaptive step size
977+
- threading
978+
979+
+ =false (default): sequential execution of the numerical integration
980+
+ =true: computations using threads (shared memory multi-threading) for stage-wise parallelization
976981

977982
### SimpleDiffEq.jl
978983

docs/src/tutorials/bvp_example.md

Lines changed: 94 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -31,32 +31,26 @@ function simplependulum!(du, u, p, t)
3131
end
3232
```
3333

34-
### Boundary Condition
35-
3634
There are two problem types available:
3735

38-
- A problem type for general boundary conditions `BVProblem` (including conditions that may be anywhere/ everywhere on the integration interval).
39-
- A problem type for boundaries that are specified at the beginning and the end of the integration interval `TwoPointBVProblem`
40-
41-
#### `BVProblem`
36+
- A problem type for general boundary conditions `BVProblem` (including conditions that may be anywhere/ everywhere on the integration interval, aka multi-points BVP).
37+
- A problem type for boundaries that are specified at the beginning and the end of the integration interval `TwoPointBVProblem`(aka two-points BVP)
4238

4339
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.
4440

41+
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.
42+
4543
```@example bvp
4644
function bc1!(residual, u, p, t)
47-
residual[1] = u[end ÷ 2][1] + pi / 2 # the solution at the middle of the time span should be -pi/2
48-
residual[2] = u[end][1] - pi / 2 # the solution at the end of the time span should be pi/2
45+
residual[1] = u(pi / 4)[1] + pi / 2 # the solution at the middle of the time span should be -pi/2
46+
residual[2] = u(pi / 2)[1] - pi / 2 # the solution at the end of the time span should be pi/2
4947
end
5048
bvp1 = BVProblem(simplependulum!, bc1!, [pi / 2, pi / 2], tspan)
5149
sol1 = solve(bvp1, MIRK4(), dt = 0.05)
5250
plot(sol1)
5351
```
5452

55-
The third argument of `BVProblem` is the initial guess of the solution, which is constant in this example.
56-
57-
<!-- add examples of more general initial conditions -->
58-
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.
59-
If you can have a good initial guess, `Shooting` method works very well.
53+
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`.
6054

6155
```@example bvp
6256
using OrdinaryDiffEq
@@ -70,14 +64,12 @@ sol3 = solve(bvp3, Shooting(Vern7()))
7064
```
7165

7266
The initial guess can also be supplied via a function of `t` or a previous solution type, this is especially handy for parameter analysis.
73-
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.
67+
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).
7468

7569
```@example bvp
7670
plot(sol3)
7771
```
7872

79-
#### `TwoPointBVProblem`
80-
8173
`TwoPointBVProblem` is operationally the same as `BVProblem` but allows for the solver
8274
to specialize on the common form of being a two-point BVP, i.e. a BVP which only has
8375
boundary conditions at the start and the finish of the time interval.
@@ -99,5 +91,89 @@ plot(sol2)
9991
Note here that `bc2a!` is a boundary condition for the first time point, and `bc2b!` is a boundary condition
10092
for the final time point. `bcresid_prototype` is a prototype array which is passed in order to know the size of
10193
`resid_a` and `resid_b`. In this case, we have one residual term for the start and one for the final time point,
102-
and thus we have `bcresid_prototype = (zeros(1), zeros(1))`. If we wanted to only have boundary conditions at the
103-
final time, we could instead have done `bcresid_prototype = (zeros(0), zeros(2))`.
94+
and thus we have `bcresid_prototype = (zeros(1), zeros(1))`.
95+
96+
## Example 2: Directly Solving with Second Order BVP
97+
98+
Suppose we want to solve the second order BVP system which can be formulated as
99+
100+
```math
101+
\begin{cases}
102+
u_1''(x)= u_2(x),\\
103+
\epsilon u_2''(x)=-u_1(x)u_2'(x)- u_3(x)u_3'(x),\\
104+
\epsilon u_3''(x)=u_1'(x)u_3(x)- u_1(x) u_3 '(x)
105+
\end{cases}
106+
```
107+
108+
with initial conditions:
109+
110+
```math
111+
u_1(0) = u_1'(0)= u_1(1)=u_1'(1)=0,u_3(0)=
112+
-1, u_3(1)=1
113+
```
114+
115+
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.
116+
117+
```@example bvp
118+
function f!(ddu, du, u, p, t)
119+
ϵ = 0.1
120+
ddu[1] = u[2]
121+
ddu[2] = (-u[1] * du[2] - u[3] * du[3]) / ϵ
122+
ddu[3] = (du[1] * u[3] - u[1] * du[3]) / ϵ
123+
end
124+
function bc!(res, du, u, p, t)
125+
res[1] = u(0.0)[1]
126+
res[2] = u(1.0)[1]
127+
res[3] = u(0.0)[3] + 1
128+
res[4] = u(1.0)[3] - 1
129+
res[5] = du(0.0)[1]
130+
res[6] = du(1.0)[1]
131+
end
132+
u0 = [1.0, 1.0, 1.0]
133+
tspan = (0.0, 1.0)
134+
prob = SecondOrderBVProblem(f!, bc!, u0, tspan)
135+
sol = solve(prob, MIRKN4(; jac_alg = BVPJacobianAlgorithm(AutoForwardDiff())), dt = 0.01)
136+
```
137+
138+
## Example 3: Semi-Explicit Boundary Value Differential-Algebraic Equations
139+
140+
Consider a semi-explicit boundary value differential-algebraic equation formulated as
141+
142+
```math
143+
\begin{cases}
144+
x_1'=(\epsilon+x_2-p_2(t))y+p_1'(t) \\
145+
x_2'=p_2'(t) \\
146+
x_3'=y \\
147+
0=(x_1-p_1(t))(y-e^t)
148+
\end{cases}
149+
```
150+
151+
with boundary conditions
152+
153+
```math
154+
x_1(0)=0,x_3(0)=1,x_2(1)=\sin(1)
155+
```
156+
157+
We need to choose the Ascher methods for solving BVDAEs.
158+
159+
```@example bvp
160+
function f!(du, u, p, t)
161+
e = 2.7
162+
du[1] = (1 + u[2] - sin(t)) * u[4] + cos(t)
163+
du[2] = cos(t)
164+
du[3] = u[4]
165+
du[4] = (u[1] - sin(t)) * (u[4] - e^t)
166+
end
167+
function bc!(res, u, p, t)
168+
res[1] = u[1]
169+
res[2] = u[3] - 1
170+
res[3] = u[2] - sin(1.0)
171+
end
172+
u0 = [0.0, 0.0, 0.0, 0.0]
173+
tspan = (0.0, 1.0)
174+
fun = BVPFunction(f!, bc!, mass_matrix = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 0])
175+
prob = BVProblem(fun, u0, tspan)
176+
sol = solve(prob,
177+
Ascher4(zeta = [0.0, 0.0, 1.0], jac_alg = BVPJacobianAlgorithm(AutoForwardDiff())),
178+
dt = 0.01)
179+
```

docs/src/types/bvp_types.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
```@docs
44
SciMLBase.BVProblem
5+
SciMLBase.SecondOrderBVProblem
56
```
67

78
## Solution Type

0 commit comments

Comments
 (0)