New Package PseudoArcLengthContinuation

[rveltz]

Dear All,

I have accumulated, over the last years, some methods for performing pseudo-arclength continuation of solutions of PDE or large scale problems. I decided to package these methods and release it publicly in case others find it useful. I tried my best to write something customizable where one can easily specify a custom linear solver (GMRES, , \ on GPU…) and eigen solver to adapt the method to the specificity of the problem under study. It works reasonably well although I was not very careful about allocations (premature optimization is...)

So the package can perform continuation of solutions, detection of codim 1 bifurcation (Branch point, Fold point, Hopf point) and continue them as function of another parameter. In the context of matrix free methods, there aren’t so many codes which does this (pde2path, trilinos, ?)

I did not implement branch switching because I could not be bothered implementing the different normal forms (PR suggested ) BUT I implemented a very powerful method described by P Farrell which largely makes up for this and allows you to discover many more solutions than with branch switching.

Finally, I also provide some methods for computing periodic orbits and continuation of them. There are Matrix Free methods and one suitable for Sparse problems.

I did my best not to rely on AbstractVector so one can use Matrices as a state space, or use ApproxFun or ArrayFire…

Please have a look at PseudoArcLengthContinuation.jl and at the examples. Feel free to suggest improvements, design choices and submit PR

Feature	Matrix Free	Custom state
Newton	Y	Y
Newton + Deflation	Y	Y
Continuation (Natural, Secant, Tangent)	Y	Y
Branching point detection	Y	Y
Fold detection	Y	Y
Hopf detection	Y	Y
Fold continuation	Y	N
Hopf continuation	Y	N
Periodic Orbit Newton	Y	N
Periodic Orbit continuation	Y	N

Best regards,

[ChrisRackauckas]

This looks amazing! Can we get this setup so it can handle ODEProblem/SteadyStateProblem inputs? We can link over to this in the DiffEq docs for bifurcation analysis.

[rveltz]

Sure! You only need to provide the vector field and the jacobian (compute with AD in your package) to my package. It should be relatively easy. You can even pass “your” linearsolvers to the continuation function.

[dawbarton]

Looks fantastic! You beat me to producing a package on this (looking at my current work rate, probably by a few years…). I look forward playing around with it!

I’ve only had a brief look at your code (looks pretty clean ) but it looks like each problem (regular continuation versus fold continuation, etc) is implemented quite discretely, that is, it’s not possible to easily compose different problems. For example, to track a periodic orbit and an equilibrium state at the same time (e.g., a homoclinic problem). Is that the case?

I’ve been looking at developing something that is composable, along the lines of COCO (Matlab-based), since Julia seems ideal for such things but, again, time has got the better of me.

[rveltz]

@dawbarton

Looks fantastic! You beat me to producing a package on this (looking at my current work rate, probably by a few years…). I look forward playing around with it!

Thank you!

@dawbarton

For example, to track a periodic orbit and an equilibrium state at the same time (e.g., a homoclinic problem). Is that the case?

You are right, there is design tension for now and I hope it will be solved later. For now, I focused on something that works but is less elegant than what you suggest.

For example, I should call newton on [F, N] where N is the arclength constraint to perform continuation but I call a specific solver newtonPsArcLength for now which does not rely on newton.

It is on my todo list to provide a way to combine functionals. I thought about providing a way of appending functionals (and Jacobian ) like when you want to specify a phase condition. I have a newtonBordered function which does this but it is not pushed.

This said, in Julia it is very simple to do what you suggest. For example, I would “loosely” give to newton a function like (uper,u) -> BorderedVector(poPb(uper) , F(u)) where the problem for Periodic orbit is encoded here.

[fredrikekre]

@rveltz

I can’t manage to deploy docs, so that’s why I referenced the .md file in the link above.

See Fix documentation. by fredrikekre · Pull Request #1 · bifurcationkit/BifurcationKit.jl · GitHub. It would be interesting to hear what parts of Documenters manual that was not sufficient.

[rveltz]

I had forgotten the .jl in the name of the github rep. Not sure how it impacted the build of the docs. Thank you a lot for your fixes!

[tkf]

@rveltz

I did not implement branch switching because I could not be bothered implementing the different normal forms (PR suggested )

I wrote some codim 2 detection a while ago Modified Morris-Lecar model · Bifurcations.jl I probably should try to convert that into a PR to your code.

@rveltz

I also provide some methods for computing periodic orbits and continuation of them.

Cool! Any plan for heteroclinic orbit?

@ChrisRackauckas

handle ODEProblem/SteadyStateProblem inputs

I suggest using Setfield.Lens for specifying the bifurcation parameters. It’s a very handy way to make virtually any Julia object “differentiable” using ForwardDiff. Or maybe I can just pull out that part from my code.

[rveltz]

@tkf

I wrote some codim 2 detection a while ago https://tkf.github.io/Bifurcations.jl/latest/examples/morris_lecar/#Modified-Morris-Lecar-model-1 I probably should try to convert that into a PR to your code.

If you can, this would be great!

As for the heteroclinic orbits, I still have to improve the case of periodic orbits first. Right now, it works but it is not the best way to do this in large dimensions. Although there are some more urgent to do on my README…

Thank you for the suggestion about Lens

As for Bifurcations.jl, I am sorry but I did not noticed this before. It looks GREAT. You should make an announcement on discourse to make it known. Also, you should add a link to your docs on your README. I thought in Julia, we ``only’’ had PyDSTool.jl but we now have a serious work for ODE written in Julia.

I do think there is room for multiple packages doing the same thing. My goal is to do PDE where the problem is (partly) to perform linear solves efficiently. In ODE, you can do all sort of things that are impossible for PDEs. So we can have a specialized package for this. It would be possible to do it with my package by specialising the linear solves for example but we have Bifurcations.jl…

[ChrisRackauckas]

@tkf

I suggest using Setfield.Lens for specifying the bifurcation parameters. It’s a very handy way to make virtually any Julia object “differentiable” using ForwardDiff. Or maybe I can just pull out that part from my code.

I’m not sure what you mean here. Example?

But we should definitely get these two packages setup to take a problem type in, and flesh them out as the options in the bifurcations page. PyDSTool isn’t that great so there’s no reason it should stay front and center in our docs

[tkf]

@rveltz

As for Bifurcations.jl, I am sorry but I did not noticed this before. It looks GREAT. You should make an announcement on discourse to make it known.

Thanks! Well, it’s my fault not advertising it I was kind of hiding it since there are tons of things I wanted to fix. (And then I caught up in something else…) Anyway, I added a link to the docs in Github.

@rveltz

I do think there is room for multiple packages doing the same thing. My goal is to do PDE where the problem is (partly) to perform linear solves efficiently. In ODE, you can do all sort of things that are impossible for PDEs.

Yeah, my main usecase is for (low dimensional) ODE. I’d imagine continuation in PDE is a much harder beast.

[rveltz]

@fredrikekre

It would be interesting to hear what parts of Documenters manual that was not sufficient.

Despite my commits, it does not seem to push the new docs. The end of the build is as follows

[ Info: Populate: populating indices.
[ Info: RenderDocument: rendering document.
┌ Info: Deployment criteria:
│ - ✘ ENV["TRAVIS_REPO_SLUG"]="rveltz/PseudoArcLengthContinuation.jl" occurs in repo="github.com/rveltz/PseudoArcLengthContinuation.git"
│ - ✔ ENV["TRAVIS_PULL_REQUEST"]="false" is "false"
│ - ✔ ENV["TRAVIS_TAG"]="" is (i) empty or (ii) a valid VersionNumber
│ - ✔ ENV["TRAVIS_BRANCH"]="master" matches devbranch="master (if tag is empty)
│ - ✔ ENV["DOCUMENTER_KEY"] exists
└ Deploying: ✘
The command "julia --project=docs/ docs/make.jl" exited with 0.

[kristoffer.carlsson]

│ - ✘ ENV["TRAVIS_REPO_SLUG"]="rveltz/PseudoArcLengthContinuation.jl" occurs in repo="github.com/rveltz/PseudoArcLengthContinuation.git"

https://github.com/rveltz/PseudoArcLengthContinuation.jl/blob/master/docs/make.jl#L19 needs to be updated to match the repo name.

[rveltz]

Thank you!

I added an example making use of ApproxFun (maybe @dlfivefifty has a better way of doing this ).

[rveltz]

Dear All,

I am writing this post to report a couple of improvements to PseudoArcLengthContinuation.jl.

1. The first concerns the computation of periodic orbits for systems where the Jacobian can be expressed with a sparse matrix. The issue in the code was a slow conversion from BlockArray to a sparse matrix. This has been fixed by the function blockToSparse here, see issue and quick answer by @dlfivefifty . This allows to study much larger PDE problems and compute associated periodic orbits. The third example in the docs shows an application.

2. The second improvement allows the computation of branches entirely on the GPU using Matrix-Free methods thanks to KrylovKit.jl !!
   I have to thank @juthohaegeman for his help on this issue. Hence, I propose another example described in the docs where I show how to do this for a nontrivial problem.

I hope you will find it useful,

If you think you can improve my code, advice me… do not hesitate!

Best regards