Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master'
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Feb 3, 2018
2 parents 6e974a5 + 344beb7 commit fc7118b
Show file tree
Hide file tree
Showing 5 changed files with 65 additions and 30 deletions.
6 changes: 3 additions & 3 deletions docs/src/basics/common_solver_opts.md
Original file line number Diff line number Diff line change
Expand Up @@ -176,9 +176,9 @@ explanations of the timestepping algorithms, see the
* `callback`: Specifies a callback. Defaults to a callback function which
performs the saving routine. For more information, see the
[Event Handling and Callback Functions manual page](../../features/callback_functions.html).
* `isoutofdomain`: Specifies a function `isoutofdomain(t,u)` where, when it
returns false, it will reject the timestep. Defaults to always false.
* `unstable_check`: Specifies a function `unstable_check(dt,t,u)` where, when
* `isoutofdomain`: Specifies a function `isoutofdomain(u,p,t)` where, when it
returns true, it will reject the timestep. Disabled by default.
* `unstable_check`: Specifies a function `unstable_check(dt,u,p,t)` where, when
it returns true, it will cause the solver to exit and throw a warning. Defaults
to `any(isnan,u)`, i.e. checking if any value is a NaN.
* `verbose`: Toggles whether warnings are thrown when the solver exits early.
Expand Down
4 changes: 2 additions & 2 deletions docs/src/features/callback_library.md
Original file line number Diff line number Diff line change
Expand Up @@ -280,9 +280,9 @@ for the saved values, and then call the solver with the callback.

```julia
using DiffEqCallbacks, OrdinaryDiffEq
prob = ODEProblem((du,u,p,t)->du.=u,rand(4,4),(0.0,1.0))
prob = ODEProblem((du,u,t)->du.=u,rand(4,4),(0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,p,t,integrator)->(trace(u),norm(u)), saved_values)
cb = SavingCallback((u,t,integrator)->(trace(u),norm(u)), saved_values)
sol = solve(prob, Tsit5(), callback=cb)

print(saved_values.saveval)
Expand Down
2 changes: 1 addition & 1 deletion docs/src/features/performance_overloads.md
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ is necessary. For the DAE
G(du,u,p,t) = res
```

The Jacobian should be given in the form `dG/du + gamma*dG/d(du)` where `gamma`
The Jacobian should be given in the form `dG/d(du) + gamma*dG/du` where `gamma`
is given by the solver. This means that the signature is:

```julia
Expand Down
56 changes: 34 additions & 22 deletions docs/src/tutorials/bvp_example.md
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ using BoundaryValueDiffEq
const g = 9.81
L = 1.0
tspan = (0.0,pi/2)
function simplependulum(du,u,p,t)
function simplependulum!(du,u,p,t)
θ = u[1]
= u[2]
du[1] =
Expand All @@ -31,51 +31,63 @@ end

### Boundary Condition

And here is where the `Boundary` comes in. We need to write a function that calculate the residual in-place from the problem solution, such that the residual is $\vec{0}$ when the boundary condition is satisfied.
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`

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.

```julia
function bc1(residual, u, p)
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
end
bvp1 = BVProblem(simplependulum, bc1, [pi/2,pi/2], tspan)
bvp1 = BVProblem(simplependulum!, bc1!, [pi/2,pi/2], tspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt=0.05)
plot(sol1)
```

![BVP Example Plot1](../assets/bvp_example_plot1.png)

We need to use `GeneralMIRK4` or `Shooting` methods to solve `BVProblem`. We have boundary conditions at the beginning and the ending of the time span in common cases. We can use the `TwoPointBVProblem` problem type for such cases.

```julia
function bc2(residual, ua, ub, p) # ua is the beginning of the time span, and ub is the ending
residual[1] = ua[1] + pi/2 # the solution at the beginning of the time span should be -pi/2
residual[2] = ub[1] - pi/2 # the solution at the end of the time span should be pi/2
end
bvp2 = TwoPointBVProblem(simplependulum, bc2, [pi/2,pi/2], tspan)
sol2 = solve(bvp2, MIRK4(), dt=0.05) # we need to use the MIRK4 solver for TwoPointBVProblem
plot(sol2)
```

![BVP Example Plot2](../assets/bvp_example_plot2.png)

We have used the mono-implicit Runge–Kutta (MIRK) method to solve the BVP, but we can always use reduce a BVP to an IVP and a root-finding problem, which is the `Shooting` method. If you can have a good initial guess, shooting method works very well.
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 `GeneralMIRK4` or `Shooting` methods to solve `BVProblem`. `GeneralMIRK4` 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.

```julia
using OrdinaryDiffEq
u₀_2 = [-1.6, -1.7] # the initial guess
function bc3(residual, sol, p)
function bc3!(residual, sol, p)
residual[1] = sol(pi/4)[1] + pi/2 # use the interpolation here, since indexing will be wrong for adaptive methods
residual[2] = sol(pi/2)[1] - pi/2
end
bvp3 = BVProblem(simplependulum, bc3, u₀_2, tspan)
bvp3 = BVProblem(simplependulum!, bc3!, u₀_2, tspan)
sol3 = solve(bvp3, Shooting(Vern7()))
```

The initial guess can also be supplied via a function of `t` or a previous solution type, this is espacially 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 of 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.

```julia
plot(sol3)
```

![BVP Example Plot3](../assets/bvp_example_plot3.png)

#### `TwoPointBVProblem`

Defining a similar problem as `TwoPointBVProblem` is shown in the following example. At the moment `MIRK4` is the only solver for `TwoPointBVProblem`s.

```julia
function bc2!(residual, u, p, t) # u[1] is the beginning of the time span, and u[end] is the ending
residual[1] = u[1] + pi/2 # the solution at the beginning of the time span should be -pi/2
residual[2] = u[end] - pi/2 # the solution at the end of the time span should be pi/2
end
bvp2 = TwoPointBVProblem(simplependulum!, bc2!, [pi/2,pi/2], tspan)
sol2 = solve(bvp2, MIRK4(), dt=0.05) # we need to use the MIRK4 solver for TwoPointBVProblem
plot(sol2)
```
Note that `u` is a tuple of `( u[1], u[end] )` just like `t` is `( t[1], t[end] )` and `p` holds the parameters of the given problem.

![BVP Example Plot2](../assets/bvp_example_plot2.png)

27 changes: 25 additions & 2 deletions docs/src/tutorials/ode_example.md
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,29 @@ plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")

where the pieces are described below.

#### Note

If that code errors, it is most likely due to versioning issues.

```julia
Pkg.status("DifferentialEquations")
```

should return a version higher than `v4.0.0`. If not, `Pkg.update()` to
get the latest versions. If that does not work, then there is some other
package upper bounding the allowed version of DifferentialEqautions.jl.
The changes are described at length here:

http://juliadiffeq.org/2018/01/24/Parameters.html

If you don't wish to update right now, simply do

```julia
Pkg.pin("DifferentialEquations",v"3.1.0")
```

and use the sidebar to change the documentation to v3.2.0.

### Step 1: Defining a Problem

To solve this numerically, we define a problem type by giving it the equation,
Expand Down Expand Up @@ -169,13 +192,13 @@ Convenience features are also included. We can build an array using a
comprehension over the solution tuples via:

```julia
[t+u for (u,p,t) in tuples(sol)]
[t+u for (u,t) in tuples(sol)]
```

or more generally

```julia
[t+2u for (u,p,t) in zip(sol.u,sol.t)]
[t+2u for (u,t) in zip(sol.u,sol.t)]
```

allows one to use more parts of the solution type. The object that is returned by
Expand Down

0 comments on commit fc7118b

Please sign in to comment.