Maybe there's a way to re-imagine them (using symbols or a series of types/value types?), but we need some setup of retcodes for the solution types.

This means that parameter estimation won't pass on keyword arguments that are in the solver. This isn't good.

We need to create linear operator types which will act like matrices M*u but also have a call M(t,u)=M*u so that way they can be used in place of functions. The general interface should be comparable to LinearMap. Concrete types would be akin to one like a FunctionMap, but then also ones that can represent things like the operator for the upwind-scheme, or non-constant mass matrices.

I always use Juno which has special rendering and looks great. But trying out different versions of Julia, I saw the REPL...

julia> using DifferentialEquations

julia> f(t,u) = 1.01*u
f (generic function with 1 method)

julia> u0=1/2
0.5

julia> tspan = (0.0,1.0)
(0.0, 1.0)

julia> prob = ODEProblem(f,u0,tspan)
DiffEqBase.ODEProblem{Float64,Float64,false,#f,Void,UniformScaling{Int64}}(f, 0.5, (0.0, 1.0), nothing, UniformScaling{Int64}
1*I)

julia> sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
DiffEqBase.ODESolution{Float64,1,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Float64,Float64,false,#f,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}([0.5, 0.506305, 0.52193, 0.543053, 0.569506, 0.602174, 0.641202, 0.687147, 0.740325, 0.801221, 0.870275, 0.948019, 1.03501, 1.1319, 1.23937, 1.3582, 1.3728], nothing, nothing, [0.0, 0.0124078, 0.0425009, 0.0817804, 0.128873, 0.184097, 0.246274, 0.314792, 0.388595, 0.46686, 0.548714, 0.633432, 0.720359, 0.808954, 0.898761, 0.989412, 1.0], Array{Float64,1}[[Inf], [0.505, 0.506019, 0.507074, 0.510728, 0.51124, 0.511368, 0.511368], [0.511368, 0.513871, 0.516476, 0.52555, 0.526831, 0.527151, 0.52715], [0.52715, 0.530517, 0.534033, 0.546313, 0.548052, 0.548487, 0.548483], [0.548483, 0.552683, 0.557081, 0.572476, 0.574663, 0.57521, 0.575201], [0.575201, 0.580367, 0.585789, 0.60482, 0.607532, 0.608211, 0.608196], [0.608196, 0.614345, 0.620814, 0.643571, 0.646824, 0.647638, 0.647614], [0.647614, 0.65483, 0.662437, 0.689249, 0.693091, 0.694054, 0.694018], [0.694018, 0.702347, 0.711143, 0.742198, 0.746659, 0.747777, 0.747728], [0.747728, 0.757244, 0.767308, 0.802891, 0.808014, 0.809296, 0.809233], [0.809233, 0.820004, 0.831408, 0.871778, 0.877599, 0.879057, 0.878977], [0.878977, 0.891086, 0.903919, 0.949387, 0.955951, 0.957596, 0.957499], [0.957499, 0.971034, 0.985387, 1.03628, 1.04364, 1.04548, 1.04536], [1.04536, 1.06042, 1.0764, 1.1331, 1.1413, 1.14335, 1.14322], [1.14322, 1.15991, 1.17763, 1.24053, 1.24963, 1.25191, 1.25176], [1.25176, 1.27021, 1.2898, 1.35936, 1.36943, 1.37195, 1.37178], [1.37178, 1.37414, 1.37659, 1.38505, 1.38623, 1.38653, 1.38653]], DiffEqBase.ODEProblem{Float64,Float64,false,#f,Void,UniformScaling{Int64}}(f, 0.5, (0.0, 1.0), nothing, UniformScaling{Int64}
1*I), OrdinaryDiffEq.Tsit5(), OrdinaryDiffEq.InterpolationData, true, 0, :Success)

julia> using RecursiveArrayTools

julia> VectorOfArray(sol.u)
17-element Array{Float64,1}:
 0.5
 0.506305
 0.52193
 0.543053
 0.569506
 0.602174
 0.641202
 0.687147
 0.740325
 0.801221
 0.870275
 0.948019
 1.03501
 1.1319
 1.23937
 1.3582
 1.3728

julia> DiffEqArray(sol.t,sol.u)
17-element Array{Float64,1}:
 0.0
 0.0124078
 0.0425009
 0.0817804
 0.128873
 0.184097
 0.246274
 0.314792
 0.388595
 0.46686
 0.548714
 0.633432
 0.720359
 0.808954
 0.898761
 0.989412
 1.0

I think we can start with the DiffEqArray in RecursiveArrayTools since it matches the solution object quite well. It has "time" and a solution array over time: can we get that to show well?

DiffEqBase.jl is well covered by the other libraries, but it would be good to increase the coverage here for things like problem creation.

I was trying to use a sparse u (and du) but the solver fails due to a lack of recursivecopy! method for sparse vectors. Is there an issue with supporting sparse objects?

get_du was exported so get_du! should also be exported.

If xscale is set to e.g. :log10, logspace should be used instead of
linspace. This is necessary for visualizing the example from the DAE tutorial
well (SciML/DiffEqDocs.jl#113).

So that way the users can call it when they want after discontinuities or during initialization

I've encountered some strange error messages if tType == Any in ODEProblem:

julia> using DifferentialEquations
julia> f(t,u) = u
julia> prob1 = ODEProblem(f, u0, (0.,1))
julia> solve(prob1)
ERROR: MethodError: no method matching OrdinaryDiffEq.DEOptions{uEltype,uEltypeNoUnits,tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsType,ECType,SType}(::Int64, ::Int64, ::Bool, ::Bool, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Rational{Int64}, ::OrdinaryDiffEq.#ODE_DEFAULT_NORM, ::Void, ::DataStructures.BinaryHeap{Any,DataStructures.LessThan}, ::DataStructures.BinaryHeap{Any,DataStructures.LessThan}, ::DataStructures.BinaryHeap{Any,DataStructures.LessThan}, ::Void, ::Bool, ::Int64, ::String, ::OrdinaryDiffEq.#ODE_DEFAULT_PROG_MESSAGE, ::Bool, ::Bool, ::Float64, ::Float64, ::Float64, ::Bool, ::Bool, ::DiffEqBase.CallbackSet{Tuple{},Tuple{}}, ::OrdinaryDiffEq.#ODE_DEFAULT_ISOUTOFDOMAIN, ::OrdinaryDiffEq.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool)
Closest candidates are: ...

julia> using DifferentialEquations
julia> f(t,u) = u
julia> prob2 = ODEProblem(f, u0, (0,1.))
julia> solve(prob2)
ERROR: InexactError()
 in convert at ./rational.jl:68 [inlined]
 in Int64(::Rational{Int64}) at ./sysimg.jl:53
 in #init#45(::Int64, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Int64, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Rational{Int64}, ::OrdinaryDiffEq.#ODE_DEFAULT_NORM, ::OrdinaryDiffEq.#ODE_DEFAULT_ISOUTOFDOMAIN, ::OrdinaryDiffEq.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::OrdinaryDiffEq.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::DiffEqBase.CallbackSet{Tuple{},Tuple{}}, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Float64,Any,false,#f,DiffEqBase.CallbackSet{Tuple{},Tuple{}}}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ...

Since these errors occur not in the construction of the ODEProblem but in solve, it took me some time to fix it. Thus, it might make sense to test whether tType == Any in ODEProblem. I could start implement such tests in the constructors in ode_problems.jl.

We might at least want to print some warnings if tType == Any.

I think it would be reasonable to document an advanced timestepping method for fixed timestep methods as follows:

• If dt is chosen, every step will be of size dt unless in the interval there is a value of tstops, which it will then walk to that value
• If dt is not chosen and tstops is empty, automatic choice of dt will occur
• If dt is not chosen and tstops is not empty, the algorithm will walk to each value of tstops

@mauro3 does that make sense? This would use the tstops vector instead of special casing dt into a vector for this case.

Let's start with a brief overview. We started extending ODEProblem to things like SplitODEProblem by giving every single one of them a different type. They had exactly the same type structure as ODEProblem, but there were tons of them. This made maintenance a mess, because then a change to noise would need to change 10 different SDEProblem extras, all of which are quite niche.

The solution there was to instead just have everything be an ODEProblem. ODEProblem((f1,f2),u0,tspan) means SplitODEProblem, and this is interpreted by the solver. These definitions are implicit in the setup and are only enforced by the solver algorithm: i.e. it will error if it's wrong, and that's about it. This works okay, but it makes conversions between different "subproblem types" neigh impossible: if you want to convert that into a standard ODEProblem via adding the functions, it's not easy to write or maintain code to do that given that ODEProblem((f1,f2),u0,tspan) can mean different things. It also then becomes hard to document what the "Refined ODE Problems" mean.

Fast forward to today with discussions for JuliaLang/julia#23089 . With the change of ODEProblem to allow for the user to set if it's inplace via ODEProblem{iip}(...), I wanted to replicate that in the SecondOrderODEProblem "helper function". For reference:

struct SecondOrderODEProblem{iip} end
function SecondOrderODEProblem(f,u0,du0,tspan;kwargs...)
  iip = isinplace(f,4)
  SecondOrderODEProblem{iip}(_f,_u0,tspan;kwargs...)
end
function SecondOrderODEProblem{iip}(f,u0,du0,tspan;kwargs...)
  if iip
    f1 = function (t,u,v,du)
      du .= v
    end
  else
    f1 = function (t,u,v)
      v
    end
  end
  _f = (f1,f)
  _u0 = (u0,du0)
  ODEProblem{iip}(_f,_u0,tspan;kwargs...)
end

This actually isn't its own type, and just builds an ODEProblem for maintainability (in partitioned form, to be exact). But the interesting thing here is that we have struct SecondOrderODEProblem{iip} end which is just thrown away.

What about using that kind of struct as a trait on the ODEProblem for designating different problem subtypes? These traits can then be applied to things like @YingboMa's recent work in DiffEqPhysics, and a trait hierarchy can be built. All of these will use the same underlying ODEProblem, so issues with handling mass matrices and other things like that will still be uniform across the types automatically, keeping the maintenance cost down. This will give to the user a constructor function that they need to use, and the solvers can check/convert to different forms.

The idea is that it is often useful to get MonteCarlo intermediate results in a form of statistic summary MonreCarloSummary. I think it could be implemented in two ways (let's say we need to obtain MonreCarloSummary each n iterations of num_monte):

1. Via reduction function which allows to obtain (print, save, upload to a DB) statistic parameters like MonreCarloSummary each n iterations but doesn't reduce the solution data as it will be needed for future MonreCarloSummary calculations
2. Via reduction function which calculates statistic parameters like MonreCarloSummary and reduces (deletes) all the solution data each n iterations. Here as we delete the solution data and save only MonreCarloSummary the reduction function should use some "online" algorithms to calculate and approximate statistic parameters each n iterations using, for example, OnlineStats.jl

Would it be possible to opt-out from the update warning once I have seen it?

I have to admit that I don't have a suggestion for the implementation though. Eg by setting an environment variable?

Personally, I like it when the axes' limits do not change during an animation. Would you be interested in applying default xlims and ylims in the @recipe to make animate have only the series changing, but not the limits of the axes? What do you think? Mine is of course just a taste...

The tag name "0.11.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @ChrisRackauckas

It should just be a parameterically-typed field where we check for it being a test if it's not void. It's silly to double up the entire infrastructure when this solution is just as efficient.

At the moment tuples(sol) returns tuples of the form (u,t) for ODEProblems. In contrast, in the case of a MonteCarloSimulation tuples(timeseries_point_mean(sim,ts)) returns tuples of the form (t,u).

Continuing from #1 , we have the types, we just need conversions and promotions. The conversions are mostly trivial, except @mauro3 found that handling the Jacobians is non-trivial. We should work to get everything except the Jacobians working first, then go back to working on those.

As discussed in SciML/Roadmap#5, we have set out for the DAE types something like DAEProblem for ODEs defined by an implicit function f(t,u,du,resid) where resid is the residual (does anyone have a better common name for that arg? Or is that good?), and MassMatrixODEProblem for problems with a mass matrix (extends the ODE problem to have a spot for the mass matrix). What should be name that field? mass_matrix? We would then implement a promotion structure:

• ODEProblem -> MassMatrixODEProblem by making the mass_matrix = I the uniform scaling operator (allows one to do a type check as well to turn off extra calculations, and acts nicely in the common ways mass matrices are used).
• MassMatrixODEProblem -> DAEProblem by defining the residual problem out = M*du - f(t,u,du) (as real code, of course this won't work)? In this promotion we also need to provide du0 and that can be determined by an evaluation of f.
• Is there an intermediate common case we are missing?

Ideally you'd then just write your algorithm at the top level with enough promotion information exposed so that way you don't lose performance. So say for my Rosenbrock algorithms it can solve a MassMatrixODEProblem, so what I write those in a manner such that, if someone does

solve(prob,Rosenbrock23;kwargs)

and prob is an ODEProblem, it will be promoted to a MassMatrixODEProblem with mass_matrix = I and given to my Rosenbrock routine.

I don't know how to do promotions really at all, so could someone explain this @mauro3?

Note that some starter code is in AlgebraicDiffEq.jl, but this is probably easier to just re-write given what we have from ODEProblem in DiffEqBase since a lot has changed. For each of these we should also have an appropriate TestProblem and build_..._problem function for testing.

I can do this over the weekend if needed, but maybe someone else who's more interested in DAEs should take a stab at this to get to know the infrastructure well and identify potential problems. @pwl or @mauro3?

Once this is together, I'll make a common interface dispatch for Sundials' idasol to move that code out of AlgebraicDiffEq.jl. One that's done, AlgebraicDiffEq.jl will be a shell package, which is fine because it will be ready for people to add methods if they so please (I'll be leaving my Rosenbrock functions and creating new ones in OrdinaryDiffEq.jl since I tend to consider mass matrix problems more ODEs. That can be a discussion though).

Was running some benchmarks again and noticed some regressions.

Before elementwise tolerances:

After elementwise tolerances:

While it seems odd to get rid of ! on a mutating function, @ccontrer 's example brings up a big issue:

callback = ContinuousCallback(condtion,affect_bounce!,rootfind,save_positions,affect_neg! = affect_negative!)

That works.

callback = ContinuousCallback(condtion,affect_bounce!,rootfind,save_positions,affect_neg!=affect_negative!)

That errors. See the difference? Julia parses affect_neg!=affect_negative! as affect_neg != affect_negative!, and errors. That's not nice behavior, and so I think it might be best to just take off the !.

It would be nice to have a DEStats type which stores useful data about how well the method worked. It should count things like:

• Number of function calls
• Number of Jacobian calculations
• Number of linear solves
• Number of accepted steps
• Number of rejected steps
• Number of interpolated points

Does broadcasting a DataArray return a DataArray? I think in the long run we should do something like

https://github.com/JuliaDiffEq/MultiScaleArrays.jl/blob/master/src/math.jl#L9

It would be good to be able to define an "initialization" function on a callback which takes in an integrator. For example, sometimes you need to set some values in the callbacks to match or use the first timepoint in the integrator. Right now this can be done hackily in the user script, but decreases the generality of the callback.

It will be required for expanding jump callbacks to the signature (t,u,integrator) since at the start the integrator doesn't exist, so its setup is currently limited to (t,u). AutoAbstol also needs an initialization with the initial maximum of u.

Instead of having isinplace as a keyword argument which then is not able to be inferred (even after kwarg dispatch is a thing), instead it should be a two-stage process, so that way an advanced user might want to

ODEProblem{true}(f,u0,tspan)

to skip the auto-detection.

When simulating systems, one often needs to keep track of the input and onput signals. Assume we have the following representation:

xdot = A*x + B*u(t,x)
y = C*x + D*u(t,x)

Then, a solution similar to

type TimeResponse{T} <: DEDataArray{T}
  x::Array{T,1}
  u::Array{T,1}
  y::Array{T,1}
end

function f(t::Real,x::TimeResponse,dx::AbstractVector,u)
  ucalc = u(t, x.x) # callable object on (::Real,::AbstractVector)
  isa(ucalc, AbstractVector) || DomainError()
  x.u = ucalc
  dx[:] = A*x.x + B*x.u # for some A, B
  x.y = C*x.x + D*x.u
end

would handle everything. However, in this way, interpolation on u and y are note available in the current design.

Moreover, u object above can be a discrete-time system object on its own, with some samplingtime(u) defined. Then, this object will calculate the new values at each samplingtime(u) and hold the previous value at all times. See, e.g.,

function (u::System)(t::Real, x::AbstractVector)
  if t >= u.tprev + samplingtime(u)
    u.tprev = t
    u.uprev = ... # calculate new value
  end
  u.uprev
end

this approach can also run easily using tstops keyword.

However, if there is gonna be some change in DEDataArray type to support the algebraic equations as those above, one needs to think about this issue. Should it be DEDataArrays responsibility to take this discrete behaviour in mind (interpolation would yield the stairs), or should it be the object's job to cache some values and not do any interpolation when asked by the DifferentialEquations ecosystem packages? Would there exist some common interface between the two?

Cc: @neveritt

Instead of using trait dispatch which is inherently unsound after 265, we should instead create a function-holder type like Optim does. This type would just hold a bunch of functions, and then have the dispatches which we currently look for in order to help with version compatibility. To help users, we'll auto-convert the dispatch-based form and throw depwarns for a bit when they put a function with dispatches into the ODEProblem etc.

I think isinplace should be a trait and not a type parameter. Thoughts?

I have a function with explicitly stated constrained generic arguments:

function equation!{T<:Real, U <:AbstractUnsaturatedModel{T}, V <: Inlet{T}}(t :: T, s :: Vector{T}, d :: Darcy{T,U,V}, Ds :: Vector{T})

Creating a parameterized function as

ParameterizedFunction(equation!, d)

Results in error in numargs (utils.jl) function that seems unable to calculate correctly the number of arguments of such a generic function. The problem is in the call

m.sig.parameters

where m is a method of equation!: there is no field parameters for my function.

For something simpler it works:

function lz(t, u, du)
	σ = 10.0
	ρ = 28.0
	β = 8/3.0
	du[1] = σ*(u[2]-u[1])
	du[2] = u[1]*(ρ-u[3]) - u[2]
	du[3] = u[1]*u[2] - β*u[3]
end

first(methods(lz)).sig.parameters

svec(DarcyExample.#lz, Any, Any, Any)

However the change to explicit argument type results in error:

function lz1{T<:Real}(t::T, u::Vector{T}, du::Vector{T})
	σ = 10.0
	ρ = 28.0
	β = 8/3.0
	du[1] = σ*(u[2]-u[1])
	du[2] = u[1]*(ρ-u[3]) - u[2]
	du[3] = u[1]*u[2] - β*u[3]
end

first(methods(lz1)).sig.parameters

type UnionAll has no field parameters

@mauro3 brought up a good point that the scope of DiffEqBase may not include PDEs. I think this is a good point because:

1. The interface on ODEs (ODE, SDE, DAE, DDE) is pretty clear given what we laid out. Sure there are a few other related things (IMEX), but if we stick to ODEs, DiffEqBase could be something which is quite stable in no time.
2. The set of PDEs will be constantly changing and growing. This is by design. The hope is that people begin to start contributing PDE problems and algorithms for solving the PDEs, and then we have a large library of readily available PDEs with solvers. This is very "batteries included" and antithetical to the purpose of DiffEqBase.
3. Other sets of tools make sense in a Base library for PDEs; meshing tools is the big one.

These facts seem to point that PDEs should be put in a separate DiffEqPDEBase (name pending), which would likely be more lively and user-focused, i.e. developed for giving users "the common interface" on tons of problems, and not for giving developers a common interface on the basic problems / types like DiffEqBase.

You are overloading display in order to customize output. This is the wrong way to do it — see the manual on custom pretty printing.

To customize REPL output, you should override Base.show(io::IO, x::MyType) and/or Base.show(io::IO, ::MIME"text/plain", x::MyType), as explained in the manual.

Overriding display directly will break IJulia, Juno, and other environments that have a custom display mechanism.

This came up over here:

SciML/OrdinaryDiffEq.jl#28

There are a lot of times where a user may want to plot not just the values, but a function of the values. Complex numbers are a good example because Plots.jl can't handle complex numbers, and so one thing that we could do is allow the user to choose a function to be applied to the number they wish to plot. This would allow someone to plot the norms of the complex numbers, or the real/imaginary parts.

I think the DiscreteProblem is not tracking the user defined types properly.

Minimal working example below:

A = diagm([0.3,0.6,0.9])
B = [1 2 3].'
C = [1/3 1/3 1/3]

type SimType{T} <: DEDataArray{T}
    x::Vector{T}
    y::Vector{T}
    u::Vector{T}
end

function mysystem(t,x,dx,u)
    ucalc = u(t,x)
    x.u = ucalc
    x.y = C*x.x
    dx[:] = A*x.x + B*x.u
end

input = (t,x)->(1*one(t)≤t≤2*one(t)?[one(t)]:[zero(t)])
prob = DiscreteProblem((t,x,dx)->mysystem(t,x,dx,input), SimType(zeros(3), zeros(1), zeros(1)), (0.,4.))
sln = solve(prob, FunctionMap(scale_by_time=false), dt = 0.1)

u1 = [sln[idx].u for idx in 1:length(sln)]
u2 = [sln(t).u for t in linspace(0,4,41)]
any(x->x[1]>0, u1)
any(x->x[1]>0, u2)

Neither getindex(sln::SolutionType, idx::Int) nor (sln::SolutionType)(t::Real) seem to be doing what they are supposed to be, if I am not doing any other stupid mistake in calling the functions or defining my objects :( The states (us) seem to be working properly.

Will do with the next Juno.jl tag for the defaultrepr function.

The tag name "1.0.1" is not of the appropriate SemVer form (vX.Y.Z).
cc: @ChrisRackauckas

When I try to plot a solution of an ODEProblem I do get the following error

UndefVarError: tmp_lab not defined

 in add_labels!(::Array{String,1}, ::Tuple{DiffEqBase.##118#119,Int64,Int64}, ::Int64, ::DiffEqBase.ODESolution{Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,false,#f},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}) at C:\Users\thiem\.julia\v0.5\DiffEqBase\src\solutions\solution_interface.jl:231
 in solplot_vecs_and_labels(::Int64, ::Array{Tuple,1}, ::Array{Array{Float64,2},1}, ::Array{Float64,1}, ::DiffEqBase.ODESolution{Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,false,#f},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}, ::Bool) at C:\Users\thiem\.julia\v0.5\DiffEqBase\src\solutions\solution_interface.jl:295
 in macro expansion at C:\Users\thiem\.julia\v0.5\DiffEqBase\src\solutions\solution_interface.jl:84 [inlined]
 in apply_recipe(::Dict{Symbol,Any}, ::DiffEqBase.ODESolution{Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,false,#f},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}) at C:\Users\thiem\.julia\v0.5\RecipesBase\src\RecipesBase.jl:238
 in _process_userrecipes(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODESolution{Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,false,#f},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}) at C:\Users\thiem\.julia\v0.5\Plots\src\pipeline.jl:73
 in _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODESolution{Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,false,#f},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}) at C:\Users\thiem\.julia\v0.5\Plots\src\plot.jl:171
 in plot(::DiffEqBase.ODESolution{Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},DiffEqBase.ODEProblem{Array{Float64,2},Float64,false,#f},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}) at C:\Users\thiem\.julia\v0.5\Plots\src\plot.jl:46

The example code I am using is the following:

using DifferentialEquations,Plots
A = [1. 0 0 -5
     4 -2 4 -3
     -4 0 0 1
     5 -2 2 3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(t,u) = A * u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
plot(sol)

I find it very restrictive to pass only u[1] and u[end] to the residual function of TwoPointBVProblem therefore one cannot access the tspan.

One possibility is to pass solution types that only hold u[1] and u[end] as u. The other possibility I see is to explicitly pass tspan.

I am finding that I need information about what parts of the common API are implemented for different algorithms. For example, in DiffEqParamEstim.jl I want to be able to generate the numerical solution at specific times. For a given algorithm should I use saveat or tstops? LSODAAlg() from LSODA.jl does not implement tstops due to issues with the dependency, while some algorithms without interpolations do not implement saveat. It would be nice to be able to query a trait to:

1. Use saveat if possible
2. Else fallback to tstops (since it's usually slower)
3. Else throw an error: there must be some way to force the solver to hit the point t!

I assume that traits for other things like progressbars would be nice because then it could throw informative warnings / errors. The traits could just default to "false" meaning "not implemented", and then manually set to true what's implemented.

If problem types were immutable, it would be free to make them. Is there ever a good reason to mutate them? Adding new callbacks or anything like that isn't possible because that changes the typing and it's strictly typed. If the u0 is mutable you can already just mutate that, but in many cases you want to change the type so instead creating a new object is really what's the right way to do it. tspan is a tuple so it can't be modified, and neither can the function. So I am not sure if the mutability is even useful.

Does anyone use the mutability of the problem type? @Datseris @devmotion @ranocha

If the units-test do their work properly, then we can remove Parameter.jl and Ranges.jl. I don't know how easily RecursiveArrayTools.jl could be removed. We're stuck with SimpleTraits.jl and RecipesBase.jl. Although the latter could be moved to a plots-package (but I'd rather not).

Ref: luchr/ODEInterface.jl#9

As we use function overloading for many extra functions, why not use it for the analytic solution as well? (Now that we're using traits dispatch would work too.) Would seem more consistent to me.

It seems that the interface for the other solvers is truly defined by the tuple length for the functions and the tuple length for the inputs. In this case, it's redundant to have separate type names. Instead, the refined ODE/SDE/DAE is just an interface on the tuple information in the standard type. I am thinking we can remove all of the extra names before they become standardized, and just write the interface via tuples.

@honghuzi's example pointed out that simple progress monitoring on the user side will not work in all cases because IO is not threadsafe.

SciML/DifferentialEquations.jl#190

We should setup a monte_progress=true at this higher level (since progress=true is for individual DE progressbars) which is threadsafe for monitoring the progress through batches.

It would be nice to have a type for 2nd order ODEs

u'' = f(t,y)

for things like Runge-Kutta-Nystrom methods. I know there are a bunch of methods which are specifically designed for these problems. Is there a way to also automatically have them convert to a 1st order ODE?

reinit would set the problem back to the initial condition and tspan[1]. This would be a great way to cheaply resolve with different parameters. Also, since one can change integrator.prob.u0 and integrator.prob.tspan, this would be a great way to solve with slightly different setups but reuse the internal caches.

Description

When using plotting functionality with vars keyword, the recipe recognizes the use of the (independent) time variable t and labels the x-axis properly. However, when the place of the time variable is switched, the recipe cannot recognize it anymore.

Not sure if such a behaviour should be needed, but there is this minor inconsistency.

Actual Behaviour

using DiffEqBase, OrdinaryDiffEq
using Plots
gr()

f(t,x,dx) = (dx[1] = -x[1]; nothing)
prob = ODEProblem(f, [1.], (0.,10.))
sln = solve(prob, Tsit5())

plot(sln, vars=[(0,1)]) # plots x[1] vs t, having x-label as `t`
plot(sln, vars=[(1,0)]) # changes the axes of above; but no `y-label` as `t`

Desired Behaviour

using DiffEqBase, OrdinaryDiffEq
using Plots
gr()

f(t,x,dx) = (dx[1] = -x[1]; nothing)
prob = ODEProblem(f, [1.], (0.,10.))
sln = solve(prob, Tsit5())

plot(sln, vars=[(0,1)]) # plots x[1] vs t, having x-label as `t`
plot(sln, vars=[(1,0)]) # changes the axes of above; AND `y-label` as `t`

Hi! Is it possible to do this? I am guessing that it would mess-up the method that "automatically detects a solver".