Pass Jacobian function to NLsolve about steadystatediffeq.jl HOT 9 CLOSED

Yeah we'd need to modify https://github.com/SciML/SteadyStateDiffEq.jl/blob/master/src/solve.jl#L30 to handle this.

from steadystatediffeq.jl.

Yes that's what I thought. I can try a PR :)

Side-ish question: Is there a reason why this __solve method dispatches on the abstract type alg::SteadyStateDiffEqAlgorithm and then checks if typeof(alg) <: SSRootfind?

Is that better than dispatching directly on alg::SSRootfind and having a fallback on the abstract type, i.e., something like __solve(..., alg::SteadyStateDiffEqAlgorithm, ...) = error("Algorithm not recognized")?

from steadystatediffeq.jl.

I thought it would evolve to support multiple libraries. Instead, in 2020 it should just use NonlinearProblem and take in a solver for that.

from steadystatediffeq.jl.

I am confused 😅. I was talking about the dispatch on the algorithm, but you are mentioning the problem type. Does that mean that the AbstractSteadyStateProblems are going to be deprecated and you want that replaced by a AbstractNonLinearProblem?

from steadystatediffeq.jl.

No, I want SSRootfind to be a method that builds a NonlinearProblem and solves it, which would then work with all packages that can solve a NonlinearProblem.

from steadystatediffeq.jl.

Ok I think I see. I will try to write something up taking this into account.

from steadystatediffeq.jl.

Actually I'm still a bit confused. Sorry if I need much hand-holding, I'm trying! I'll split this post in two, API and Jacobian use, for clarity.

API

AFAIU, SSRootfind does not build a problem but just wraps a function that solves a problem given f, u0, and abstol:

So maybe I can try to outline an API and you can tell me where I'm confusing things 😅 I think I want to do something like

prob = SteadyStateProblem(f, u0, p) # where `f` is an `ODEfunction that may contain a jacobian`
sol = solve(prob, SSRootfind(alg))

My understanding is that you want this to be recast into

prob = NonlinearProblem(f, u0, p)
sol = solve(prob, alg)

So solve(prob, alg::SSRootfind) should be recast the SteadyStateProblem as a NonlinearProblem and do something like

function solve(prob, alg::SSRootfind)
    NLprob = NonlinearProblem(prob.f, prob.u0, prob.p)
    return solve(NLprob, alg.nlsolve)
end

Is that correct?

Use of the Jacobian

I guess that in place of

df = isnothing(f.jac) ? f : OnceDifferentiable(f.f, f.jac, ...)
u = alg.nlsolve(df, u0, abstol)

to let nlsolve know about the jacobian if it is provided, but I'm not sure... What do you think?

from steadystatediffeq.jl.

Is that correct?

Yup exactly.

we could do something like (omitting the ! for readbility)

We should expand the whole NonlinearProblem to be like ODEProblem: have a NonlinearFunction with a spot for the Jacobian and all of that. Then in this algorithm you just set all of those pieces. We'll want to at least have dispatches for nonlinear solving to Sundials (Kinsol) and NLsolve.jl, and other packages could follow after. This is a bit more work of course, but getting this part of the interface complete and inline with the differential equations will have some pretty nice advantages in the near future.

from steadystatediffeq.jl.