SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.85k stars 226 forks source link

Lagrangian Descriptors #899

Closed rmsrosa closed 2 years ago

rmsrosa commented 2 years ago

I suggested implementing this idea of Lagrangian descriptors described in Painting the Phase Portrait of a Dynamical System with the Computational Tool of Lagrangian Descriptors, by S. Wiggins & V. J. García-Garrido (AMS Notices June/July 2022) and @Datseris suggested I open an issue somewhere so we could track this.

We could implement it at JuliaDynamics/DynamicalSystems.jl or somewhere here at SciML.

rmsrosa commented 2 years ago

So, here is my initial idea:

  1. One builds a problem prob of a given type, say an ODEProblem for du/dt = f(u, p, t).
  2. Then we pass it to ldprob = LDProblem(prob, L, uu), where L = L(u, p, t) is the descriptor, which is a scalar function, e.g. L(u, p, t) = norm(f(u, p, t)), and uu is some iterator with a collection of initial conditions (e.g. an Array for a mesh in phase space or in a sub-manifold of the phase space).
  3. LDProblem uses prob.f and prob.tspan to create, via ComponentArrays, a new ODEProblem for an augmented system of the form

$$ \begin{cases} \displaystyle \frac{du}{dt} = f(u, p, t) \ \displaystyle \frac{dv}{dt} = -f(v, p, 2t_0-t) \ \displaystyle \frac{dL_f}{dt} = L(u,p,t) \ \displaystyle \frac{dL_b}{dt} = L(v,p,t) \ \end{cases} $$

  1. Notice $v$ solves the system backwards. If tspan = (t0, tf), then $v$ solves it backwards in the interval (2t0 - tf, t0), since 2t0 - tf = t0 - (tf - t0). So, we solve the system forwards and backwards at the same time.
  2. $L_f$ and $L_b$ are the forward and backward integrations of the Lagrangian descriptor.
  3. Then, solving an LDProblem works via an EnsembleProblem, where at each new solve, a new initial condition is picked.
  4. At the end of each of those solves, we only need to save the values of Lf[end] and Lb[end].
  5. We can visualize the Lagrangian descriptor using a heatmap of Lf, Lb or Lf + Lb.
ChrisRackauckas commented 2 years ago

I think this either makes sense in DynamicalSystems.jl or in its own package (which could live in SciML). I am a bit at a loss as to what other package it could fit into.

rmsrosa commented 2 years ago

Yeah, it actually makes it easier for development to start in its own package, and then we think about whether to include it in DynamicalSystems.jl or not

rmsrosa commented 2 years ago

I just to mention another and simpler way of doing it. Instead of augmenting, just solve the original system and integrate the descriptor over the solution.

Datseris commented 2 years ago

Yeah this fits well in DynamicalSystems.jl, but I'd still make it its own module/package initially that you can host in JuliaDynamics if you want. I've invited you to the org already! Simpler ways to do things are definitely preferred in my mind!

rmsrosa commented 2 years ago

Thanks @ChrisRackauckas and @Datseris! Let's make it a JuliaDynamics/LagrangianDescriptors.jl, then!

PS: The simpler versions is tempting, but this is a demanding computation since we'll repeat solving the systems so many times. We will play with both versions and see if the more involved one has a clear performance advantage or not.

rmsrosa commented 2 years ago

I started the project at JuliaDynamics/LagrangianDescriptors.jl (it has a README at least! 😆 ), so I guess I will just close this issue.