# New Package PseudoArcLengthContinuation
**URL:** https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028
**Category:** Package Announcements
**Tags:** announcement
**Created:** [January 24, 2019, 10:56am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028 "2019-01-24T10:56:29Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 24, 2019, 10:56am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/1 "2019-01-24T10:56:29Z")
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](http://www.pefarrell.org/) 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](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/) 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,
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 24, 2019, 11:36am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/2 "2019-01-24T11:36:54Z")
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.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 24, 2019, 11:59am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/3 "2019-01-24T11:59:13Z")
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.
---
### Author:  [@dawbarton](https://discourse.julialang.org/u/dawbarton)
#### Post date: [January 24, 2019, 12:29pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/4 "2019-01-24T12:29:34Z")
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](https://sourceforge.net/projects/cocotools/?source=directory) (Matlab-based), since Julia seems ideal for such things but, again, time has got the better of me.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 24, 2019, 1:10pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/5 "2019-01-24T13:10:50Z")
> [@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](https://github.com/rveltz/PseudoArcLengthContinuation/blob/master/src/periodicorbit/PeriodicOrbitFD.jl#L28).
---
### Author:  [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)
#### Post date: [January 24, 2019, 1:14pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/6 "2019-01-24T13:14:41Z")
> [@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](https://github.com/rveltz/PseudoArcLengthContinuation/pull/1). It would be interesting to hear what parts of Documenters manual that was not sufficient.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 24, 2019, 2:57pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/7 "2019-01-24T14:57:19Z")
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!
---
### Author:  [@tkf](https://discourse.julialang.org/u/tkf)
#### Post date: [January 25, 2019, 4:27am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/8 "2019-01-25T04:27:09Z")
> [@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](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.
> [@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`](https://jw3126.github.io/Setfield.jl/latest/index.html#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.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 25, 2019, 8:00am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/9 "2019-01-25T08:00:17Z")
> [@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](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](https://github.com/tkf/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](https://github.com/JuliaDiffEq/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](https://github.com/tkf/Bifurcations.jl)…
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 25, 2019, 8:38am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/10 "2019-01-25T08:38:48Z")
> [@tkf](#):
>
> I suggest using [`Setfield.Lens`](https://jw3126.github.io/Setfield.jl/latest/index.html#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 🙂
---
### Author:  [@tkf](https://discourse.julialang.org/u/tkf)
#### Post date: [January 25, 2019, 9:30am UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/11 "2019-01-25T09:30:21Z")
> [@rveltz](#):
>
> As for [Bifurcations.jl](https://github.com/tkf/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.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 25, 2019, 7:35pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/12 "2019-01-25T19:35:42Z")
> [@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
```julia
[ 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.
```
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [January 25, 2019, 7:46pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/13 "2019-01-25T19:46:56Z")
```julia
│ - ✘ 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](https://github.com/rveltz/PseudoArcLengthContinuation.jl/blob/master/docs/make.jl#L19) needs to be updated to match the repo name.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [January 25, 2019, 8:30pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/14 "2019-01-25T20:30:39Z")
Thank you!
I added an [example](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/#Example-4:-nonlinear-pendulum-with-ApproxFun-1) making use of `ApproxFun` (maybe @dlfivefifty has a better way of doing this ).
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [March 21, 2019, 9:34pm UTC](https://discourse.julialang.org/t/new-package-pseudoarclengthcontinuation/20028/15 "2019-03-21T21:34:05Z")
Dear All,
I am writing this post to report a couple of improvements to [PseudoArcLengthContinuation.jl](https://github.com/rveltz/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](https://github.com/rveltz/PseudoArcLengthContinuation.jl/blob/master/src/periodicorbit/PeriodicOrbitFD.jl#L125), see [issue](https://github.com/JuliaArrays/BlockArrays.jl/issues/78) 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](https://github.com/Jutho/KrylovKit.jl/issues/15). Hence, I propose another example described in the [docs](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/#Example-5:-the-Swift-Hohenberg-equation-on-the-GPU-1) 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