bifurcationkit / BifurcationKit.jl

A Julia package to perform Bifurcation Analysis
https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable
Other
288 stars 34 forks source link

(potentially long-term) Complete integration with ModelingTookit's nonlinear problem #148

Open TorkelE opened 3 months ago

TorkelE commented 3 months ago

It is currently possible to define a set of non-linear equations in ModelingTookit, and use these as input to solve a bifurcation problem (https://docs.sciml.ai/ModelingToolkit/stable/tutorials/bifurcation_diagram_computation/). This depends on an extension which converts MTK's NonlinearProblem to a f and a J function that BifurcationKit can interpret internally. Would it make sense to actually store the NonlinearProblem all the way through BifurcationKit's internal process?

I could see some advantages with this. Some really mundane, like parameter's "names" being stored, so these can easier for labeling plot x-axises with parameters' names. Furthermore, if the full bifurcation path is stored, SymbolicIndexingInterface could be used to get the output path's values for a specific variable, e.g.

br = continuation(prob, ...)
br[:X] #gives the full path of `X` values.
br[:Y] #gives the full path of `Y` values.

Also, NonlinearSolve.jl has recently seen a lot of work (https://arxiv.org/abs/2403.16341). Here, BifurcationKit could directly leverage the Newton (and other nonlinear solvers) implemented in NonlinearSolve.jl (it would be even cooler if NonlinearSolve.jl could utilise features implemented in BifurcationKit only, like deflated Newton).

It might require some work, but I think it is worth a thought, I think it could turn out really nice.

rveltz commented 3 months ago

You can already have this feature br.X when X is in record_from_solution

rveltz commented 3 months ago

As for NonlinearSolve, I am thinking about this with @Datseris My main blocking point for now is that I want to allow a vector inferface which is not an AbstractArray like KrylovKit do.

TorkelE commented 3 months ago

You mean defining values in plain vector form like

u = [1.0, 2.0, 3.0]

instead of a map

u = [X => 1.0, Y => 2.0, Z => 3.0]

?

rveltz commented 3 months ago

No, define

prob = BifurcationProblem(..., record_from_solution = (x,p) -> (X = x[1], Y = x[2]))
TorkelE commented 3 months ago

Right, that makes sense. Yeah, that would require some modification (or supporting alternative approaches simultaneously, but that also seems messy).

rveltz commented 3 months ago

The biggest advantage would be to tackle parameter derivatives not with finite differences

TorkelE commented 3 months ago

Actually, I checked again and vector interface should be supported for NonlinearSystems. Presumable, in your example, x would end up a SciML integrator object. For these, both integrator[X] and integrator.u[1] should be supported (Furthermore, integrator[1] would be easy to support, if you'd need that).

rveltz commented 3 months ago

Or really?? you mean this https://jutho.github.io/KrylovKit.jl/stable/ and VectorInterface.jl ?

TorkelE commented 3 months ago

Sorry, I think I am missunderstanding something. Not familiar with those packages, so there are probably stuff going on that I don't really understand.