Skip to content
This repository has been archived by the owner on Sep 9, 2024. It is now read-only.

Drafted some general API documentation #28

Merged
merged 1 commit into from
Mar 11, 2014
Merged

Drafted some general API documentation #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 33 additions & 0 deletions doc/api.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
#General API for all solvers
Please open pull requests or issues to propose changes or clarifications.

##Basic interface
The general interface for all ODE solvers is:

```julia
t_out, y_out = odeXX(F, y0, tspan; kwargs...)
```

Each (adaptive) solver accepts 3 arguments

- `F`: the RHS of the ODE `dy/dt = F(t,y)`, which is a function of `t` and `y(t)` and returns `dy/dt::typeof(y./t)`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why y./t rather than y/t? Division by a scalar should always be defined (though we could also write it as y*inv(t)), and I think we should specify that t will be a scalar.

(In principle, you could have an interface where a vector of times are passed, but this is more important in languages that rely on vectorization for performance. It greatly simplifies the implementation of F if you know that t will always be a scalar.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because of JuliaLang/julia#5810. If you think "Division by a scalar should always be defined", you must raise your concern there. (Or maybe I'm just plaint wrong, in which case you should tell me here)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A/α is still defined. JuliaLang/julia#5810 only deprecated α/A.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I was confused in that regard. Thanks for the clarification. I'll change this.

- `y0`: initial value for `y`. The type of `y0` determines the element type of the `y_out` vector (`y_out::Vector{typeof(y0)}`)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Technically, this may not be typeof(y0) due to promotion. Maybe typeof(y0*one(t))?

We should also specify that y0 should be of a type that supports +, -, * by scalars (we could also require / by scalars, but this is not strictly needed), as well as the norm function below.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good points.

I was not sure what to do with the promotion issues, because I apparently lost the discussion in #26. I'll let others decide how that should be specified.

I want to be explicit about the operations the ODE solver variables should support, but I was not sure I could provide a complete list, so I left it out for now.

- `tspan`: Any iterable of sorted `t` values at which the solution (`y`) is requested. Most solvers will only consider `tspan[1]` and `tspan[end]`, and intermediary points will be interpolated. If `tspan[1] > tspan[end]` the integration is performed backwards.

The solver returns two arrays

- `tout`: Vector of points at which solutions were obtained (also see keyword `points`)
- `yout`: solutions at times `tout`; if `points=:all | :specified` is used, `yout::Vector{typeof(y0)}`. Note that if `y0` is a vector, you can get a matlab-like matrix with `hcat(yout...)`.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought we had agreed that yout would always be the same type regardless of the keywords, for type stability?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought that too. I have removed :final as a valid value, but it seems like the bullet point was left in a horrible broken state afterwards. I'll fix that.


Each solver might implement its own keywords, but the following keywords have a standardized interpretation across all solvers. Solvers should raise an error if a unrecognized keyword argument is used.

- `norm`: user-supplied norm for determining the error `E` (default `Base.norm`)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would suggest defaulting to base.vecnorm, similar to quadgk in Julia 0.3. This has two advantages:

  • It is defined for any iterable type, including multidimensional arrays. (norm is only defined for vectors and matrices).
  • norm(A) is slow for matrices (where it requires the SVD).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like a reasonable choice. (If we should not just leave the default unspecified, and let the specific algorithm decide the most reasonable norm.)

- `abstol` and/or `reltol`: an integration step is accepted if `E <= abstol || E <= reltol*abs(y)` (ideally we want both criteria for all solvers, **done** in #13)
- `points=:all | :specified`: controls the type of output according to
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If there are only two possible values of points, wouldn't it be cleaner to just have allpoints::Bool=true or false?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are we sure that there will never be a need to extend points to allow other specifications?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At some point we considered :event, but we probably don't need that anymore? MATLAB has an option called Refine, which allows you to increases the number of output points by some factor. I guess they are using dense output to calculate those intermediate solutions. I don't know if we want something like that, but I thought I should mention it. We could use points to pass the factor to the solver.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think points = :all | :specified is a much better API, even if we never extend the choices of points beyond this. If I read someone else's code, it's much easier to, at a glance, realize what points = :specified does, than allpoints=false - the second option doesn't convey any information about what the alternative to allpoints=true is.

* `points==:all` (default) output is given for each value in `tspan` as well as for each intermediate point the solver used.
* `points==:specified` output is given only for each value in `tspan`.
- `maxstep`, `minstep` and `initstep`: determine the maximal, minimal and initial integration step.
- `retries = 0` Sometimes an integration step takes you out of the region where `F(t,y)` has a valid solution and F might throw `DomainError` or other exceptions. `retries` sets a limit to the number of times the solver might try with a smaller step.

##Iterator interface
Under construction #27
10 changes: 10 additions & 0 deletions doc/solvers.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
#Currently implemented ODE solvers

##ODE45
* Status?
* Special considerations


##ODE23
* Status
* Special considerations