diff --git a/.gitignore b/.gitignore new file mode 100644 index 000000000..a6b2f497b --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +build +docs/build +docs/site diff --git a/.travis.yml b/.travis.yml index f00154aa1..a28bc851d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,8 +3,7 @@ os: - osx - linux julia: - - 0.3 - - release + - 0.5 - nightly git: depth: 999999 @@ -12,6 +11,9 @@ notifications: email: false script: - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi - - julia -e 'Pkg.clone(pwd()); Pkg.build("ODE"); Pkg.test("ODE"; coverage=true)'; + - julia -e 'Pkg.clone("https://github.com/JuliaODE/ODE.jl.git"); Pkg.build("ODE"); Pkg.test("ODE"; coverage=true)'; + after_success: - julia -e 'cd(Pkg.dir("ODE")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder()); Codecov.submit(process_folder())'; + - julia -e 'Pkg.add("Documenter")' + - julia -e 'cd(Pkg.dir("ODE")); include(joinpath("docs", "make.jl"))' diff --git a/README.md b/README.md index 323b44c71..39f99e199 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,16 @@ Various basic Ordinary Differential Equation solvers implemented in Julia. -[![Build Status](https://travis-ci.org/JuliaLang/ODE.jl.svg?branch=master)](https://travis-ci.org/JuliaLang/ODE.jl) -[![Coverage Status](https://img.shields.io/coveralls/JuliaLang/ODE.jl.svg)](https://coveralls.io/r/JuliaLang/ODE.jl) -[![ODE](http://pkg.julialang.org/badges/ODE_0.3.svg)](http://pkg.julialang.org/?pkg=ODE&ver=0.3) -[![ODE](http://pkg.julialang.org/badges/ODE_0.4.svg)](http://pkg.julialang.org/?pkg=ODE&ver=0.4) +[![Join the chat at https://gitter.im/pwl/ODE.jl](https://badges.gitter.im/pwl/ODE.jl.svg)](https://gitter.im/pwl/ODE.jl?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) +[![Build Status](https://travis-ci.org/JuliaODE/ODE.jl.svg?branch=master)](https://travis-ci.org/JuliaODE/ODE.jl) +[![Coverage Status](https://img.shields.io/coveralls/JuliaODE/ODE.jl.svg)](https://coveralls.io/r/JuliaODE/ODE.jl) +[![ODE](http://pkg.julialang.org/badges/ODE_0.5.svg)](http://pkg.julialang.org/?pkg=ODE&ver=0.5) + +[![](https://img.shields.io/badge/docs-latest-blue.svg)](https://JuliaODE.github.io/ODE.jl/latest) Pull requests are always highly welcome to fix bugs, add solvers, or anything else! # API discussions + There are currently discussions about how the Julian API for ODE solvers should look like, and the current documentation is more like a wishlist than a documentation. The API has changed considerably since the initial v0.1 release, so be carefull when you upgrade to v0.2 or later versions. # Current status of the project diff --git a/REQUIRE b/REQUIRE index 8e8352dde..ffc6c5b9f 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,3 +1,4 @@ -julia 0.3 +julia 0.5- Polynomials +ForwardDiff Compat 0.4.1 diff --git a/doc/api.md b/doc/api.md deleted file mode 100644 index ae94aad44..000000000 --- a/doc/api.md +++ /dev/null @@ -1,35 +0,0 @@ -#General API for all solvers -**This is a working draft for how the API might look like in the future** - -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)` -- `y0`: initial value for `y`. The type of `y0`, promoted as necessary according to the numeric type used for the times, determines the element type of the `yout` vector (`yout::Vector{typeof(y0*one(t))}`) -- `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 times are promoted as necessary to a common floating-point type. - -The solver returns two arrays - -- `tout`: Vector of points at which solutions were obtained (also see keyword `points`) -- `yout`: solutions at times `tout`, stored as a vector `yout` as described above. Note that if `y0` is a vector, you can get a matlab-like matrix with `hcat(yout...)`. - -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.vecnorm`) -- `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 - * `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 diff --git a/doc/solvers.md b/doc/solvers.md deleted file mode 100644 index 735018de8..000000000 --- a/doc/solvers.md +++ /dev/null @@ -1,10 +0,0 @@ -#Currently implemented ODE solvers - -##ODE45 -* Status? -* Special considerations - - -##ODE23 -* Status -* Special considerations \ No newline at end of file diff --git a/docs/.documenter.enc b/docs/.documenter.enc new file mode 100644 index 000000000..12d46f201 Binary files /dev/null and b/docs/.documenter.enc differ diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 000000000..3df633a3c --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,25 @@ +using Documenter, ODE + +makedocs( + format = Documenter.Formats.HTML, + modules = [ODE], + clean = false, + sitename = "ODE.jl", + pages = [ + "Home" => "index.md", + + "Manual" => [ + "Basics" => "man/basics.md", + "Base" => "man/base.md" + ] + ] + + ) + +deploydocs( + repo = "github.com/JuliaODE/ODE.jl.git", + target = "build", + julia = "0.5", + deps = Deps.pip("pygments", "python-markdown-math"), + make = nothing + ) diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 000000000..ea3c6fdeb --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,218 @@ +```@contents +Pages = [ + "tutorials/euler_integrator.md", + "man/basics.md", + "man/base.md" + ] +``` + +# ODE.jl + +## Top level interface + +If you are looking to getting a solution without the additional hussle +of handling an iterator we provide the wrappers `ODE.odeXX`. They +provide a simplistic way of handling explicit differential equations +of the form `y'=F(t,y)` with `y` being either a `Number` or an +`AbstractArray` of any dimension. Below we solve a simple initial +value problem given by `y'=y` with initial data at `t0=0.0` and +`y0=1.0` on the interval `[t0,1]`. + +```@example ode +using ODE +tspan = [0.0,1.0] +y0 = 1.0 +F(t,y) = y +(t,y) = ODE.ode(F,y0,tspan) +``` + +The vectors `t` and `y` store the time and solution values at the +corresponding times. + +## Solve interface + +`ODE.ode` only supports explicit differential equations defined as +`y'=F(t,y)`, for more advenced uses consider using `ODE.solve`, which +was designed to work with a variety of other types of initial value +problems and is optimized for better performance. First we have to +define an initial value problem, in our case this is an explicit +differential equation `y'=y` with inital data `y0=[1.0]` given at the +time `t0=0.0`. + +```@example solve +using ODE +t0 = 0.0 +y0 = [1.0] +F!(t,y,dy) = dy[1]=y[1] +ode = ODE.ExplicitODE(t0,y0,F!) +``` + +Note that unlike in `ODE.ode` we now have to supply an in place +function `F!` instead of an explicit function `F`. We can solve the +ODE problem `ode` by simply calling + +```@example solve +sol = ODE.solve(ode, tstop = 1) +``` + +This returns a `Solution` type, which stores the solution. You +probably noticed that we passed a keyword argument `tstop`, this is +the final time of integration which we have to specify because `tstop` +defaults to `Inf` and the integration would carry on forever. You can +access the solution with + +```@example solve +(t,y) = sol.t, sol.y +``` + +You can change the default algorithm (Runge-Kutta (4,5)) by passing an +optional argument `solver` + +```@example solve +sol = ODE.solve(ode, tstop = 1, solver = ODE.RKIntegratorAdaptive{:dopri5}) +``` + +For other options accepted by `solve` see [Options](/Options/) below. + +You might still find this interface limiting. First of all, it stores +all the results, so if you are only interested in the final value of +`y` it still stores all the intermediate steps. Secondly, you cannot +process the results on the fly (e.g. plot the current state of a +solution). If you need more control you should consider using the +iterator interface. + +## Iterator interface + +To offeset the limitations of the `ODE.ode` interface we implemented a +general. We use the same problem as before as an example + +```@example iterate +using ODE +t0 = 0.0 +y0 = [1.0] +F!(t,y,dy) = dy[1]=y[1] +ode = ODE.ExplicitODE(t0,y0,F!) +``` + +Now we have full flow control over the solver, we can analyze the +intermediate results or interrupt the integration at any point. + +```@example iterate +for (t,y) in ODE.iterate(ode) + @show (t,y) + if t > 1 + break + end +end +``` + +Note that we had to break the loop because `sol` would keep producing +the results. To set the final integration time and other parameters +of the integrator `integ` we can pass optional arguments to +`ODE.solver`. + +```@example iterate +for (t,y) in ODE.iterate(ode; tstop = 1) + @show (t,y) +end +``` + +This approach has the added benefit of the solution never exceeding +the final time. Both `ODE.iterate` and `ODE.solve` support the same +options, so you can easily change the method of integration with the +keyword `solver`. + +Apart from the time and value `(t,y)` the `ODE.solve` also returns the +derivative, you can retrive it as the third argument in the returned +tuple. In the following example we use it compute the absolute +residual error (zero in this case). + +```@example iterate +for (t,y,dy) in ODE.iterate(ode; tstop = 1) + err = norm(y-dy) + @show err +end +``` + +With `tstop` specified we can also get all results at once using +`collect` and other constructs working on iterators +(e.g. generators). For example + +```@example iterate +iter = ODE.iterate(ode; tstop = 1) +solution = collect(iter) +``` + +returns a vector of triples `(t,y,dy)`. Or if you only wan the first +component of a solution you could simply use + +```@example iterate +y1 = collect(y[1] for (t,y) in iter) +``` + +There are, however, several caveats that you should take into account: + +1. Each time the iterator is collected the differential equation is + actually solved, which has potentially high computational cost and + might be inefficient. + +2. The `length` is undefined for the result of `ODE.iterate`, because + we don't know a priori how many steps the integration will require + (especially in the case of adaptive solvers). This means that the + functions requireing `length` might not work. For the same reason + there are no `getindex` methods. + +## Options + +Both `ODE.ode` and `ODE.solve` accept the following keyword arguments. + +- `integ`: the type of integrator to use, defaults to a adaptive + Runge-Kutta method of order 4/5. To see the list of available + integrators see [`Integrators`](@ref). + +- `initstep`: The initial step size, defaults to `eps(T)^(1/3)`. + +- `tstop`: The final integration time, never exceeded by the + integrator. In case of `ODE.ode(F,y0,tspan)` this option defaults + to the last element of `tspan` if it is a vector. In `ODE.solve` + the default is `tstop=Inf`. If `tstop` is smaller then `t0` the + integration runs backwards in time. + +Apart from these general options, each integrator has its own keyword +arguments. In particular all integrators with adaptive step size +can be cotrolled with + +- `reltol`, `abstol`: The relative and absolute error tolerances. The + solution guarantees that at each step we have + `norm((y-yc)*reltol.+abstol)<=1`, where `yc` is a true solution to + and IVP. Defaults are `reltol=eps(T)^(1/3)/10`, + `abstol=eps(T)^(1/2)/10`. + +- `norm`: The norm used to measure error in the formula above, + defaults to `y->Base.vecnorm(y,Inf)`. You can specify it to assign + different weights to different components of `y`. + +- `minstep`, `maxstep`: Minimal and maximal stepsize for the + integrator. If at any point the stepsize exceeds these limits the + integrator will yield an error and cease producing + results. Deafaults are `minstep=10*eps(T)` and `maxstep=1/minstep`. + +- `maxiters`: The number of iterations before the integrator ceases to + work, defaults to `Inf`. Useful as a safeguard from iterator + continuing ad infinitum. + +- `isoutofdomain`: Applied to each component of `y`, if + `isoutofdomain(y[i])==true` the integrator stops. Defaults to + `Base.isnan`. + +Apart from these, each integrator may support additional options. + +## Integrators + +### Explicit Runge-Kutta integrators + +### Rosenbrock methods + +### Backwards differential formula (BDF) methods + +### ??? diff --git a/docs/src/man/base.md b/docs/src/man/base.md new file mode 100644 index 000000000..535f2443e --- /dev/null +++ b/docs/src/man/base.md @@ -0,0 +1,46 @@ +```@meta +CurrentModule = ODE +``` + +# Base + +The file `base.jl` implements the most basic iterator infrastructure +for solvers and the definitions of the types representing general IVP +(initial value problem) and solvers. + +## General functions for solving initial value problems + +```@docs +solve +iterate +``` + +## Predefined types of initial value problems + +```@docs +AbstractIVP +IVP +ExplicitODE +ImplicitODE +``` + +## Solver architecture + +```@docs +AbstractSolver +AbstractIntegrator +AbstractState +Solution +``` + +The fallback constructor for `AbstractSolver(ivp::IVP;opts...)` ensures +that an error is thrown if a solver is constructed for an unsupported +type of the given IVP. + +## Fallback functions for solvers + +```@docs +Base.length(::AbstractSolver) +output(::AbstractState) +Base.eltype{T,Y}(::Type{AbstractIVP{T,Y}}) +``` diff --git a/docs/src/man/basics.md b/docs/src/man/basics.md new file mode 100644 index 000000000..18248a870 --- /dev/null +++ b/docs/src/man/basics.md @@ -0,0 +1,3 @@ +# Basic usage + +Consider an ODE $y'=y$ diff --git a/docs/src/tutorials/euler_integrator.md b/docs/src/tutorials/euler_integrator.md new file mode 100644 index 000000000..09d1d5214 --- /dev/null +++ b/docs/src/tutorials/euler_integrator.md @@ -0,0 +1,165 @@ +Below you will find the simplest implementation of a reasonable +generic solver that finds solutions to explicit ODE equations. + +```julia +type EulerIntegrator <: AbstractIntegrator + # nothing in here yet +end + +type EulerState + t + y + dy +end + +output(state::EulerState) = state.t, state.y, state.dy + +# see we don't need a fancy constructor +function solver(ode::ExplicitODE, + ::Type{EulerIntegrator}; + options...) + return Problem(ode,EulerIntegrator()) +end + +function init(ode::ExplicitODE, integ::EulerIntegrator) + t0, y0 = ode.t0, ode.y0 + dy0 = similar(ode.dy0) + ode.F!(t0,y0,dy0) # fill in the values of the derivative + EulerState(t0,y0,dy0) +end + +function onestep!(ode::ExplicitODE, integ::EulerIntegrator, state::EulerState) + t, y, dy = output(state) + dt = 0.1 + y += dt*dy + t += dt + ode.F!(t0,y0,dy0) # update the derivative + return cont +end +``` + +There are several problems with the above implementation. First of +all, it has a constant prescribed step size. This could easily be +fixed by changing the type definition and the `solver` to + + +```julia +type EulerIntegrator <: AbstractIntegrator + initstep +end + +function solver(ode::ExplicitODE, + ::Type{EulerIntegrator}; + initstep = 0.1, + options...) + return Problem(ode,EulerIntegrator(initstep)) +end +``` + +we should also change the line `dt = 0.1` in the `onestep!` function +to `dt = stepper.initstep`. Now we can run our integrator with a +custom step size! + +```julia +sol = solver(ode,EulerIntegrator,initstep = 0.01) +for (t,y) in sol + if t > 1 + print(t,y) + break + end +end +``` + +Another issue is type stability, to make `EulerIntegrator` perform +better we should explicitly annotate the fields in both +`EulerIntegrator` and `EulerState` like this + +```julia +type EulerIntegrator{T,Y} <: AbstractIntegrator{T,Y} + initstep::T +end + +type EulerState{T,Y} + t::T + y::Y + dy::Y +end + +function solver(ode::ExplicitODE{T,Y}, + ::Type{EulerIntegrator}; + initstep::T = T(0.1), + options...) + return Problem(ode,EulerIntegrator{T,Y}(initstep)) +end +``` + +But the `EulerState{T,Y}` is exactly the same as `Step` from +`base.jl`, so we can simplify it a bit more + +```julia +type EulerState{T,Y} + step::Step{T,Y} +end +``` + +Once we do that, in the `init` we should change +`EulerState(t0,y0,dy0)` to `EulerState(Step(t0,y0,dy0))` and redefine +`output` to `output(state::EulerState)=output(state.step)` +(`output(::Step)` is already defined in `base.jl`). + +One could even replace `EulerState` with `Step` completely, but this +won't allow us to extend the state with some additional variables and +storage space in the future. + +The last thing is that our stepper will continue the integration +forever: it doesn't have a stopping condition. We could add one as an +option. + +```julia +type EulerIntegrator{T,Y} <: AbstractIntegrator{T,Y} + initstep::T + tstop::T +end + +function solver(ode::ExplicitODE{T,Y}, + ::Type{EulerIntegrator}; + initstep::T = T(0.1), + tstop::T = T(Inf) + options...) + return Problem(ode,EulerIntegrator{T,Y}(initstep,tstop)) +end + +function onestep!(ode::ExplicitODE, integ::EulerIntegrator, state::EulerState) + t, y, dy = output(state) + + if t > integ.tstop + return finished + end + + dt = integ.initstep + y += dt*dy + t += dt + ode.F!(t0,y0,dy0) # update the derivative + return cont +end +``` + +As a final improvement, we can (although this is not necessary) use a +structure `FixedOptions` from `options.jl` to keep our options in one +structure. A corresponding options type for adaptive solver is +`AdaptiveOptions`. This way we can use the standarized defaults for +most options and keep our solver in line with the standard +keywords. Naturally, we have to update `onestep!` to use the subtype. + +```julia +type EulerIntegrator{T,Y} <: AbstractIntegrator{T,Y} + options::FixedOptions{T,Y} +end + +function solver(ode::ExplicitODE{T,Y}, + ::Type{EulerIntegrator}; + options...) + options = FixedOptions{T}(;options...) + return Problem(ode,EulerIntegrator{T,Y}(options)) +end +``` diff --git a/examples/custom_integrator.jl b/examples/custom_integrator.jl new file mode 100644 index 000000000..af2f136e9 --- /dev/null +++ b/examples/custom_integrator.jl @@ -0,0 +1,127 @@ +include("../src/ODE.jl") + +""" + +Here we demonstrate how to implement an integrator so that it is +compatible with `ODE.jl`. The implementation can be contained in an +external package, but if the API conforms to our backend one could use +this integrator to solve an explicit ODE with `solve(...)`. + +""" + +module MyIntegrator + +using ODE +# We have to define these methods on our integrator +import ODE: init, onestep! +import ODE: AbstractIntegrator, AbstractState, ExplicitODE, Step + +# The options are stored in the integrator type. Integrator is +# immutable: we can't change the options once we start the +# integration. +immutable EulerIntegrator{T} <: AbstractIntegrator{T} + tstop::T + initstep::T +end + +# This is necessary for the `solve` to construct the iterator. Every +# iterator has to have a constructor of the form +# `Iterator(::IVP;opts...)`, here we implement an interator that only +# works with explicit differential equations, hence we restrict the +# constructor to `ExplicitODE`. If someone tries to use our iterator +# to solve an equation of unsupported type `solve` will throw an error +# (there is a fallback constructor for `AbstractIntegrator` that +# throws an error). +function EulerIntegrator{T}(ode::ExplicitODE{T}; + tstop = T(Inf), + initstep = T(1//10), + opts...) + EulerIntegrator{T}(tstop,initstep) +end + +# The state of the integrator, it stores the current values for time +# `t` and the solution `y` itself, along with its derivative, +# `dy`. For convenience we used the already available type to store +# this kind of information, `Step`, that has three fields: `t`, `y` +# and `dy`. See the end of the module for an alternative definition +# of a state. +type EulerState{T,Y} <: AbstractState{T,Y} + step::ODE.Step{T,Y} +end +output(state::EulerState) = output(state.step) + +# Generates the state given an ODE and an integrator, this method has +# to be specialized for `EulerIntegrator`, we don't have to specialize +# it for `ExplicitODE` as the type of an ODE is already filtered by +# the specialized the construtor, but we do it here for clarity. +function init(ode::ExplicitODE, integ::EulerIntegrator) + EulerState(ODE.Step(ode.t0,copy(ode.y0),copy(ode.dy0))) +end + +function onestep!(ode::ExplicitODE, integ::EulerIntegrator, state::EulerState) + # as mentioned before, this function unpacks the variables + # `t,y,dy` from the current state. It is not necessary, you could + # access the fields directly but we use it here for convenience. + t, y, dy = ODE.output(state) + + tdir = sign(integ.tstop-ode.t0) + + # the only stop condition our solver has. Note the use of `abs`, + # which enables integration backward in time. + if tdir*t >= tdir*integ.tstop + # this flag finalizes the iterator + return ODE.finish + else + # trim the stepsize to match the `tstop`, prevents + # overshooting + dt = tdir*min(integ.initstep,abs(integ.tstop-t)) + + # update the time, + state.step.t += dt + # the function (this is the `dy` from the previous step) + state.step.y += dt*dy + # and its derivative + ode.F!(t,y,dy) + # return a flag to continue the integration + return ODE.cont + end +end + +# OPTIONAL: Define properties of this integrator: order, name and +# whether it is adaptive or not. At this point the information +# supplied here is not used but it might be a good idea to implement +# these methods for future use. +order{T}(::Type{EulerIntegrator{T}}) = 1 +name{T}(::Type{EulerIntegrator{T}}) = "My own Euler integrator" +isadaptive{T}(::Type{EulerIntegrator{T}}) = false + +# OPTIONAL: +# Another possiblity to implement state would be to declare +type EulerState2{T,Y} <: ODE.AbstractState{T,Y} + t::T + y::Y + dy::Y +end +# but then we have to define the method `output` +output(state::EulerState2) = (state.t, state.y, state.dy) + +end + +# Usage example +using ODE +# import ODETests: test_integrator +using ODETests +using MyIntegrator + +integ = MyIntegrator.EulerIntegrator + +# declare the ODE as usual +ode =ODE.ExplicitODE(0.0,[1.0],(t,y,dy)->copy!(dy,y)) +# solve the `ode` with our integrator, note that we can pass options to `solve` +sol =ODE.solve(ode,integ;tstop=1.0,initstep=0.001) +# # print the last step of the solution +collect(sol)[end] + +# test the integrator +# TODO: for now I can't figure out why it fails. +ODETests.test_integrator(integ) diff --git a/examples/ex_iter.jl b/examples/ex_iter.jl new file mode 100644 index 000000000..3847d2a2b --- /dev/null +++ b/examples/ex_iter.jl @@ -0,0 +1,65 @@ +using ODE + +# pick your solver +integ = [ODE.RKIntegratorAdaptive{:rk45}, + ODE.ModifiedRosenbrockIntegrator][2] + +# Define IVP-instance which holds the mathematical problem definition: +t0 = 0.0 +y0 = [1.0] +ode = ODE.ExplicitODE(t0,y0,(t,y,dy)->dy[1]=y[1]) + + +### Forward time integration + +# options for the solver +opts = Dict(:initstep=>0.1, + :tstop=>1.0, + :reltol=>1e-5, + :abstol=>1e-5) + +# create a Problem instance +prob = ODE.solve(ode,integ;opts...) + + # iterate over the solution +println("t, y, err") +for (t,y) in prob + println((t,y[1],abs(y[1]-e.^t))) +end + +# or collect it +println(collect(prob)) + +### Reverse time integration, rest as above +t0 = 1.0 +y0 = [1.0] +ode = ODE.ExplicitODE(t0,y0,(t,y,dy)->dy[1]=y[1]) +opts = Dict(:initstep=>0.1, + :tstop=>0.0, + :reltol=>1e-5, + :abstol=>1e-5) + +prob = ODE.solve(ode,integ;opts...) + +println("t, y, err") +for (t,y) in prob # iterate over the solution + println((t,y[1],abs(y[1]-e.^(t-1)))) +end + +println(collect(prob)) + + +### Dense output +opts = Dict(:initstep=>0.1, + :tstop=>0.0, + :reltol=>1e-5, + :abstol=>1e-5, + :tout=>[t0:-0.1:0;]) + +prob = ODE.solve(ode,ODE.DenseOutput{integ};opts...) +println("t, y, err") +for (t,y) in prob # iterate over the solution + println((t,y[1],abs(y[1]-e.^(t-1)))) +end + +println(collect(prob)) diff --git a/examples/iterator-implementation.jl b/examples/iterator-implementation.jl new file mode 100644 index 000000000..68c469d66 --- /dev/null +++ b/examples/iterator-implementation.jl @@ -0,0 +1,173 @@ +type EulerIntegrator <: AbstractIntegrator + # nothing in here yet +end + +type EulerState + t + y + dy +end + +output(state::EulerState) = state.t, state.y, state.dy + +# see we don't need a fancy constructor +function solver(ode::ExplicitODE, + ::Type{EulerIntegrator}; + options...) + return Problem(ode,EulerIntegrator()) +end + +function init(ode::ExplicitODE, integ::EulerIntegrator) + t0, y0 = ode.t0, ode.y0 + dy0 = similar(ode.dy0) + ode.F!(t0,y0,dy0) # fill in the values of the derivative + EulerState(t0,y0,dy0) +end + +function onestep!(ode::ExplicitODE, integ::EulerIntegrator, state::EulerState) + t, y, dy = output(state) + dt = 0.1 + y += dt*dy + t += dt + ode.F!(t0,y0,dy0) # update the derivative + return cont +end + + +""" + +There are several problems with the above implementation. First of +all, it has a constant prescribed step size. This could easily be +fixed by changing the type definition and the `solver` to + +""" + +type EulerIntegrator <: AbstractIntegrator + initstep +end + +function solver(ode::ExplicitODE, + ::Type{EulerIntegrator}; + initstep = 0.1, + options...) + return Problem(ode,EulerIntegrator(initstep)) +end + +""" + +we should also change the line `dt = 0.1` in the `onestep!` function +to `dt = stepper.initstep`. Now we can run our integrator with a +custom step size! + +""" + +sol = solver(ode,EulerIntegrator,initstep = 0.01) +for (t,y) in sol + if t > 1 + print(t,y) + break + end +end + +""" + +Another issue is type stability, to make `EulerIntegrator` perform +better we should explicitly annotate the fields in both +`EulerIntegrator` and `EulerState` like this + +""" + +type EulerIntegrator{T,Y} <: AbstractIntegrator{T,Y} + initstep::T +end + +type EulerState{T,Y} + t::T + y::Y + dy::Y +end + +function solver(ode::ExplicitODE{T,Y}, + ::Type{EulerIntegrator}; + initstep::T = T(0.1), + options...) + return Problem(ode,EulerIntegrator{T,Y}(initstep)) +end + + +""" + +But the `EulerState{T,Y}` is exactly the same as `Step` from +`base.jl`, so we can simplify it a bit more + +""" + +type EulerState{T,Y} + step::Step{T,Y} +end + +""" + +Once we do that, in the `init` we should change +`EulerState(t0,y0,dy0)` to `EulerState(Step(t0,y0,dy0))` and redefine +`output` to `output(state::EulerState)=output(state.step)` +(`output(::Step)` is already defined in `base.jl`). + +One could even replace `EulerState` with `Step` completely, but this +won't allow us to extend the state with some additional variables and +storage space in the future. + +The last thing is that our stepper will continue the integration +forever: it doesn't have a stopping condition. We could add one as an +option. + +""" + +type EulerIntegrator{T,Y} <: AbstractIntegrator{T,Y} + initstep::T + tstop::T +end + +function solver(ode::ExplicitODE{T,Y}, + ::Type{EulerIntegrator}; + initstep::T = T(0.1), + tstop::T = T(Inf) + options...) + return Problem(ode,EulerIntegrator{T,Y}(initstep,tstop)) +end + +function onestep!(ode::ExplicitODE, integ::EulerIntegrator, state::EulerState) + t, y, dy = output(state) + + if t > integ.tstop + return finished + end + + dt = integ.initstep + y += dt*dy + t += dt + ode.F!(t0,y0,dy0) # update the derivative + return cont +end + +""" + +As a final improvement, we can (although this is not necessary) use a +structure `FixedOptions` from `options.jl` to keep our options in one +structure. A corresponding options type for adaptive solver is +`AdaptiveOptions`. This way we can use the standarized defaults for +most options and keep our solver in line with the standard +keywords. Naturally, we have to update `onestep!` to use the subtype. + +""" + +type EulerIntegrator{T,Y} <: AbstractIntegrator{T,Y} + options::FixedOptions{T,Y} +end + +function solver(ode::ExplicitODE{T,Y}, + ::Type{EulerIntegrator}; + options...) + options = FixedOptions{T}(;options...) + return Problem(ode,EulerIntegrator{T,Y}(options)) +end diff --git a/examples/test.jl b/examples/test.jl new file mode 100644 index 000000000..fd078100e --- /dev/null +++ b/examples/test.jl @@ -0,0 +1,34 @@ +include("../src/ODE.jl") + +module Test + +using ODE + +T = Float64 +Y = Vector{T} +t0 = zero(T) +y0 = T[one(T)] + +solvers = [# ODE.RKIntegratorAdaptive{:rk45}, + # ODE.RKIntegratorFixed{:feuler}, + ODE.DenseOutput] + +for st in solvers + ode = ODE.ExplicitODE(t0,y0,(t,y,dy)->dy[1]=y[1]) + opts = Dict(:initstep=>0.1, + :tspan=>[0.,0.5,1.], + :points=>:specified, + :reltol=>1e-5, + :abstol=>1e-5) + + prob = ODE.solve(ode,st;opts...) + + println("Raw iterator") + for (t,y) in prob + println((t,y,norm(y-[exp(t)]))) + end + + println(collect(prob)) +end + +end diff --git a/src/ODE.jl b/src/ODE.jl index e8ccc55c1..73de9a648 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -1,428 +1,78 @@ -isdefined(Base, :__precompile__) && __precompile__() # Ordinary Differential Equation Solvers -module ODE - -using Polynomials -using Compat - -## minimal function export list -# adaptive non-stiff: -export ode23, ode45, ode78 -# non-adaptive non-stiff: -export ode4, ode4ms -# adaptive stiff: -export ode23s -# non-adaptive stiff: -export ode4s - -## complete function export list: see runtests.jl - -############################################################################### -## Coefficient Tableaus -############################################################################### - -# Butcher Tableaus, or more generally coefficient tables -# see Hairer & Wanner 1992, p. 134, 166 - -abstract Tableau{Name, S, T<:Real} -# Name is the name of the tableau/method (a symbol) -# S is the number of stages (an int) -# T is the type of the coefficients -# -# TODO: have a type parameter which specifies adaptive vs non-adaptive -# -# For all types of tableaus it assumes fields: -# order::(Int...) # order of the method(s) -# -# For Runge-Kutta methods it assumes fields: -# a::Matrix{T} # SxS matrix -# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive) -# c::Vector{T} # S -# -# For a tableau: -# c1 | a_11 .... a_1s -# . | a_21 . . -# . | a_31 . . -# . | .... . . -# c_s | a_s1 ....... a_ss -# -----+-------------------- -# | b_1 ... b_s this is the one used for stepping -# | b'_1 ... b'_s this is the one used for error-checking - -Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T -order(b::Tableau) = b.order -# Subtypes need to define a convert method to convert to a different -# eltype with signature: -Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert method for concrete Tableau types") - -############################################################################### -## HELPER FUNCTIONS -############################################################################### - -# estimator for initial step based on book -# "Solving Ordinary Differential Equations I" by Hairer et al., p.169 -function hinit(F, x0, t0, tend, p, reltol, abstol) - # Returns first step, direction of integration and F evaluated at t0 - tdir = sign(tend-t0) - tdir==0 && error("Zero time span") - tau = max(reltol*norm(x0, Inf), abstol) - d0 = norm(x0, Inf)/tau - f0 = F(t0, x0) - d1 = norm(f0, Inf)/tau - if d0 < 1e-5 || d1 < 1e-5 - h0 = 1e-6 - else - h0 = 0.01*(d0/d1) - end - # perform Euler step - x1 = x0 + tdir*h0*f0 - f1 = F(t0 + tdir*h0, x1) - # estimate second derivative - d2 = norm(f1 - f0, Inf)/(tau*h0) - if max(d1, d2) <= 1e-15 - h1 = max(1e-6, 1e-3*h0) - else - pow = -(2. + log10(max(d1, d2)))/(p + 1.) - h1 = 10.^pow - end - return tdir*min(100*h0, h1, tdir*(tend-t0)), tdir, f0 -end - -# isoutofdomain takes the state and returns true if state is outside -# of the allowed domain. Used in adaptive step-control. -isoutofdomain(x) = isnan(x) - -function make_consistent_types(fn, y0, tspan, btab::Tableau) - # There are a few types involved in a call to a ODE solver which - # somehow need to be consistent: - # - # Et = eltype(tspan) - # Ey = eltype(y0) - # Ef = eltype(Tf) - # - # There are also the types of the containers, but they are not - # needed as `similar` is used to make containers. - # Tt = typeof(tspan) - # Ty = typeof(y0) # note, this can be a scalar - # Tf = typeof(F(tspan(1),y0)) # note, this can be a scalar - # - # Returns - # - Et: eltype of time, needs to be a real "continuous" type, at - # the moment a AbstractFloat - # - Eyf: suitable eltype of y and f(t,y) - # --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1])) - # - Ty: container type of y0 - # - btab: tableau with entries converted to Et - - # Needed interface: - # On components: /, - - # On container: eltype, promote_type - # On time container: eltype - - Ty = typeof(y0) - Eyf = typeof(y0[1]/(tspan[end]-tspan[1])) - - Et = eltype(tspan) - @assert Et<:Real - if !(Et<:AbstractFloat) - Et = promote_type(Et, Float64) - end - - # if all are Floats, make them the same - if Et<:AbstractFloat && Eyf<:AbstractFloat - Et = promote_type(Et, Eyf) - Eyf = Et - end - - !isleaftype(Et) && warn("The eltype(tspan) is not a concrete type! Change type of tspan for better performance.") - !isleaftype(Eyf) && warn("The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.") - - btab_ = convert(Et, btab) - return Et, Eyf, Ty, btab_ -end - -############################################################################### -## NON-STIFF SOLVERS -############################################################################### - -include("runge_kutta.jl") - -# ODE_MS Fixed-step, fixed-order multi-step numerical method -# with Adams-Bashforth-Moulton coefficients -function ode_ms(F, x0, tspan, order::Integer) - h = diff(tspan) - x = Array(typeof(x0), length(tspan)) - x[1] = x0 - - if 1 <= order <= 4 - b = ms_coefficients4 - else - b = zeros(order, order) - b[1:4, 1:4] = ms_coefficients4 - for s = 5:order - for j = 0:(s - 1) - # Assign in correct order for multiplication below - # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) - p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1]))) - b[s, s - j] = ((-1)^j / factorial(j) - / factorial(s - 1 - j) * polyval(p_int, 1)) - end - end - end - - # TODO: use a better data structure here (should be an order-element circ buffer) - xdot = similar(x) - for i = 1:length(tspan)-1 - # Need to run the first several steps at reduced order - steporder = min(i, order) - xdot[i] = F(tspan[i], x[i]) - - x[i+1] = x[i] - for j = 1:steporder - x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)] - end - end - return vcat(tspan), x -end - -# Use order 4 by default -ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4) -ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5) +""" +Coding conventions: -############################################################################### -## STIFF SOLVERS -############################################################################### +- use `t,y,dy`, use type/function parameters `T` and `Y` +- `p::Problem`, use parameter `P` +- `ivp::IVP`, use parameter `O` + - if referring to a ODE or DAE, use `ode` or `dae` instead +- `integ::AbstractIntegrator`, use parameter `I` +- `opts::AbstactOptions`, , use parameter `OP` -# Crude forward finite differences estimator of Jacobian as fallback +Variables and Type variables: +- T -> t::T +- Y -> y::Y -# FIXME: This doesn't really work if x is anything but a Vector or a scalar -function fdjacobian(F, x::Number, t) - ftx = F(t, x) - # The 100 below is heuristic - dx = (x .+ (x==0))./100 - dFdx = (F(t,x+dx)-ftx)./dx +""" +module ODE - return dFdx -end +using Polynomials +using Compat +import Compat.String +using ForwardDiff -function fdjacobian(F, x, t) - ftx = F(t, x) - lx = max(length(x),1) - dFdx = zeros(eltype(x), lx, lx) - for j = 1:lx - # The 100 below is heuristic - dx = zeros(eltype(x), lx) - dx[j] = (x[j] .+ (x[j]==0))./100 - dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j] - end - return dFdx -end +import Base: start, next, done, collect, show, convert -# ODE23S Solve stiff systems based on a modified Rosenbrock triple -# (also used by MATLAB's ODE23s); see Sec. 4.1 in +# Core infrastructure # -# [SR97] L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22 +# When wrapping a new solver it will need to use and conform to +# methods and types within these files. # -# supports keywords: points = :all | :specified (using dense output) -# jacobian = G(t,y)::Function | nothing (FD) -function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, - jacobian=nothing, - points=:all, - norm=Base.norm, - minstep=abs(tspan[end] - tspan[1])/1e18, - maxstep=abs(tspan[end] - tspan[1])/2.5, - initstep=0.) - - - # select method for computing the Jacobian - if typeof(jacobian) == Function - jac = jacobian - else - # fallback finite-difference - jac = (t, y)->fdjacobian(F, y, t) - end - - # constants - const d = 1/(2 + sqrt(2)) - const e32 = 6 + sqrt(2) - - - # initialization - t = tspan[1] - - tfinal = tspan[end] - - h = initstep - if h == 0. - # initial guess at a step size - h, tdir, F0 = hinit(F, y0, t, tfinal, 3, reltol, abstol) - else - tdir = sign(tfinal - t) - F0 = F(t,y0) - end - h = tdir * min(abs(h), maxstep) - - y = y0 - tout = Array(typeof(t), 1) - tout[1] = t # first output time - yout = Array(typeof(y0), 1) - yout[1] = deepcopy(y) # first output solution - +# Note, if we go the MathProgBase.jl route, then these files would go +# into ODEBase.jl. +include("base.jl") +include("tableaus.jl") +include("options.jl") +include("helpers.jl") - J = jac(t,y) # get Jacobian of F wrt y +# Dense output wrapper +include("dense.jl") - while abs(t - tfinal) > 0 && minstep < abs(h) - if abs(t-tfinal) < abs(h) - h = tfinal - t - end +# Particular solvers +include("integrators/ode23s.jl") +include("integrators/runge-kutta.jl") +include("integrators/adams-bashford-moulton.jl") +include("integrators/rosenbrock.jl") - if size(J,1) == 1 - W = one(J) - h*d*J - else - # note: if there is a mass matrix M on the lhs of the ODE, i.e., - # M * dy/dt = F(t,y) - # we can simply replace eye(J) by M in the following expression - # (see Sec. 5 in [SR97]) +# User interface to solvers +include("top-interface.jl") - W = lufact( eye(J) - h*d*J ) - end - # approximate time-derivative of F - T = h*d*(F(t + h/100, y) - F0)/(h/100) +""" - # modified Rosenbrock formula - k1 = W\(F0 + T) - F1 = F(t + 0.5*h, y + 0.5*h*k1) - k2 = W\(F1 - k1) + k1 - ynew = y + h*k2 - F2 = F(t + h, ynew) - k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T ) +This submodule contains simple test functions for solvers/integrators +compatible with ODE.jl. You can use it to test your custom solvers, +for examples of how to use these functions see our tests in `test/` +directory. - err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate - delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error +""" +module ODETests - # check if new solution is acceptable - if err <= delta +import Base.Test: @test, @test_approx_eq_eps - if points==:specified || points==:all - # only points in tspan are requested - # -> find relevant points in (t,t+h] - for toi in tspan[(tspan.>t) & (tspan.<=t+h)] - # rescale to (0,1] - s = (toi-t)/h +using ODE - # use interpolation formula to get solutions at t=toi - push!(tout, toi) - push!(yout, y + h*( k1*s*(1-s)/(1-2*d) + k2*s*(s-2*d)/(1-2*d))) - end - end - if (points==:all) && (tout[end]!=t+h) - # add the intermediate points - push!(tout, t + h) - push!(yout, ynew) - end +import ODE: AbstractIntegrator, AbstractState, ExplicitODE - # update solution - t = t + h - y = ynew +include("tests/minimal_types.jl") +include("tests/test_cases.jl") - F0 = F2 # use FSAL property - J = jac(t,y) # get Jacobian of F wrt y - # for new solution - end +# Test function for the iterator interface +include("tests/integrators.jl") - # update of the step size - h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) ) - end - - return tout, yout -end - - -#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method -# with provided coefficients. -function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing) - - if typeof(jacobian) == Function - G = jacobian - else - G = (t, x)->fdjacobian(F, x, t) - end - - h = diff(tspan) - x = Array(typeof(x0), length(tspan)) - x[1] = x0 - - solstep = 1 - while solstep < length(tspan) - ts = tspan[solstep] - hs = h[solstep] - xs = x[solstep] - dFdx = G(ts, xs) - # FIXME - if size(dFdx,1) == 1 - jac = 1/gamma/hs - dFdx[1] - else - jac = eye(dFdx)/gamma/hs - dFdx - end - - g = Array(typeof(x0), size(a,1)) - g[1] = (jac \ F(ts + b[1]*hs, xs)) - x[solstep+1] = x[solstep] + b[1]*g[1] - - for i = 2:size(a,1) - dx = zero(x0) - dF = zero(x0/hs) - for j = 1:i-1 - dx += a[i,j]*g[j] - dF += c[i,j]*g[j] - end - g[i] = (jac \ (F(ts + b[i]*hs, xs + dx) + dF/hs)) - x[solstep+1] += b[i]*g[i] - end - solstep += 1 - end - return vcat(tspan), x end -# Kaps-Rentrop coefficients -const kr4_coefficients = (0.231, - [0 0 0 0 - 2 0 0 0 - 4.452470820736 4.16352878860 0 0 - 4.452470820736 4.16352878860 0 0], - [3.95750374663 4.62489238836 0.617477263873 1.28261294568], - [ 0 0 0 0 - -5.07167533877 0 0 0 - 6.02015272865 0.1597500684673 0 0 - -1.856343618677 -8.50538085819 -2.08407513602 0],) - -ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian) - -# Shampine coefficients -const s4_coefficients = (0.5, - [ 0 0 0 0 - 2 0 0 0 - 48/25 6/25 0 0 - 48/25 6/25 0 0], - [19/9 1/2 25/108 125/108], - [ 0 0 0 0 - -8 0 0 0 - 372/25 12/5 0 0 - -112/125 -54/125 -2/5 0],) - -ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian) - -# Use Shampine coefficients by default (matching Numerical Recipes) -const ode4s = ode4s_s - -const ms_coefficients4 = [ 1 0 0 0 - -1/2 3/2 0 0 - 5/12 -4/3 23/12 0 - -9/24 37/24 -59/24 55/24] - - end # module ODE diff --git a/src/base.jl b/src/base.jl new file mode 100644 index 000000000..0a154d51c --- /dev/null +++ b/src/base.jl @@ -0,0 +1,484 @@ +# The main types: +# - IVP -- holds the mathematical aspects of a IVP +# - AbstractIntegrator -- an integrator/solver (maybe AbstractIntegrator?) +# - Problem -- holds IVP + Integrator (maybe ProblemSpec, Problem, Spec?) +# - AbstractState -- holds the iterator state +# - Step -- holds the state at one time +# - + + +""" + AbstractIVP{T,Y} + +The abstract supertype of all IVPs (initial value problems). The type +parameters `T` and `Y` correspond to the types of time and state +variable respectively. + +""" +abstract AbstractIVP{T,Y} + +""" +The "elements" of AbstractIVP are `t,y,dy` so the `eltype` returns `Tuple{T,Y,Y}` +""" +Base.eltype{T,Y}(::Type{AbstractIVP{T,Y}}) = Tuple{T,Y,Y} + +""" + + IVP{T,Y,F,G,J} <: AbstractIVP{T,Y} + +Defines the mathematical part of an IVP (initial value problem) +specified in the general form: + +`F(t, y) = G(t, y, dy)` with `y(t0)= y0` + +Depending on the combination of the parameters this type can represent +a wide range of problems, including ODE, DAE and IMEX. Nevertheless +not all solvers will support any combinations of `F` and `G`. Note +that not specifying `G` amounts to `G=dy/dt`. + + +- `tspan` -- tuple `(start_t,end_t)` +- `y0` -- initial condition +- `F!` -- in-place `F` function `F!(t,y,res)`. If `F=0` set to `nothing`. +- `G!` -- in-place `G` function `G!(t,y,dy,res)`. If `G=dy/dt` then + set to `nothing` (or `dy` if the solver supports this). Can + also be a mass matrix for a RHS `M dy/dt` +- `J!` -- in-place Jacobian function `J!(t,y,dy,res)`. + +TODO: how to fit the sparsity pattern in J? + +""" +type IVP{T,Y,F,G,J} <: AbstractIVP{T,Y} + t0 ::T + y0 ::Y + dy0 ::Y + F! ::F + G! ::G + J! ::J +end + + +""" + + typealias ExplicitODE{T,Y,F,J} IVP{T,Y,F,Void,J} + +Can be constructed by calling + + ODE.ExplicitODE(t0,y0,F!;J!=jacobian)) + +Explicit ODE representing the problem + +`dy = F(t,y)` with `y(t0)=y0` + +- t0, y0: initial conditions +- F!: in place version of `F` called by `F!(t,y,dy)` +- J!: (optional) computes `J=dF/dy` in place, called with `J!(t,y,J)` + +""" +typealias ExplicitODE{T,Y,F,J} IVP{T,Y,F,Void,J} +@compat function (::Type{ExplicitODE}){T,Y}(t0::T, + y0::Y, + F!::Function; + J!::Function = forward_jacobian!(F!,copy(y0)), + kargs...) + # precompute y' + dy0 = copy(y0) + F!(t0,y0,dy0) + IVP(t0,y0,dy0,F!,nothing,J!) +end + +""" + +Implicit ODE representing the problem + +`G(t,y,dy)=0` with `y(t0)=y0` and optionally `y'(t0)=dy0` + +- t0, y0: initial conditions +- G!: in place version of `G` called by `G!(res,t,y,dy)`, + returns residual in-place in `res`. +- J!: (optional) computes `J=dF/dy+a*dF/dy'` for prescribed `a`, called with `J!(out,t,y,dy,a)`. + Returns Jacobian in-place in `out`. + +""" +typealias ImplicitODE{T,Y,G,J} IVP{T,Y,Void,G,J} +@compat function (::Type{ImplicitODE}){T,Y}(t0::T, + y0::Y, + G!::Function; + J!::Function = forward_jacobian_implicit!(G!,similar(y0)), + dy0::Y = zero(y0), + kargs...) + IVP(t0,y0,dy0,nothing,G!,J!) +end + +""" + + AbstractSolver + +The supertype of anything which can get you to a solution of a IVP. +Subtypes include: `AbstractIntegrator`s but also `DenseOutput` + +""" +abstract AbstractSolver + +Base.length(s::AbstractSolver) = error("`length` is not defined for $(typeof(s)).") + +""" + +The fallback generator of an abstract solver, throws an error if the +solver cannot be generated for the given initial value problem type. + +""" +@compat (::Type{S}){S<:AbstractSolver}(ivp;opts...) = + error("The solver $S doesn't support IVP of form $(typeof(ivp))") + + +""" + +The abstract type of the actual algorithm to solve an IVP. + +""" +abstract AbstractIntegrator{T} <: AbstractSolver + + +""" + +AbstractState keeps the temporary data (state) for the iterator +Problem{::AbstractIntegrator}. + +""" +abstract AbstractState{T,Y} + +""" +Returns variables returned during iterations. + +output(st::AbstractState) = t,y,dy +""" +output(st::AbstractState) = st.step.t, st.step.y, st.step.dy + + +""" + +Holds a value of a function and its derivative at time t. This is +usually used to store the solution of an IVP at particular times. + +""" +type Step{T,Y} + t ::T + y ::Y + dy::Y +end + +output(s::Step) = s.t, s.y, s.dy + +function show(io::IO, state::Step) + println("t =$(state.t)") + println("y =$(state.y)") + println("dy =$(state.dy)") +end + + +""" + +This is an iterable type, each call to next(...) produces a next step +of a numerical solution to an IVP. + +- ivp: is the prescrived ivp, along with the initial data +- solver: the algorithm used to produce subsequent values from the ivp + +""" +immutable Problem{O<:AbstractIVP,S<:AbstractSolver} + ivp ::O + solver ::S +end + +Base.length(prob::Problem) = length(prob.solver) + +Base.eltype{O,S}(::Type{Problem{O,S}}) = eltype(O) + +""" +Makes some generic operations on iterators work, like +generator comprehensions: + + tgen=(t for (t,y) in sol) + tout=collect(tgen) + +or + errgen=(y-[exp(t)] for (t,y) in sol) + errout=collect(errgen) +""" +Base.iteratorsize{O,S}(::Type{Problem{O,S}}) = Base.SizeUnknown() + +""" + iterate(ivp::IVP; solver=RKIntegratorAdaptive{:rk45}, opts...) + +Iterate creates an iterable `Problem` instance from an `IVP` instance +(specifying the math) and from a `Type{AbstractSolver}` (the numerical +integrator). The simplest use case is + + for (t,y,dy) in iterate(...) + # do something with t, y an dy + end + +If the integration interval, defined by the keyword argument `tstop`, +is finite you can request all the results at once by calling + + collect(iterate(...)) # => Vector{Tuple{T,Y,Y}} + +Notes: + +- usually solvers require the ivp to be in a certain form, say an + `ExplicitODE`. +- the second argument is the *Type* of the solver and not an instance. + The instance of the solve can only be created together with the + `ivp` as their type parameters need to match. + +Input: + +- `ivp::IVP` +- `S::Type{AbstractSolver}` + +Output: + +- `::Problem` + +""" +function iterate{S<:AbstractSolver}(ivp::IVP; + solver::Type{S} = RKIntegratorAdaptive{:rk45}, + opts...) + Problem(ivp, solver(ivp; opts...)) +end + +""" + solve(ivp::IVP; solver=RKIntegratorAdaptive{:rk45}, opts...) + +Solve the initial value problem `ivp` using an algorithm `solver` +(defaults to Runge-Kutta (4,5) integrator). One can pass additional +options to the `solver` via keyword arguments to `solve` (here denoted +as `options`). The output is a `Solution` type (currently simply a +tuple of vectors `(Vector{T},Vector{Y})`, where `T,Y=eltype(ivp)`). + +""" + +function solve(ivp::IVP; opts...) + + # TODO: perhaps there is a more + # graceful way to treat these + # cases. We only care about + # infinite `tstop` but if + # `tstop` is unspecified it + # defaults to `tout[end]`, with + # `tout` defaulting to + # `[Inf]`. Maybe we should add + # `tout(::Problem)` to fix this? + # Or maybe we could store + # `tstop` in `Problem`. + + dopts = Dict(opts) + if !in(:tstop,keys(dopts)) & !in(:tout,keys(dopts)) + error("Neither `tstop` nor `tout` was specified.") + end + if in(:tstop,keys(dopts)) & !isfinite(dopts[:tstop]) + error("Trying to integrate over an infinite time span, try specifying `|tstop|sol.t[end]) | (t t>ti, sol.t) + theta = (t-sol.t[i])/(sol.t[i+1]-sol.t[i]) + return sol.y[i]+theta*(sol.y[i+1]-sol.y[i]) + end +end + + +# In Julia 0.5 the collect needs length to be defined, we cannot do +# that for a Problem but we can implement our own collect +function collect(prob::Problem) + pairs = Array(eltype(prob),0) + for (t,y,dy) in prob + push!(pairs,(t,copy(y),copy(dy))) + end + return pairs +end + + +# Iteration: take one step on a IVP `Problem` +# +# Defines: +# start(iter) -> state +# next(iter, state) -> output(state), state +# done(iter, state) -> bool +# +# Perhaps unintuitively, the next step is computed in `done`. Such +# implementation allows to decide if the iterator is exhausted in case +# when the next step was computed but it was deemed incorrect. In +# such situation `done` returns `false` after computing the step and +# the failed step never sees the light of the day (by not being +# returned by `next`). +# +# TODO: this implementation fails to return the zeroth step (t0,y0) +# +# TODO: store the current Step outside of the actual state + +Base.start(prob::Problem) = init(prob.ivp, prob.solver) + +function Base.done(prob::Problem, st) + # Determine whether the next step can be made by calling the + # stepping routine. onestep! will take the step in-place. + status = onestep!(prob.ivp, prob.solver, st) + if status==cont + return false + elseif status==finish + return true + else #if status==abort + warn("aborting") + return true + # else + # error("unsported Status: $status") + end +end + +function Base.next(prob::Problem, st) + # Output the step (we know that `done` allowed it, so we are safe + # to do it) + return output(st), st +end + +""" + +Holds the solver status, used inside of `onestep!`. + +Values: + +- cont -- continue integration +- abort -- abort integration +- finish -- integration reached the end + +Statuses can be combined with &: +- cont&cont == cont +- finish&cont == finish +- abort&cont == abort +- abort&finish = abort + +""" +@enum Status cont=1 abort=0 finish=-1 +# The values of Statuses are chose to turn & into a *: +@compat Base.:&(s1::Status, s2::Status) = Status(Int(s1)*Int(s2)) + +##### +# Interface to implement by solvers to hook into iteration +##### +# +# See runge_kutta.jl and rosenbrock.jl for example implementations. + +# A stepper has to implement +# - init +# - output +# and either +# - onestep! +# - trialstep!, errorcontrol! and accept! + +""" + +Take a step, modifies `state` in-place. This is the core function to +be implemented by a solver. However, if possible solvers should opt +to implement the sub-step functions `trialstep!`, `errorcontrol!` and +`accept!`, instead of directly `onestep!`. + +Input: + +- prob::Problem, state::AbstractState + +Output: + +- Status + +substeps. + +""" +function onestep!(ivp::IVP, integ::AbstractIntegrator, state::AbstractState) + opt = integ.opts + while true + status = trialstep!(ivp, integ, state) + err, status_err = errorcontrol!(ivp, integ, state) + status &= status_err + if err<=1 + # a successful step + status &= accept!(ivp, integ, state) + return status + elseif status==abort || status==finish + return status + end + # if we get here: try step again with updated state (step + # size, order). + end +end + + +""" + +Advances the solution by trying to compute a single step. The new +step is kept in the `state` in work arrays so that `errorcontrol!` can +compute the magnitude of its error. If the error is small enough +`accept!` updates `state` to reflect the state at the new time. + +Returns `Status`. + +""" +trialstep!{I<:AbstractIntegrator}(::IVP, ::I, ::AbstractState) = + error("Function `trialstep!` and companions (or alternatively `onestep!`) need to be implemented for adaptive integrator $I") + +""" + +Estimates the error (such that a step is accepted if err<=1). +Depending on the stepper it may update the state, e.g. by computing a +new dt or a new order (but not by computing a new solution!). + +Returns `(err,Status)`. + +If the `status==abort` then the integration is aborted, status values +of `cont` and `finish` are ignored. + +""" +errorcontrol!{I<:AbstractIntegrator}(::IVP, ::I, ::AbstractState) = + error("Function `errorcontrol!` and companions (or alternatively `onestep!`) need to be implemented for adaptive integrator $I") + +""" + +Accepts (in-place) the computed step. Called if `errorcontrol!` gave +a small enough error. + +Returns `Status`. + +""" +accept!{I<:AbstractIntegrator}(::IVP, ::I, ::AbstractState) = + error("Function `accept!` and companions (or alternatively `onestep!`) need to be implemented for adaptive integrator $I") diff --git a/src/dense.jl b/src/dense.jl new file mode 100644 index 000000000..dcb3df18b --- /dev/null +++ b/src/dense.jl @@ -0,0 +1,233 @@ +# A higher level solver, defined as a wrapper around an integrator. + +""" + +Dense output options: + +- tout ::Vector{T} output times + +""" + +immutable DenseOptions{T<:Number,TO<:AbstractVector} <: Options{T} + tout::TO + + # Planned options: + # - points ::Symbol which points are returned: `:specified` only the + # ones in tspan or `:all` which includes also the step-points of the solver. + # - stopevent Stop integration at a zero of this function + # - roottol +end + +@compat function (::Type{DenseOptions{T}}){T}(; + tstop = T(1//0), + tout::AbstractVector{T} = T[tstop], + # points::Symbol= :all, + # stopevent::S = (t,y)->false, + # roottol = eps(T)^T(1//3), + kargs...) + DenseOptions{T,typeof(tout)}(tout) +end + + +""" + +A solver specialized in dense output; it wraps an integrator. It +stores the subsequent steps generated by `Problem` and interpolates +the results (currently this means at the output times stored in +`opts.tout`). + +""" +immutable DenseOutput{I<:AbstractSolver,OP<:DenseOptions} <: AbstractSolver + integ::I + opts::OP +end + +Base.length(dense::DenseOutput) = length(dense.opts.tout) + +@compat function (::Type{DenseOutput{I}}){T,I}(ivp::IVP{T}; + tstop = T(1//0), + tout::AbstractVector{T} = T[tstop], + opts...) + if all(tout.>=ivp.t0) + tout = sort(tout) + elseif all(tout.<=ivp.t0) + tout = reverse(sort(tout)) + else + error("Elements of tout should all be either to the right or to the left of `t0`.") + end + + # normalize `tstop` to the last element of `tout` + tstop = tout[end] + + # create integrator + integ = I(ivp; tstop=tstop, opts...) + + # create dense solver + dense_opts = DenseOptions{T}(; tout=tout, opts...) + dense_solver = DenseOutput(integ, dense_opts) + return dense_solver +end + +""" + +The state of the dense solver `DenseOutput`. + +""" +type DenseState{St<:AbstractState,T,Y} <: AbstractState{T,Y} + tout_i::Int + step_prev::Step{T,Y} + step_out::Step{T,Y} + integrator_state::St +end + + +output(ds::DenseState) = output(ds.step_out) + + +function init(ivp::IVP, + solver::DenseOutput) + integrator_state = init(ivp, solver.integ) + dy0 = copy(ivp.y0) + ivp.F!(ivp.t0,ivp.y0,dy0) + step_prev = Step(ivp.t0,copy(ivp.y0),dy0) + step_out = Step(ivp.t0,copy(ivp.y0),copy(ivp.y0)) + return DenseState(1,step_prev,step_out,integrator_state) +end + + +function onestep!(ivp::IVP, + solver::DenseOutput, + dstate::DenseState) + i = dstate.tout_i + if i > length(solver.opts.tout) + return finish + end + + # the underlying integrator + integ = solver.integ + + # our next output time + ti = solver.opts.tout[i] + + istate = dstate.integrator_state + + + # try to get a new set of steps enclosing `ti`, if all goes + # right we end up with t∈[t1,t2] with + # t1,_=output(dstate.step_prev) + # t2,_=output(dstate.integrator_state) + status = next_interval!(ivp, integ, istate, dstate.step_prev, ti) + if status == abort + # we failed to get enough steps + warn("Iterator was exhausted before the dense output could produce the output.") + return abort + else + # we got the steps, proceed with the interpolation, this fills + # the dstate.step_out with y(ti) and y'(ti) according to an + # interpolation algorithm specific for a method (defaults to + # hermite O(3)). + interpolate!(istate, dstate.step_prev, ti, dstate.step_out) + + # increase the counter + dstate.tout_i += 1 + return cont + end +end + +""" + +Takes steps using the underlying integrator until it reaches a first +step such that `t>=tout`. It fills the `steps` variable with +(Step(t1,y(t1),dy(t1)),Step(t2,y(t2),dy(t2))), where `t1` is is the +step before `tout` and `t2` is `>=tout`. In other words +`tout∈[t1,t2]`. +""" +function next_interval!(ivp, integ, istate, step_prev, tout) + while true + # get the current time + t1 = step_prev.t + t2,_ = output(istate) + + # in case we are integrating backwards in time reverse the + # time interval + if t2 < t1 + t1, t2 = t2, t1 + end + + if t1 <= tout <= t2 + # we found the enclosing times + return cont + end + + # save the current state of solution + t, y, dy = output(istate) + step_prev.t = t + copy!(step_prev.y,y) + copy!(step_prev.dy,dy) + + # try to perform a single step: + status = onestep!(ivp, integ, istate) + + if status != cont + return status + end + end + + # this will never happen + return abort +end + + +""" +Makes dense output + +interpolate!(istate::AbstractState, step_prev::Step, tout, step_out::Step) + +Input: + +- `istate::AbstractState` state of the integrator +- `step_prev` the previous step, part of `dstate` +- tout -- time of requested output +- step_out::Step -- inplace output step + +Output: nothing + +TODO: output dy too + +TOOD: provide arbitrary order dense output. Maybe use work of @obiajulu on A-B-M methods. +""" +function interpolate! end + +""" +Make dense output using Hermite interpolation of order O(3), should +work for most integrators and is used as default. This only needs y +and dy at t1 and t2. + +Ref: Hairer & Wanner p.190 +""" +function interpolate!(istate::AbstractState, + step_prev::Step, + tout, + step_out::Step) + t1,y1,dy1 = output(step_prev) + t2,y2,dy2 = output(istate) + if tout==t1 + copy!(step_out.y,y1) + elseif tout==t2 + copy!(step_out.y,y2) + else + dt = t2-t1 + theta = (tout-t1)/dt + for i=1:length(y1) + step_out.y[i] = + (1-theta)*y1[i] + + theta*y2[i] + + theta*(theta-1) * + ( (1-2*theta)*(y2[i]-y1[i]) + + (theta-1)*dt*dy1[i] + + theta*dt*dy2[i]) + end + end + step_out.t = tout + return nothing +end diff --git a/src/helpers.jl b/src/helpers.jl new file mode 100644 index 000000000..28afc685f --- /dev/null +++ b/src/helpers.jl @@ -0,0 +1,51 @@ +## TODO: reactivate when introducing events/rootfinding +# """ + +# A simple bisection algorithm for finding a root of a solution f(x)=0 +# starting within the range x∈rng, the result is a point x₀ which is +# located within the distance eps from the true root of f(x)=0. For +# this algorithm to work we need f(rng[1]) to have a different sign then +# f(rng[2]). + +# """ +# function findroot(f,rng,eps) +# xl, xr = rng +# fl, fr = f(xl), f(xr) + +# if fl*fr > 0 || xl > xr +# error("Inconsistent bracket") +# end + +# while xr-xl > eps +# xm = (xl+xr)/2 +# fm = f(xm) + +# if fm*fr > 0 +# xr = xm +# fr = fm +# else +# xl = xm +# fl = fm +# end +# end + +# return (xr+xl)/2 +# end + +# generate a jacobian using ForwardDiff +function forward_jacobian(F,y0::AbstractArray) + (t,y)->ForwardDiff.jacobian(y->F(t,y),y) +end + +function forward_jacobian(F,y0) + (t,y)->ForwardDiff.derivative(y->F(t,y),y) +end + +function forward_jacobian!(F!,tmp) + jac!(t,y,J)=ForwardDiff.jacobian!(J,(dy,y)->F!(t,y,dy),tmp,y) + return jac! +end + +function forward_jacobian_implicit!(F!,tmp) + error("Not implemented yet") +end diff --git a/src/integrators/adams-bashford-moulton.jl b/src/integrators/adams-bashford-moulton.jl new file mode 100644 index 000000000..5226dbd5a --- /dev/null +++ b/src/integrators/adams-bashford-moulton.jl @@ -0,0 +1,52 @@ +# ODE_MS Fixed-step, fixed-order multi-step numerical method +# with Adams-Bashforth-Moulton coefficients +function ode_ms{Ty,T}(F, x0::Ty, tspan::AbstractVector{T}, order::Integer; kargs...) + + if !isleaftype(T) + error("The output times have to be of a concrete type.") + elseif !(T <:AbstractFloat) + error("The time variable should be a floating point number.") + end + + if !isleaftype(Ty) & !isleaftype(eltype(Ty)) + error("The initial data has to be of a concrete type (or an array)") + end + + h = diff(tspan) + x = Array(typeof(x0), length(tspan)) + x[1] = x0 + + if 1 <= order <= 4 + b = ms_coefficients4 + else + b = zeros(order, order) + b[1:4, 1:4] = ms_coefficients4 + for s = 5:order + for j = 0:(s - 1) + # Assign in correct order for multiplication below + # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) + p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1]))) + b[s, s - j] = ((-1)^j / factorial(j) + / factorial(s - 1 - j) * polyval(p_int, 1)) + end + end + end + + # TODO: use a better data structure here (should be an order-element circ buffer) + xdot = similar(x) + for i = 1:length(tspan)-1 + # Need to run the first several steps at reduced order + steporder = min(i, order) + xdot[i] = F(tspan[i], x[i]) + + x[i+1] = x[i] + for j = 1:steporder + x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)] + end + end + return vcat(tspan), x +end + +# Use order 4 by default +ode4ms(F, x0, tspan; kargs...) = ode_ms(F, x0, tspan, 4; kargs...) +ode5ms(F, x0, tspan; kargs...) = ODE.ode_ms(F, x0, tspan, 5; kargs...) diff --git a/src/integrators/ode23s.jl b/src/integrators/ode23s.jl new file mode 100644 index 000000000..35414e276 --- /dev/null +++ b/src/integrators/ode23s.jl @@ -0,0 +1,176 @@ +# ODE23S Solve stiff systems based on a modified Rosenbrock triple +# (also used by MATLAB's ODE23s); see Sec. 4.1 in +# +# [SR97] L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22 + +immutable ModifiedRosenbrockIntegrator{T<:Number} <: AbstractIntegrator + opts::AdaptiveOptions{T} + const_d::T + const_e::T +end + +function ModifiedRosenbrockIntegrator{T,Y<:AbstractVector}(ode::ExplicitODE{T,Y};opts...) + const_d = 1/(2+sqrt(T(2))) + const_e = 6+sqrt(T(2)) + + ModifiedRosenbrockIntegrator( AdaptiveOptions{T}(;opts...), const_d, const_e ) +end + +order(::ModifiedRosenbrockIntegrator) = 2 +name(::ModifiedRosenbrockIntegrator) = "Modified Rosenbrock Integrator" +isadaptive(::ModifiedRosenbrockIntegrator) = true +tdir(ode::ExplicitODE, integ::ModifiedRosenbrockIntegrator) = sign(integ.opts.tstop - ode.t0) + +""" +The state for the Rosenbrock integrator + +- step: Last successful step +- F1,F2: Work arrays for storing the intermediate values of y' +- J: Jacobian +- iters: Number of successful steps made + +""" +type RosenbrockState{T,Y,J<:AbstractMatrix} <: AbstractState + step ::Step{T,Y} + dt ::T + F1 ::Y + F2 ::Y + k1 ::Y + k2 ::Y + k3 ::Y + ynew ::Y + dtold::T + jac ::J + iters::Int +end + + +# for debugging +function show(io::IO, state::RosenbrockState) + show(io,state.step) + println("dt =$(state.dt)") + println("F1 =$(state.F1)") + println("F2 =$(state.F2)") + println("jac =$(state.jac)") +end + + +function init{T}(ode::ExplicitODE{T}, + integ::ModifiedRosenbrockIntegrator) + t = ode.t0 + dt = integ.opts.initstep + y = ode.y0 + dy = zero(y) + + jac = Array(eltype(y),length(y),length(y)) + + step = Step(t,copy(y),copy(dy)) + state = RosenbrockState(step, + dt, + zero(y), # F1 + zero(y), # F2 + zero(y), # k1 + zero(y), # k2 + zero(y), # k3 + zero(y), # ynew + zero(dt), # dtnew + jac, # jac + 0) # iters + + # initialize the derivative and the Jacobian + ode.F!(t,y,step.dy) + ode.J!(t,y,state.jac) + + return state +end + + +function trialstep!(ode::ExplicitODE, + integ::ModifiedRosenbrockIntegrator, + state::RosenbrockState) + # unpack + step = state.step + opts = integ.opts + F1, F2, jac = state.F1, state.F2, state.jac + k1,k2,k3,ynew = state.k1, state.k2, state.k3, state.ynew + t, dt, y, dy = step.t, state.dt, step.y, step.dy + F! = ode.F! + F0 = dy + + td = tdir(ode,integ) + + # see whether we're done + if td*t >= td*opts.tstop + # nothing left to integrate + return finish + end + + # increment iteration counter + state.iters += 1 + if state.iters > opts.maxiters + warn("Reached maximum number of iterations $(opts.maxiters)") + return abort + end + + W = lufact!( eye(jac) - dt*integ.const_d*jac ) + + # Approximate time-derivative of F, we are using F1 as a + # temporary array + F!(t+dt/100,y,F1) + tder = 100*integ.const_d*(F1-F0) + + # modified Rosenbrock formula + # TODO: update k1,k2,k3 in-place + k1[:] = W \ (F0 + tder) + F!(t+dt/2, y+dt*k1/2, F1) + k2[:] = W \ (F1 - k1) + k1 + for i=1:length(y) + ynew[i] = y[i] + dt*k2[i] + end + F!(t+dt, ynew, F2) + k3[:] = W \ (F2 - integ.const_e*(k2 - F1) - 2*(k1 - F0) + tder ) + + return cont +end + +function errorcontrol!(ode::ExplicitODE, + integ::ModifiedRosenbrockIntegrator, + state::RosenbrockState) + + step = state.step + opts = integ.opts + k1,k2,k3 = state.k1, state.k2, state.k3 + k1,k2,k3,ynew = state.k1, state.k2, state.k3, state.ynew + t, dt, y, dy = step.t, state.dt, step.y, step.dy + + td = tdir(ode,integ) + + # allowable error + delta = max(opts.reltol*max(opts.norm(y), opts.norm(ynew)),opts.abstol) + + # error estimate + err = (abs(dt)/6)*(opts.norm(k1 - 2*k2 + k3))/delta + + # new step-size + dtnew = td*min(opts.maxstep, abs(dt)*(8//10)*err^(-1//3) ) + + # trim in case newdt > dt + dtnew = td*min(abs(dtnew), abs(opts.tstop-(t+dt))) + + state.dtold = dt + state.dt = dtnew + return err, cont +end + +function accept!(ode::ExplicitODE, + integ::ModifiedRosenbrockIntegrator, + state::RosenbrockState) + step = state.step + # update the state + step.t = step.t+state.dtold + copy!(step.y, state.ynew) + copy!(step.dy, state.F2) + ode.J!(step.t,step.y,state.jac) + + return cont +end diff --git a/src/integrators/rosenbrock.jl b/src/integrators/rosenbrock.jl new file mode 100644 index 000000000..8bcba8e6e --- /dev/null +++ b/src/integrators/rosenbrock.jl @@ -0,0 +1,89 @@ +#ODEROSENBROCK Solve stiff differential equations, Rosenbrock method +# with provided coefficients. +function oderosenbrock{Ty,T}(F, x0::Ty, tspan::AbstractVector{T}, + gamma, a, b, c; + jacobian = forward_jacobian(F,x0), + kargs...) + + if !isleaftype(T) + error("The output times have to be of a concrete type.") + elseif !(T <:AbstractFloat) + error("The time variable should be a floating point number.") + end + + if !isleaftype(Ty) & !isleaftype(eltype(Ty)) + error("The initial data has to be of a concrete type (or an array)") + end + + h = diff(tspan) + x = Array(typeof(x0), length(tspan)) + x[1] = x0 + + solstep = 1 + while solstep < length(tspan) + ts = tspan[solstep] + hs = h[solstep] + xs = x[solstep] + dFdx = jacobian(ts, xs) + # FIXME + if size(dFdx,1) == 1 + jac = 1/gamma/hs - dFdx[1] + else + jac = eye(dFdx)/gamma/hs - dFdx + end + + g = Array(typeof(x0), size(a,1)) + g[1] = (jac \ F(ts + b[1]*hs, xs)) + x[solstep+1] = x[solstep] + b[1]*g[1] + + for i = 2:size(a,1) + dx = zero(x0) + dF = zero(x0/hs) + for j = 1:i-1 + dx += a[i,j]*g[j] + dF += c[i,j]*g[j] + end + g[i] = (jac \ (F(ts + b[i]*hs, xs + dx) + dF/hs)) + x[solstep+1] += b[i]*g[i] + end + solstep += 1 + end + return vcat(tspan), x +end + + +# Kaps-Rentrop coefficients +const kr4_coefficients = (0.231, + [0 0 0 0 + 2 0 0 0 + 4.452470820736 4.16352878860 0 0 + 4.452470820736 4.16352878860 0 0], + [3.95750374663 4.62489238836 0.617477263873 1.28261294568], + [ 0 0 0 0 + -5.07167533877 0 0 0 + 6.02015272865 0.1597500684673 0 0 + -1.856343618677 -8.50538085819 -2.08407513602 0],) + +ode4s_kr(F, x0, tspan; kargs...) = oderosenbrock(F, x0, tspan, kr4_coefficients...; kargs...) + +# Shampine coefficients +const s4_coefficients = (0.5, + [ 0 0 0 0 + 2 0 0 0 + 48/25 6/25 0 0 + 48/25 6/25 0 0], + [19/9 1/2 25/108 125/108], + [ 0 0 0 0 + -8 0 0 0 + 372/25 12/5 0 0 + -112/125 -54/125 -2/5 0],) + +ode4s_s(F, x0, tspan; kargs...) = oderosenbrock(F, x0, tspan, s4_coefficients...; kargs...) + +# Use Shampine coefficients by default (matching Numerical Recipes) +const ode4s = ode4s_s + +const ms_coefficients4 = [ 1 0 0 0 + -1/2 3/2 0 0 + 5/12 -4/3 23/12 0 + -9/24 37/24 -59/24 55/24] diff --git a/src/integrators/runge-kutta.jl b/src/integrators/runge-kutta.jl new file mode 100644 index 000000000..6b1a11053 --- /dev/null +++ b/src/integrators/runge-kutta.jl @@ -0,0 +1,524 @@ +# This file contains the implementation of explicit Runkge-Kutta +# solver from (Hairer & Wanner 1992 p.134, p.165-169). + +########################################### +# Tableaus for explicit Runge-Kutta methods +########################################### +immutable TableauRKExplicit{T} <: Tableau{T} + order::(@compat(Tuple{Vararg{Int}})) # the order of the methods + a::Matrix{T} + # one or several row vectors. First row is used for the step, + # second for error calc. + b::Matrix{T} + c::Vector{T} + isFSAL::Bool + s::Int + name::String + function TableauRKExplicit(name,order,a,b,c) + s = length(c) + @assert c[1]==0 + @assert istril(a) + @assert s==size(a,1)==size(a,2)==size(b,2) + @assert size(b,1)==length(order) + @assert maxabs(sum(a,2)-c'')= td*integ.opts.tstop + # nothing left to integrate + return finish + end + + dof = length(step.y) + b = integ.tableau.b + dt = td*min(abs(state.dt),abs(integ.opts.tstop-step.t)) + + copy!(work.ynew,step.y) + + for k=1:length(b) + calc_next_k!(work, k, ode, integ.tableau, step, dt) + for d=1:dof + work.ynew[d] += dt * b[k]*work.ks[k][d] + end + end + step.t += dt + copy!(step.y,work.ynew) + return cont +end + + +######################## +# Adaptive step method # +######################## + + +const timeout_const = 5 + +# `trialstep!` ends with a step computed for the stepsize `state.dt` +# and stores it in `work.y`, so `work.y` contains a candidate for +# `y(t+dt)` with `dt=state.dt`. +function trialstep!(ode::ExplicitODE, integ::RKIntegratorAdaptive, state::RKState) + work = state.work + step = state.step + tableau = integ.tableau + opts = integ.opts + + td = tdir(ode,integ) + + # use the proposed step size to perform the computations + state.dt = state.newdt + dt = state.dt + + if td*step.t >= td*opts.tstop + # nothing left to integrate + return finish + end + + if abs(dt) < opts.minstep + warn("Minimum step size reached") + return abort + end + + # work.y and work.yerr and work.ks are updated after this step + rk_embedded_step!(work, ode, tableau, step, dt) + + return cont +end + +# computes the error for the candidate solution `y(t+dt)` with +# `dt=state.dt` and proposes a new time step +function errorcontrol!(ode::ExplicitODE, + integ::RKIntegratorAdaptive, + state::RKState) + work = state.work + step = state.step + tableau = integ.tableau + timeout = state.timeout + opts = integ.opts + err, state.newdt, state.timeout = + stepsize_hw92!(work, step, tableau, state.dt, state.timeout, opts) + + td = tdir(ode,integ) + + # trim in case newdt > dt + state.newdt = td*min(abs(state.newdt), abs(opts.tstop-(state.step.t+state.dt))) + + if err > 1 + # The error is too large, the step will be rejected. We reset + # the timeout and set the new stepsize + state.timeout = timeout_const + end + + return err, cont +end + +# Here we assume that trialstep! and errorcontrol! have already been +# called, that is `work.y` holds `y(t+dt)` with `dt=state.dt`, and +# error was small enough for us to keep `y(t+dt)` as the next step. +function accept!(ode::ExplicitODE, + integ::RKIntegratorAdaptive, + state::RKState) + work = state.work + step = state.step + tableau = integ.tableau + + # preload ks[1] for the next step + if tableau.isFSAL + copy!(work.ks[1],work.ks[end]) + else + ode.F!(step.t+state.dt, work.ynew, work.ks[1]) + end + + # Swap bindings of y and ytrial, avoids one copy + step.y, work.ynew = work.ynew, step.y + # state.dt holds the size of the last successful step + step.t += state.dt + + return cont +end + + +########################## +# Lower level algorithms # +########################## + +function rk_embedded_step!(work ::RKWorkArrays, + ode ::ExplicitODE, + tableau ::Tableau, + last_step ::Step, + dt) + # Does one embedded R-K step updating work.ynew, work.yerr and work.ks. + # Assumes that work.ks[:,1] is already calculated! + # Modifies work.y, work.ynew and work.yerr only + + y = last_step.y + dof = length(y) + b = tableau.b + + fill!(work.ynew, zero(eltype(y))) + fill!(work.yerr, zero(eltype(y))) + + for s=1:lengthks(tableau) + # we skip the first step beacause we assume that work.ks[1] is + # already computed + if s > 1 + calc_next_k!(work, s, ode, tableau, last_step, dt) + end + for d=1:dof + work.ynew[d] += b[1,s]*work.ks[s][d] + work.yerr[d] += b[2,s]*work.ks[s][d] + end + end + + for d=1:dof + work.yerr[d] = dt*(work.ynew[d]-work.yerr[d]) + work.ynew[d] = y[d] + dt*work.ynew[d] + end + + return nothing +end + + +function stepsize_hw92!{T}(work, + last_step ::Step, + tableau ::TableauRKExplicit, + dt ::T, + timeout ::Int, + opts ::Options{T}) + # Estimates the error and a new step size following Hairer & + # Wanner 1992, p167 (with some modifications) + # + # If timeout>0 no step size increase is allowed, timeout is + # decremented in here. + # + # Returns the error, newdt and the number of timeout-steps + # + # TODO: + # - allow component-wise reltol and abstol? + + ord = T(minimum(order(tableau))) + timout_after_nan = 5 + # fac = T[0.8, 0.9, (0.25)^(1/(ord+1)), (0.38)^(1/(ord+1))][1] + fac = T(8//10) + facmax = T(5) # maximal step size increase. 1.5-5 + facmin = 1/facmax # maximal step size decrease. ? + dof = length(last_step.y) + + # in-place calculate yerr./tol + for d=1:dof + + # TODO: for some reason calling opts.isoutofdomain + # generates a lot of allocations + + # if opts.isoutofdomain(work.y[d])::Bool + if isnan(work.y[d]) + return T(10), dt*facmin, timout_after_nan + end + + y0 = last_step.y[d] + y1 = work.ynew[d] # the approximation to the next step + sci = (opts.abstol + opts.reltol*max(norm(y0),norm(y1))) + work.yerr[d] ./= sci # Eq 4.10 + end + + err = opts.norm(work.yerr) # Eq. 4.11 + newdt = sign(dt)*min(opts.maxstep, abs(dt)*clamp(fac*(1/err)^(1/(ord+1)),facmin,facmax)) # Eq 4.13 modified + + if timeout > 0 + newdt = sign(dt)*min(abs(newdt), abs(dt)) + timeout -= 1 + end + + return err, newdt, timeout +end + + +# For clarity we pass the RKWorkArrays part of the state separately, +# this is the only part of state that can be changed here +function calc_next_k!(work ::RKWorkArrays, + i ::Int, + ode ::ExplicitODE, + tableau ::Tableau, + last_step ::Step, + dt) + dof = length(last_step.y) + t, a, c = last_step.t, tableau.a, tableau.c + + copy!(work.y,last_step.y) + for j=1:i-1 + for d=1:dof + work.y[d] += dt * work.ks[j][d] * a[i,j] + end + end + ode.F!(t + c[i]*dt, work.y, work.ks[i]) + return nothing +end + +################################### +## Tableaus for explicit RK methods +################################### + +# Fixed step: +const tableaus_rk_explicit = Dict{Symbol,ODE.TableauRKExplicit{Rational{Int}}}() + +tableaus_rk_explicit[:feuler] = + TableauRKExplicit{Rational{Int64}}("Forward Euler",(1,), + zeros(Int,1,1), + [1]', + [0] + ) + +tableaus_rk_explicit[:midpoint] = + TableauRKExplicit{Rational{Int64}}("Midpoint",(2,), + [0 0 + 1//2 0], + [0, 1]', + [0, 1//2] + ) + +tableaus_rk_explicit[:heun] = + TableauRKExplicit{Rational{Int64}}("Heun",(2,), + [0 0 + 1 0], + [1//2, 1//2]', + [0, 1]) + +tableaus_rk_explicit[:rk4] = + TableauRKExplicit{Rational{Int64}}("Runge-Kutta(4)",(4,), + [0 0 0 0 + 1//2 0 0 0 + 0 1//2 0 0 + 0 0 1 0], + [1//6, 1//3, 1//3, 1//6]', + [0, 1//2, 1//2, 1]) + +# Adaptive step: +# Heun Euler https://en.wikipedia.org/wiki/Runge–Kutta_methods +tableaus_rk_explicit[:rk21] = + TableauRKExplicit{Rational{Int64}}("Heun Euler",(2,1), + [0 0 + 1 0], + [1//2 1//2 + 1 0], + [0, 1]) + +# Bogacki–Shampine coefficients +tableaus_rk_explicit[:rk23] = + TableauRKExplicit{Rational{Int64}}("Bogacki-Shampine",(2,3), + [0 0 0 0 + 1/2 0 0 0 + 0 3/4 0 0 + 2/9 1/3 4/9 0], + [7/24 1/4 1/3 1/8 + 2/9 1/3 4/9 0], + [0, 1//2, 3//4, 1] + ) + +# Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method +tableaus_rk_explicit[:rk45] = + TableauRKExplicit{Rational{Int64}}("Fehlberg",(4,5), + [ 0 0 0 0 0 0 + 1//4 0 0 0 0 0 + 3//32 9//32 0 0 0 0 + 1932//2197 -7200//2197 7296//2197 0 0 0 + 439//216 -8 3680//513 -845//4104 0 0 + -8//27 2 -3544//2565 1859//4104 -11//40 0 ], +[25//216 0 1408//2565 2197//4104 -1//5 0 + 16//135 0 6656//12825 28561//56430 -9//50 2//55], +[0, 1//4, 3//8, 12//13, 1, 1//2]) + +# Dormand-Prince https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method +tableaus_rk_explicit[:dopri5] = + TableauRKExplicit{Rational{Int64}}("Dormand-Prince", (5,4), + [0 0 0 0 0 0 0 + 1//5 0 0 0 0 0 0 + 3//40 9//40 0 0 0 0 0 + 44//45 -56//15 32//9 0 0 0 0 + 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0 + 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0 + 35//384 0 500//1113 125//192 -2187//6784 11//84 0], + [35//384 0 500//1113 125//192 -2187//6784 11//84 0 + 5179//57600 0 7571//16695 393//640 -92097//339200 187//2100 1//40], + [0, 1//5, 3//10, 4//5, 8//9, 1, 1] + ) + +# Fehlberg 7(8) coefficients +# Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control". +# National Aeronautics and Space Administration. +tableaus_rk_explicit[:feh78] = + TableauRKExplicit{Rational{Int64}}("Fehlberg(7,8)", (7,8), + [ 0 0 0 0 0 0 0 0 0 0 0 0 0 + 2//27 0 0 0 0 0 0 0 0 0 0 0 0 + 1//36 1//12 0 0 0 0 0 0 0 0 0 0 0 + 1//24 0 1//8 0 0 0 0 0 0 0 0 0 0 + 5//12 0 -25//16 25//16 0 0 0 0 0 0 0 0 0 + 1//20 0 0 1//4 1//5 0 0 0 0 0 0 0 0 + -25//108 0 0 125//108 -65//27 125//54 0 0 0 0 0 0 0 + 31//300 0 0 0 61//225 -2//9 13//900 0 0 0 0 0 0 + 2 0 0 -53//6 704//45 -107//9 67//90 3 0 0 0 0 0 + -91//108 0 0 23//108 -976//135 311//54 -19//60 17//6 -1//12 0 0 0 0 + 2383//4100 0 0 -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0 + 3//205 0 0 0 0 -6//41 -3//205 -3//41 3//41 6//41 0 0 0 + -1777//4100 0 0 -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0], + [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0 + 0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0 41//840 41//840], + [0, 2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1] + ) diff --git a/src/options.jl b/src/options.jl new file mode 100644 index 000000000..6194c103e --- /dev/null +++ b/src/options.jl @@ -0,0 +1,78 @@ +abstract Options{T} + +""" + +Generic options for adaptive ODE solvers. This type has a key-word +constructor which will fill the structure with default values. + +General: + +- initstep ::T initial step size (always positive) +- tstop ::T end integration time +- reltol ::T relative tolerance (m3: could this be a vector?) +- abstol ::T absolute tolerance (m3: could this be a vector?) +- minstep ::T minimal allowed step size (always positive) +- maxstep ::T maximal allowed step size (always positive) +- norm function to calculate the norm in step control +- maxiters ::T maximum number of steps +- isoutofdomain::Function checks if the solution is outside of the allowed domain + +""" +immutable AdaptiveOptions{T,N<:Function,O<:Function} <: Options{T} + tstop::T + reltol::T + abstol::T + minstep::T + maxstep::T + initstep::T + norm::N + maxiters::T + isoutofdomain::O +end + +@compat function (::Type{AdaptiveOptions{T}}){T,N,O}(; + tout = [T(1//0)], + tstop = tout[end], + reltol = eps(T)^T(1//3)/T(10), + abstol = eps(T)^T(1//2)/T(10), + minstep = T(10)*eps(T), + maxstep = 1/minstep, + initstep = eps(T)^T(1//3), + norm::N = y->maxabs(y), + maxiters = T(1//0), + isoutofdomain::O = Base.isnan, + kargs...) + @assert minstep>=T(0) && maxstep>=T(0) && initstep>=T(0) # TODO: move to inner constructor + AdaptiveOptions{T,N,O}(tstop,reltol,abstol,minstep,maxstep,initstep,norm,maxiters,isoutofdomain) +end + +""" + +Generic options for fixed step ODE solvers. This type has a key-word +constructor which will fill the structure with default values. + +General: + +- initstep ::T initial step (always positive) +- tstop ::T end integration time + +""" +immutable FixedOptions{T} <: Options{T} + tstop::T + initstep::T +end + +@compat function (::Type{FixedOptions{T}}){T}(; + tout = [T(1//0)], + tstop = tout[end], + initstep = T(1//100), + kargs...) + @assert initstep>=0 + FixedOptions{T}(tstop,initstep) +end + +function show{T}(io::IO, opts :: Options{T}) + for name in fieldnames(opts) + @printf("%-20s = %s\n",name,getfield(opts,name)) + end +end diff --git a/src/runge_kutta.jl b/src/runge_kutta.jl deleted file mode 100644 index a080602de..000000000 --- a/src/runge_kutta.jl +++ /dev/null @@ -1,480 +0,0 @@ -# Explicit Runge-Kutta solvers -############################## -# (Hairer & Wanner 1992 p.134, p.165-169) - -########################################### -# Tableaus for explicit Runge-Kutta methods -########################################### - -immutable TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T} - order::(@compat(Tuple{Vararg{Int}})) # the order of the methods - a::Matrix{T} - # one or several row vectors. First row is used for the step, - # second for error calc. - b::Matrix{T} - c::Vector{T} - function TableauRKExplicit(order,a,b,c) - @assert isa(S,Integer) - @assert isa(Name,Symbol) - @assert c[1]==0 - @assert istril(a) - @assert S==length(c)==size(a,1)==size(a,2)==size(b,2) - @assert size(b,1)==length(order) - @assert norm(sum(a,2)-c'',Inf)<1e-10 # consistency. - new(order,a,b,c) - end -end -function TableauRKExplicit{T}(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), - a::Matrix{T}, b::Matrix{T}, c::Vector{T}) - TableauRKExplicit{name,length(c),T}(order, a, b, c) -end -function TableauRKExplicit(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), T::Type, - a::Matrix, b::Matrix, c::Vector) - TableauRKExplicit{name,length(c),T}(order, convert(Matrix{T},a), - convert(Matrix{T},b), convert(Vector{T},c) ) -end -conv_field{T,N}(D,a::Array{T,N}) = convert(Array{D,N}, a) -function Base.convert{Tnew<:Real,Name,S,T}(::Type{Tnew}, tab::TableauRKExplicit{Name,S,T}) - # Converts the tableau coefficients to the new type Tnew - newflds = () - @compat for n in fieldnames(tab) - fld = getfield(tab,n) - if eltype(fld)==T - newflds = tuple(newflds..., conv_field(Tnew, fld)) - else - newflds = tuple(newflds..., fld) - end - end - TableauRKExplicit{Name,S,Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way? -end - - -isexplicit(b::TableauRKExplicit) = istril(b.a) # Test whether it's an explicit method -isadaptive(b::TableauRKExplicit) = size(b.b, 1)==2 - - -# First same as last. Means ks[:,end]=ks_nextstep[:,1], c.f. H&W p.167 -isFSAL(btab::TableauRKExplicit) = btab.a[end,:]==btab.b[1,:] && btab.c[end]==1 # the latter is not needed really - -## Tableaus for explicit RK methods -# Fixed step: -const bt_feuler = TableauRKExplicit(:feuler,(1,), Rational{Int64}, - zeros(Int,1,1), - [1]', - [0] - ) -const bt_midpoint = TableauRKExplicit(:midpoint,(2,), Rational{Int64}, - [0 0 - 1//2 0], - [0, 1]', - [0, 1//2] - ) -const bt_heun = TableauRKExplicit(:heun,(2,), Rational{Int64}, - [0 0 - 1 0], - [1//2, 1//2]', - [0, 1]) - -const bt_rk4 = TableauRKExplicit(:rk4,(4,),Rational{Int64}, - [0 0 0 0 - 1//2 0 0 0 - 0 1//2 0 0 - 0 0 1 0], - [1//6, 1//3, 1//3, 1//6]', - [0, 1//2, 1//2, 1]) - -# Adaptive step: -# Heun Euler https://en.wikipedia.org/wiki/Runge–Kutta_methods -const bt_rk21 = TableauRKExplicit(:heun_euler,(2,1), Rational{Int64}, - [0 0 - 1 0], - [1//2 1//2 - 1 0], - [0, 1]) - -# Bogacki–Shampine coefficients -const bt_rk23 = TableauRKExplicit(:bogacki_shampine,(2,3), Rational{Int64}, - [0 0 0 0 - 1/2 0 0 0 - 0 3/4 0 0 - 2/9 1/3 4/9 0], - [7/24 1/4 1/3 1/8 - 2/9 1/3 4/9 0], - [0, 1//2, 3//4, 1] - ) - -# Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method -const bt_rk45 = TableauRKExplicit(:fehlberg,(4,5),Rational{Int64}, - [ 0 0 0 0 0 0 - 1//4 0 0 0 0 0 - 3//32 9//32 0 0 0 0 - 1932//2197 -7200//2197 7296//2197 0 0 0 - 439//216 -8 3680//513 -845//4104 0 0 - -8//27 2 -3544//2565 1859//4104 -11//40 0 ], - [25//216 0 1408//2565 2197//4104 -1//5 0 - 16//135 0 6656//12825 28561//56430 -9//50 2//55], - [0, 1//4, 3//8, 12//13, 1, 1//2]) - -# Dormand-Prince https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method -const bt_dopri5 = TableauRKExplicit(:dopri, (5,4), Rational{Int64}, - [0 0 0 0 0 0 0 - 1//5 0 0 0 0 0 0 - 3//40 9//40 0 0 0 0 0 - 44//45 -56//15 32//9 0 0 0 0 - 19372//6561 -25360//2187 64448//6561 -212//729 0 0 0 - 9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0 - 35//384 0 500//1113 125//192 -2187//6784 11//84 0], - [35//384 0 500//1113 125//192 -2187//6784 11//84 0 - 5179//57600 0 7571//16695 393//640 -92097//339200 187//2100 1//40], - [0, 1//5, 3//10, 4//5, 8//9, 1, 1] - ) - -# Fehlberg 7(8) coefficients -# Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control". -# National Aeronautics and Space Administration. -const bt_feh78 = TableauRKExplicit(:feh78, (7,8), Rational{Int64}, - [ 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2//27 0 0 0 0 0 0 0 0 0 0 0 0 - 1//36 1//12 0 0 0 0 0 0 0 0 0 0 0 - 1//24 0 1//8 0 0 0 0 0 0 0 0 0 0 - 5//12 0 -25//16 25//16 0 0 0 0 0 0 0 0 0 - 1//20 0 0 1//4 1//5 0 0 0 0 0 0 0 0 - -25//108 0 0 125//108 -65//27 125//54 0 0 0 0 0 0 0 - 31//300 0 0 0 61//225 -2//9 13//900 0 0 0 0 0 0 - 2 0 0 -53//6 704//45 -107//9 67//90 3 0 0 0 0 0 - -91//108 0 0 23//108 -976//135 311//54 -19//60 17//6 -1//12 0 0 0 0 - 2383//4100 0 0 -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0 - 3//205 0 0 0 0 -6//41 -3//205 -3//41 3//41 6//41 0 0 0 - -1777//4100 0 0 -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0], - [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0 - 0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0 41//840 41//840], - [0, 2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1] - ) - - -################################ -# Fixed step Runge-Kutta methods -################################ - -# TODO: iterator method -ode1(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_feuler) -ode2_midpoint(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_midpoint) -ode2_heun(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_heun) -ode4(fn, y0, tspan) = oderk_fixed(fn, y0, tspan, bt_rk4) - -function oderk_fixed(fn, y0, tspan, btab::TableauRKExplicit) - # Non-arrays y0 treat as scalar - fn_(t, y) = [fn(t, y[1])] - t,y = oderk_fixed(fn_, [y0], tspan, btab) - return t, vcat_nosplat(y) -end -function oderk_fixed{N,S}(fn, y0::AbstractVector, tspan, - btab_::TableauRKExplicit{N,S}) - # TODO: instead of AbstractVector use a Holy-trait - - # Needed interface: - # On components: - # On y0 container: length, deepcopy, similar, setindex! - # On time container: getindex, convert. length - - Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) - dof = length(y0) - - ys = Array(Ty, length(tspan)) - allocate!(ys, y0, dof) - ys[1] = deepcopy(y0) - - tspan = convert(Vector{Et}, tspan) - # work arrays: - ks = Array(Ty, S) - # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place - ytmp = similar(y0, Eyf, dof) - for i=1:length(tspan)-1 - dt = tspan[i+1]-tspan[i] - ys[i+1][:] = ys[i] - for s=1:S - calc_next_k!(ks, ytmp, ys[i], s, fn, tspan[i], dt, dof, btab) - for d=1:dof - ys[i+1][d] += dt * btab.b[s]*ks[s][d] - end - end - end - return tspan, ys -end - -############################## -# Adaptive Runge-Kutta methods -############################## - -ode21(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk21; kwargs...) -ode23(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk23; kwargs...) -ode45_fe(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_rk45; kwargs...) -ode45_dp(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_dopri5; kwargs...) -# Use Dormand-Prince version of ode45 by default -const ode45 = ode45_dp -ode78(fn, y0, tspan; kwargs...) = oderk_adapt(fn, y0, tspan, bt_feh78; kwargs...) - -function oderk_adapt(fn, y0, tspan, btab::TableauRKExplicit; kwords...) - # For y0 which don't support indexing. - fn_ = (t, y) -> [fn(t, y[1])] - t,y = oderk_adapt(fn_, [y0], tspan, btab; kwords...) - return t, vcat_nosplat(y) -end -function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S}; - reltol = 1.0e-5, abstol = 1.0e-8, - norm=Base.norm, - minstep=abs(tspan[end] - tspan[1])/1e18, - maxstep=abs(tspan[end] - tspan[1])/2.5, - initstep=0., - points=:all - ) - # Needed interface: - # On components: - # - note that the type of the components might change! - # On y0 container: length, similar, setindex! - # On time container: getindex, convert, length - - # For y0 which support indexing. Currently y0<:AbstractVector but - # that could be relaxed with a Holy-trait. - !isadaptive(btab_) && error("Can only use this solver with an adaptive RK Butcher table") - - Et, Eyf, Ty, btab = make_consistent_types(fn, y0, tspan, btab_) - # parameters - order = minimum(btab.order) - timeout_const = 5 # after step reduction do not increase step for - # timeout_const steps - - ## Initialization - dof = length(y0) - tspan = convert(Vector{Et}, tspan) - t = tspan[1] - tstart = tspan[1] - tend = tspan[end] - - # work arrays: - y = similar(y0, Eyf, dof) # y at time t - y[:] = y0 - ytrial = similar(y0, Eyf, dof) # trial solution at time t+dt - yerr = similar(y0, Eyf, dof) # error of trial solution - ks = Array(Ty, S) - # allocate!(ks, y0, dof) # no need to allocate as fn is not in-place - ytmp = similar(y0, Eyf, dof) - - # output ys - nsteps_fixed = length(tspan) # these are always output - ys = Array(Ty, nsteps_fixed) - allocate!(ys, y0, dof) - ys[1] = y0 - - # Option points determines where solution is returned: - if points==:all - tspan_fixed = tspan - tspan = Et[tstart] - iter_fixed = 2 # index into tspan_fixed - sizehint!(tspan, nsteps_fixed) - elseif points!=:specified - error("Unrecognized option points==$points") - end - # Time - dt, tdir, ks[1] = hinit(fn, y, tstart, tend, order, reltol, abstol) # sets ks[1]=f0 - if initstep!=0 - dt = sign(initstep)==tdir ? initstep : error("initstep has wrong sign.") - end - # Diagnostics - dts = Et[] - errs = Float64[] - steps = [0,0] # [accepted, rejected] - - ## Integration loop - islaststep = abs(t+dt-tend)<=eps(tend) ? true : false - timeout = 0 # for step-control - iter = 2 # the index into tspan and ys - while true - # do one step (assumes ks[1]==f0) - rk_embedded_step!(ytrial, yerr, ks, ytmp, y, fn, t, dt, dof, btab) - # Check error and find a new step size: - err, newdt, timeout = stepsize_hw92!(dt, tdir, y, ytrial, yerr, order, timeout, - dof, abstol, reltol, maxstep, norm) - - if err<=1.0 # accept step - # diagnostics - steps[1] +=1 - push!(dts, dt) - push!(errs, err) - - # Output: - f0 = ks[1] - f1 = isFSAL(btab) ? ks[S] : fn(t+dt, ytrial) - if points==:specified - # interpolate onto given output points - while iter-1= tdir*tend - dt = tend-t - islaststep = true # next step is the last, if it succeeds - end - elseif abs(newdt)0 no step size increase is allowed, timeout is - # decremented in here. - # - # Returns the error, newdt and the number of timeout-steps - # - # TODO: - # - allow component-wise reltol and abstol? - # - allow other norms - - # Needed interface: - # On components: isoutofdomain - # On y0 container: norm, get/setindex - - timout_after_nan = 5 - fac = [0.8, 0.9, 0.25^(1/(order+1)), 0.38^(1/(order+1))][1] - facmax = 5.0 # maximal step size increase. 1.5-5 - facmin = 1./facmax # maximal step size decrease. ? - - # in-place calculate xerr./tol - for d=1:dof - # if outside of domain (usually NaN) then make step size smaller by maximum - isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan - xerr[d] = xerr[d]/(abstol + max(norm(x0[d]), norm(xtrial[d]))*reltol) # Eq 4.10 - end - err = norm(xerr, 2) # Eq. 4.11 - newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified - if timeout>0 - newdt = min(newdt, dt) - timeout -= 1 - end - return err, tdir*newdt, timeout -end - -function calc_next_k!{Ty}(ks::Vector, ytmp::Ty, y, s, fn, t, dt, dof, btab) - # Calculates the next ks and puts it into ks[s] - # - ks and ytmp are modified inside this function. - - # Needed interface: - # On components: +, * - # On y0 container: setindex!, getindex, fn - - ytmp[:] = y - for ss=1:s-1, d=1:dof - ytmp[d] += dt * ks[ss][d] * btab.a[s,ss] - end - ks[s] = fn(t + btab.c[s]*dt, ytmp)::Ty - nothing -end - -# Helper functions: -function allocate!{T}(vec::Vector{T}, y0, dof) - # Allocates all vectors inside a Vector{Vector} using the same - # kind of container as y0 has and element type eltype(eltype(vec)). - for s=1:length(vec) - vec[s] = similar(y0, eltype(T), dof) - end -end -function index_or_push!(vec, i, val) - # Fills in the vector until there is no space, then uses push! - # instead. - if length(vec)>=i - vec[i] = val - else - push!(vec, val) - end -end -vcat_nosplat(y) = eltype(y[1])[el[1] for el in y] # Does vcat(y...) without the splatting - -function hermite_interp!(y, tquery,t,dt,y0,y1,f0,f1) - # For dense output see Hairer & Wanner p.190 using Hermite - # interpolation. Updates y in-place. - # - # f_0 = f(x_0 , y_0) , f_1 = f(x_0 + h, y_1 ) - # this is O(3). TODO for higher order. - - theta = (tquery-t)/dt - for i=1:length(y0) - y[i] = ((1-theta)*y0[i] + theta*y1[i] + theta*(theta-1) * - ((1-2*theta)*(y1[i]-y0[i]) + (theta-1)*dt*f0[i] + theta*dt*f1[i]) ) - end - nothing -end -function hermite_interp(tquery,t,dt,y0,y1,f0,f1) - # Returns the y instead of in-place - y = similar(y0) - hermite_interp!(y,tquery,t,dt,y0,y1,f0,f1) - return y -end diff --git a/src/tableaus.jl b/src/tableaus.jl new file mode 100644 index 000000000..097d59488 --- /dev/null +++ b/src/tableaus.jl @@ -0,0 +1,35 @@ +############################################################################### +## Coefficient Tableaus +############################################################################### + +# Butcher Tableaus, or more generally coefficient tables +# see Hairer & Wanner 1992, p. 134, 166 + +abstract Tableau{T<:Real} +# Name is the name of the tableau/method (a symbol) +# S is the number of stages (an int) +# T is the type of the coefficients +# +# For all types of tableaus it assumes fields: +# order::(Int...) # order of the method(s) +# +# For Runge-Kutta methods it assumes fields: +# a::Matrix{T} # SxS matrix +# b::Matrix{T} # 1 or 2 x S matrix (fixed step/ adaptive) +# c::Vector{T} # S +# +# For a tableau: +# c1 | a_11 .... a_1s +# . | a_21 . . +# . | a_31 . . +# . | .... . . +# c_s | a_s1 ....... a_ss +# -----+-------------------- +# | b_1 ... b_s this is the one used for stepping +# | b'_1 ... b'_s this is the one used for error-checking + +Base.eltype{T}(b::Tableau{T}) = T +order(b::Tableau) = b.order +# Subtypes need to define a convert method to convert to a different +# eltype with signature: +Base.convert{T<:Tableau}(::Type{T}, tab::Tableau) = error("Define convert method for concrete Tableau types") diff --git a/src/tests/integrators.jl b/src/tests/integrators.jl new file mode 100644 index 000000000..ab0e7df2f --- /dev/null +++ b/src/tests/integrators.jl @@ -0,0 +1,103 @@ +using Base.Test + +""" + test_integrator(integrator, [test_case]) + +Test a single `integrator` with either a single `test_case` or all +test cases. Test cases are defined in `ODETests.test_cases` in +`src/tests/test_cases.jl`. + +""" + +function test_integrator(integrator,test) + + ivp, sol, name, opts = test[:ivp], test[:sol], test[:name], test[:options] + + T,Y = eltype(ivp).parameters + + tol = 1//500 + @testset "$name" begin + + # 1) test the constructor + @test integrator <: AbstractIntegrator + integ=integrator(ivp;opts...) + @test typeof(integ)<:integrator + + # 2) test if the minimal backend is implemented + state=ODE.init(ivp,integ) + @test typeof(state)<:AbstractState + # output after initialization should give the initial data + @test ODE.output(state) == (ivp.t0,ivp.y0,ivp.dy0) + + # we should be able to perform the first step + @test ODE.onestep!(ivp,integ,state) == ODE.cont + # after one step the output should be updated + @test ODE.output(state) != (ivp.t0,ivp.y0,ivp.dy0) + + # 3) test the iterator interface + # pure integrator + iterator = ODE.iterate(ivp; solver=integrator, opts...) + iterator_dense = ODE.iterate(ivp; solver=ODE.DenseOutput{integrator}, opts...) + + tdir = sign(opts[:tstop]-ivp.t0) + + for iter in (iterator,iterator_dense) + for (t,y,dy) in iter + # TODO: is there a better way of doing this? We need + # a single @test statement in case of a failure for + # @testset to work properly. + if maxabs(y-sol(t)) > tol + @test maxabs(y-sol(t)) <= tol + break + end + if tdir*t > tdir*opts[:tstop] + @test tdir*t <= tdir*opts[:tstop] + break + end + end + + # generator comprehension + @test all(collect((maxabs(y-sol(t))<=tol for (t,y) in iter))) + end + + tout = opts[:tout] + @test collect((t for (t,y) in iterator_dense))==tout + + # Solution type + solution = ODE.solve(ivp; solver=integrator, opts...) + solution_dense = ODE.solve(ivp; solver=ODE.DenseOutput{integrator}, opts...) + + for s in (solution,solution_dense) + @test all(map((t,y)->maxabs(y-sol(t))<=tol,solution.t,solution.y)) + end + + # TODO: The following is not a test for a particular + # integrator but rather for the implementation of + # `ODE.solve`. It should be moved out to another test + # function. + + # If neither `tstop` nor `tout` was specified throw an error + @test_throws ErrorException ODE.solve(ivp; solver=integrator, delete!(delete!(Dict(opts),:tstop),:tout)...) + # If `tstop` was specified but is infinite + @test_throws ErrorException ODE.solve(ivp; solver=integrator, merge(Dict(opts), Dict(:tstop=>T(1//0)))...) + end +end + +function test_integrator(integrator) + @testset "$integrator" begin + for case in values(test_cases) + ODETests.test_integrator(integrator,case) + end + end +end + +# this is the functionality not yet included in test_integrator +function unused() +# 3) test the backend API + if properties != nothing + order, name, isadaptive = properties + @test ODE.order(integ) == order + @test ODE.name(integ) == name + @test ODE.isadaptive(integ) == isadaptive + end +end diff --git a/src/tests/minimal_types.jl b/src/tests/minimal_types.jl new file mode 100644 index 000000000..7f91faad4 --- /dev/null +++ b/src/tests/minimal_types.jl @@ -0,0 +1,72 @@ +import Base: <, >, >=, <= +import Base: +, -, *, /, ^ +# for the vector type +import Base: getindex, setindex!, similar +# for the scalar type +import Base: eps, convert, promote_rule, sqrt, isfinite + + +# position variable +type Position{T} <: AbstractVector{T} + x::T + y::T +end + +similar{T}(p::Position{T},::Type{T}) = Position(p.x,p.y) + +for op in (:+, :-) + @eval ($op)(p1::Position,p2::Position) = Position(($op)(p1.x,p2.x), + ($op)(p1.y,p2.y)) +end + +Base.size(::Position)=(2,) +getindex{T}(p::Position{T},i::Int) = i==1 ? p.x : p.y +setindex!{T}(p::Position{T},val::T,i::Int) = i==1 ? p.x=val : p.y=val + +# MyFloat variable. It can be constructed from any type but cannot be +# converted, so operations between MyFloat{Float64} and Float64 should +# throw a conversion error. This is to detect operations mixing high +# precision floats (like BigFloat) with lower precision constants, +# which could result in decreasing the overall precision of the +# algorithm. +immutable MyFloat <: Real + t::Float64 + # the dummy is here to prevent the use of the default constructor + dummy::Int +end + +# we need these to construct MyFloat from constants, constants are +# predefined in terms of numbers of inifinite precision such as Int or +# Rational +convert(::Type{MyFloat},s::Rational) = MyFloat(convert(Float64,s),0) +convert(::Type{MyFloat},s::Integer) = MyFloat(convert(Float64,s),0) +promote_rule{R<:Rational}(::Type{MyFloat},::Type{R}) = MyFloat +promote_rule{R<:Integer}(::Type{MyFloat},::Type{R}) = MyFloat + +eps(::Type{MyFloat}) = MyFloat(eps(Float64),0) +isfinite(x::MyFloat) = isfinite(x.t) +# necessary for the modified Rosenbrock integrator +sqrt(t::MyFloat)=MyFloat(sqrt(t.t),0) + +# binary operators +for op in (:+, :-, :*, :/, :^) + @eval ($op)(t1::MyFloat,t2::MyFloat) = MyFloat(($op)(t1.t,t2.t),0) +end + +# See #18114 +^(x::MyFloat, y::Rational) = x^(convert(MyFloat,y.num)/convert(MyFloat,y.den)) + +# unary operators +for op in (:-,) + @eval ($op)(t::MyFloat) = MyFloat(($op)(t.t),0) +end + +# comparisons +for op in (:<, :>, :>=, :<=) + @eval ($op)(t1::MyFloat,t2::MyFloat) = ($op)(t1.t,t2.t) +end + +# these are only necessary because they are used in the definition of +# the ODE test case (see test_cases.jl, :harmonic_minimal_types) +Base.sin(t::MyFloat)=MyFloat(sin(t.t),0) +Base.cos(t::MyFloat)=MyFloat(cos(t.t),0) diff --git a/src/tests/test_cases.jl b/src/tests/test_cases.jl new file mode 100644 index 000000000..5e27114a7 --- /dev/null +++ b/src/tests/test_cases.jl @@ -0,0 +1,93 @@ +# some standard test cases +const test_cases = + Dict(:constant_in_time=> + Dict(:ivp => ExplicitODE(0.0,[0.0], + (t,y,dy)->dy[:]=6.0, + J! = (t,y,dy)->dy[1]=0.0), + :sol => t->[6t], + :name => "y'=6 (vector)", + :options => Dict(:initstep => 0.1, + :tstop => 1.0, + :tout => [0.0,0.1,1.0]) + ), + + :variable_in_time=> + Dict(:ivp => ExplicitODE(0.0,[0.0], + (t,y,dy)->dy[1]=2t, + J! = (t,y,dy)->dy[1]=0.0), + :sol => t->[t^2], + :name => "y'=2t", + :options=> Dict(:tout => [0:0.001:1;], + :tstop => 1.0, + :initstep => 0.001) + ), + + :linear=> + Dict(:ivp => ExplicitODE(0.0,[1.0], + (t,y,dy)->dy[1]=y[1], + J! = (t,y,dy)->dy[1]=1.0), + :sol => t->[exp(t)], + :name => "y'=y", + :options=> Dict(:tout => [0:0.001:1;], + :tstop => 1.0, + :initstep => 0.001) + ), + + :backward_in_time=> + Dict(:ivp => ExplicitODE(1.0,[1.0], + (t,y,dy)->dy[1]=y[1], + J! = (t,y,dy)->dy[1]=1.0), + :sol => t->[exp(t-1)], + :name => "y'=y backwards", + :options=> Dict(:tout => [1:-0.001:0;], + :tstop => 0.0, + :initstep => 0.001) + ), + + :harmonic=> + Dict(:ivp => ExplicitODE(0.0,[1.0,2.0], + (t,y,dy)->(dy[1]=-y[2];dy[2]=y[1]), + J! = (t,y,dy)->copy!(dy,Float64[[0,1] [-1,0]])), + :sol => t->[cos(t)-2*sin(t), 2*cos(t)+sin(t)], + :name => "harmonic (with Jacobian)", + :options=> Dict(:tout => [0:.1:1;], + :tstop => 1.0, + :initstep => 0.001) + ), + + :harmonic_no_jac=> + Dict(:ivp => ExplicitODE(0.0,[1.0,2.0], + (t,y,dy)->(dy[1]=-y[2];dy[2]=y[1])), + :sol => t->[cos(t)-2*sin(t), 2*cos(t)+sin(t)], + :name => "harmonic (no Jacobian)", + :options=> Dict(:tout => [0:.1:1;], + :tstop => 1.0, + :initstep => 0.001) + ), + + :harmonic_minimal_types=> + Dict(:ivp => ExplicitODE(MyFloat(0), + Position(MyFloat(0),MyFloat(1)), + (t,y,dy)->(dy.x=y.y;dy.y=-y.x), + J! = (t,y,J)->(J[:]=0;J[2,1]=1;J[1,2]=-1)), + :sol => t->Position(sin(t),cos(t)), + :name => "harmonic (minimal types)", + :options => Dict(:initstep => MyFloat(1//1000), + :tstop => MyFloat(1), + :tout => [MyFloat(0),MyFloat(1//10),MyFloat(1)]) + ), + + # TODO: ForwardDiff Jacobian doesn't seem to work with custom AbstractVector type + + # :harmonic_minimal_types_no_jac=> + # Dict(:ivp => ExplicitODE(MyFloat(0.0), + # Position(MyFloat(0.0),MyFloat(1.0)), + # (t,y,dy)->(dy.x=y.y;dy.y=-y.x)), + # :sol => t->Position(sin(t),cos(t)), + # :name => "harmonic (minimal types, no Jacobian)", + # :options => Dict(:initstep => MyFloat(0.1), + # :tstop => MyFloat(1.0), + # :tout => [MyFloat(0.0),MyFloat(0.1),MyFloat(1.0)]) + # ) + + ) diff --git a/src/top-interface.jl b/src/top-interface.jl new file mode 100644 index 000000000..b72446a8d --- /dev/null +++ b/src/top-interface.jl @@ -0,0 +1,104 @@ +""" + +We assume that the initial data y0 is given at tspan[1], and that +tspan[end] is the last integration time. + +""" +function ode{T,Y,M<:AbstractSolver}(F, y0::Y, + tout::AbstractVector{T}; + solver::Type{M} = RKIntegratorAdaptive{:rk45}, + points = :all, + kargs...) + + t0 = tout[1] + + # construct a Problem + equation = explicit_ineff(t0,y0,F;kargs...) + if points == :all + prob = iterate(equation; + solver = solver, + tout = tout, + kargs...) + elseif points == :specified + prob = iterate(equation; + solver = DenseOutput{solver}, + tout = tout, + kargs...) + else + error("Unsupported points value (should be :all or :specified)") + end + + # determine if we have to unpack y + extract = Y <: Number + + to = Array(T,0) + yo = Array(Y,0) + for (t,y) in prob + push!(to,t) + push!(yo, extract ? y[1] : copy(y)) + end + + return (to,yo) +end + +""" + ODE.odeXX(F,y0,t0;kargs...) + +Solves an ODE `y'=F(t,y)` with initial conditions `y0` and `t0`. +""" + +ode23s(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = ModifiedRosenbrockIntegrator, kargs...) +ode1(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorFixed{:feuler}, kargs...) +ode2_midpoint(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorFixed{:midpoint}, kargs...) +ode2_heun(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorFixed{:heun}, kargs...) +ode4(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorFixed{:rk4}, kargs...) +ode21(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorAdaptive{:rk21}, kargs...) +ode23(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorAdaptive{:rk23}, kargs...) +ode45_fe(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorAdaptive{:rk45}, kargs...) +ode45_dp(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorAdaptive{:dopri5}, kargs...) +const ode45 = ode45_dp +ode78(F,y0,t0;kargs...) = ode_conv(F,y0,t0;solver = RKIntegratorAdaptive{:feh78}, kargs...) + + +function ode_conv{Ty,T}(F,y0::Ty,t0::AbstractVector{T};kargs...) + + if !isleaftype(T) + error("The output times have to be of a concrete type.") + elseif !(T <:AbstractFloat) + error("The time variable should be a floating point number.") + end + + if !isleaftype(Ty) & !isleaftype(eltype(Ty)) + error("The initial data has to be of a concrete type (or an array)") + end + + ode(F,y0,t0;kargs...) + +end + + + +""" + +Convert a out-of-place explicitly defined ODE function to +ExplicitODE. As the name suggests, the result is not going to be very +efficient. + +""" +function explicit_ineff{T,Y}(t0::T, y0::AbstractArray{Y}, F::Function; kargs...) + F!(t,y,dy) =copy!(dy,F(t,y)) + return ExplicitODE(t0,y0,F!; kargs...) +end + +# A temporary solution for handling scalars, should be faster then the +# previous implementation. Should be used only at the top level +# interface. This function cheats by converting scalar functions F +# and jac to vector functions F! and jac!. Still, solving this ODE +# will result in a vector of length one result, so additional external +# conversion is necessary. +function explicit_ineff{T,Y}(t0::T, y0::Y, F::Function; kargs...) + F!(t,y,dy) =(dy[1]=F(t,y[1])) + new_y0 = Array(Y,1) + new_y0[1] = y0 + return ExplicitODE(t0,new_y0,F!; kargs...) +end diff --git a/test/interface-tests.jl b/test/interface-tests.jl deleted file mode 100644 index dbb811dcc..000000000 --- a/test/interface-tests.jl +++ /dev/null @@ -1,82 +0,0 @@ -# Here are tests which test what interface the solvers require. - -################################################################################ -# This is to test a scalar-like state variable -# (due to @acroy: https://gist.github.com/acroy/28be4f2384d01f38e577) - -import Base: +, -, *, /, .+, .-, .*, ./ - -const delta0 = 0. -const V0 = 1. -const g0 = 0. - -# define custom type ... -immutable CompSol - rho::Matrix{Complex128} - x::Float64 - p::Float64 - - CompSol(r,x,p) = new(copy(r),x,p) -end - -# ... which has to support the following operations -# to work with odeX -Base.norm(y::CompSol, p::Float64) = maximum([Base.norm(y.rho, p) abs(y.x) abs(y.p)]) -Base.norm(y::CompSol) = norm(y::CompSol, 2.0) - -+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p) --(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p) -*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s) -*(s::Real, y1::CompSol) = y1*s -/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s) - -### new for PR #68 -Base.abs(y::CompSol) = norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia -Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2,2)), 0., 0.) -ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p)) - -# Because the new RK solvers wrap scalars in an array and because of -# https://github.com/JuliaLang/julia/issues/11053 these are also needed: -.+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p) -.-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p) -.*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s) -.*(s::Real, y1::CompSol) = y1*s -./(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s) - - -################################################################################ - -# define RHSs of differential equations -# delta, V and g are parameters -function rhs(t, y, delta, V, g) - H = [[-delta/2 V]; [V delta/2]] - - rho_dot = -im*H*y.rho + im*y.rho*H - x_dot = y.p - p_dot = -y.x - - return CompSol( rho_dot, x_dot, p_dot) -end - -# inital conditons -rho0 = zeros(2,2); -rho0[1,1]=1.; -y0 = CompSol(complex(rho0), 2., 1.) - -# solve ODEs -endt = 2*pi; - -t,y1 = ODE.ode45((t,y)->rhs(t, y, delta0, V0, g0), y0, [0., endt]) # used as reference -print("Testing interface for scalar-like state... ") -for solver in solvers - # these only work with some Array-like interface defined: - if solver in [ODE.ode23s, ODE.ode4s_s, ODE.ode4s_kr] - continue - end - t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, linspace(0., endt, 500)) - @test norm(y1[end]-y2[end])<0.1 -end -println("ok.") - -################################################################################ -# TODO: test a vector-like state variable, i.e. one which can be indexed. diff --git a/test/iterators.jl b/test/iterators.jl new file mode 100644 index 000000000..4dd93b837 --- /dev/null +++ b/test/iterators.jl @@ -0,0 +1,18 @@ + +const integrators = [ODE.RKIntegratorFixed{:feuler}, + ODE.RKIntegratorFixed{:midpoint}, + ODE.RKIntegratorFixed{:heun}, + ODE.RKIntegratorFixed{:rk4}, + ODE.RKIntegratorAdaptive{:rk21}, + ODE.RKIntegratorAdaptive{:rk23}, + ODE.RKIntegratorAdaptive{:rk45}, + ODE.RKIntegratorAdaptive{:dopri5}, + ODE.RKIntegratorAdaptive{:feh78}, + ODE.ModifiedRosenbrockIntegrator + ] + +@testset "Iterator interfaces" begin + for integ in integrators + ODETests.test_integrator(integ) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 9c88cd281..27be75b6f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,88 +1,8 @@ using ODE +using ODE.ODETests using Base.Test -tol = 1e-2 - -solvers = [ - ## Non-stiff - # fixed step - ODE.ode1, - ODE.ode2_midpoint, - ODE.ode2_heun, - ODE.ode4, - ODE.ode4ms, - ODE.ode5ms, - # adaptive -# ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed. - ODE.ode23, - ODE.ode45_dp, - ODE.ode45_fe, - ODE.ode78, - - ## Stiff - # fixed-step - ODE.ode4s_s, - ODE.ode4s_kr, - # adaptive - ODE.ode23s] - -for solver in solvers - println("using $solver") - # dy - # -- = 6 ==> y = 6t - # dt - t,y=solver((t,y)->6.0, 0., [0:.1:1;]) - @test maximum(abs(y-6t)) < tol - - # dy - # -- = 2t ==> y = t.^2 - # dt - t,y=solver((t,y)->2t, 0., [0:.001:1;]) - @test maximum(abs(y-t.^2)) < tol - - - # dy - # -- = y ==> y = y0*e.^t - # dt - t,y=solver((t,y)->y, 1., [0:.001:1;]) - @test maximum(abs(y-e.^t)) < tol - - t,y=solver((t,y)->y, 1., [1:-.001:0;]) - @test maximum(abs(y-e.^(t-1))) < tol - - # dv dw - # -- = -w, -- = v ==> v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t) - # dt dt - # - # y = [v, w] - t,y=solver((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;]) - ys = hcat(y...).' # convert Vector{Vector{Float}} to Matrix{Float} - @test maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol +@testset "ODE tests" begin + include("iterators.jl") + include("top-interface.jl") end - -# Test negative starting times ODE.ode23s -@assert length(ODE.ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [-5., 0])[1]) > 1 - - -# rober testcase from http://www.unige.ch/~hairer/testset/testset.html -let - println("ROBER test case") - function f(t, y) - ydot = similar(y) - ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3] - ydot[3] = 3.0e7*y[2]*y[2] - ydot[2] = -ydot[1] - ydot[3] - ydot - end - t = [0., 1e11] - t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8, - maxstep=1e11/10, minstep=1e11/1e18) - - refsol = [0.2083340149701255e-07, - 0.8333360770334713e-13, - 0.9999999791665050] # reference solution at tspan[2] - @test norm(refsol-y[end], Inf) < 2e-10 -end -include("interface-tests.jl") - -println("All looks OK") diff --git a/test/top-interface.jl b/test/top-interface.jl new file mode 100644 index 000000000..17f90f0d1 --- /dev/null +++ b/test/top-interface.jl @@ -0,0 +1,110 @@ +solvers = [ + ## Non-stiff + # fixed step + ODE.ode1, + ODE.ode2_midpoint, + ODE.ode2_heun, + ODE.ode4, + ODE.ode4ms, + ODE.ode5ms, + # adaptive + ODE.ode21, + ODE.ode23, + ODE.ode45_dp, + ODE.ode45_fe, + ODE.ode78, + + ## Stiff + # fixed-step + ODE.ode4s_s, + ODE.ode4s_kr, + # adaptive + ODE.ode23s] + +@testset "Legacy interfaces" begin + for solver in solvers + @testset "$solver" begin + tol = 1e-2 + + # dy + # -- = 6 ==> y = 6t + # dt + # we need to fix initstep for the fixed-step methods + t,y=solver((t,y)->6.0, 0., [0:.1:1;], initstep=.1) + @test maximum(abs(y-6t)) < tol + tj,yj=solver((t,y)->6.0, 0., [0:.1:1;], initstep=.1, J! = (t,y,dy)->dy[1]=0.0) + @test maximum(abs(yj-6tj)) < tol + @test norm(yj-y,Inf) y = t.^2 + # dt + t,y =solver((t,y)->2t, 0., [0:.001:1;], initstep=0.001) + @test maximum(abs(y-t.^2)) < tol + tj,yj=solver((t,y)->2t, 0., [0:.001:1;], initstep=0.001, J! = (t,y,dy)->dy[1]=0.0) + @test maximum(abs(yj-tj.^2)) < tol + @test norm(yj-y,Inf) y = y0*e.^t + # dt + t,y=solver((t,y)->y, 1., [0:.001:1;], initstep=0.001) + @test maximum(abs(y-e.^t)) < tol + tj,yj=solver((t,y)->y, 1., [0:.001:1;], initstep=0.001, J! = (t,y,dy)->dy[1]=1.0) + @test maximum(abs(yj-e.^tj)) < tol + @test norm(yj-y,Inf)y, 1., [1:-.001:0;], initstep=0.001) + @test maximum(abs(y-e.^(t-1))) < tol + tj,yj=solver((t,y)->y, 1., [1:-.001:0;], initstep=0.001, J! = (t,y,dy)->dy[1]=1.0) + @test maximum(abs(yj-e.^(tj-1))) < tol + @test norm(yj-y,Inf) v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t) + # dt dt + # + # y = [v, w] + t,y=solver((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;], initstep=0.001) + ys = hcat(y...).' # convert Vector{Vector{Float}} to Matrix{Float} + @test maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol + tj,yj=solver((t,y)->[-y[2]; y[1]], [1., 2.], [0:.001:2*pi;], initstep=0.001, J! = (t,y,dy)->copy!(dy,Float64[[0,1] [-1,0]])) + ysj = hcat(yj...).' # convert Vector{Vector{Float}} to Matrix{Float} + @test maximum(abs(ysj-[cos(tj)-2*sin(tj) 2*cos(tj)+sin(tj)])) < tol + @test norm(map(norm,yj-y),Inf)2y, 0., [0,1]) + # test typeof(y0)==Vector{Int} does not throw + @test_throws ErrorException t,y=solver((t,y)->[2y], [0], [0,1]) + # test typeof(y0)==Int does not throw + @test_throws ErrorException t,y=solver((t,y)->2y, 0, [0,1]) + # test if we can deal with a mixed case + @test_throws ErrorException t,y=solver((t,y)->2y, Number[1,1.1,BigInt(1)], Rational[0,1]) + end + end +end + +# Test negative starting times ODE.ode23s +@test length(ODE.ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [-5., 0])[1]) > 1 + +# rober testcase from http://www.unige.ch/~hairer/testset/testset.html +let + println("ROBER test case") + function f(t, y) + ydot = similar(y) + ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3] + ydot[3] = 3.0e7*y[2]*y[2] + ydot[2] = -ydot[1] - ydot[3] + ydot + end + t = [0., 1e11] + t,y = ODE.ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8, + maxstep=1e11/10, minstep=1e11/1e18) + + refsol = [0.2083340149701255e-07, + 0.8333360770334713e-13, + 0.9999999791665050] # reference solution at tspan[2] + @test norm(refsol-y[end], Inf) < 2.1e-10 +end