I’m currently doing research in gravitational physics and the problems I’m working with often end up as BVPs that are very stiff and sensitive to the initial conditions. It seems this is still a weak point for Julia’s BVP tools and I was wondering if there was any update on when the BoundaryValueDiffEq.jl will see further development?
If this has been addressed somewhere else I apologise, all I could find were a few comments that the BVP solvers will be revisited in the future.
While I’m new Julia, I’d be willing to learn and would be interested in contributing if that could be helpful here.
Yes, it’s a weak point of SciML. But we have a pretty grand plan for fixing it. Most BVP solvers specialize too much on ODE 2-point boundary problems, or lose any specialization. By mixing LinearSolve.jl, NonlinearSolve.jl, SparseDiffTools.jl, and ModelingToolkit.jl tearing you can generalize the specialization to a much larger class of problems since you then no longer require just a banded matrix structure (but can always automatically recover a banded matrix), which would make a much better MIRK method. Of course, that is a bit hard to see with the naked eye, but that’s how it’s progressing. We need some more development in those other tools before piecing it all together though, but would be interested to have any help possible.
In BifurcationKit, there are 2 ponits BVP solvers but there are a bit hidden because I focus on computing periodic orbits, they work in large dimension. THere is even the state of the art one based on orthogonal collocation. If the problem is stiff, you can use an homotopy with BifurcationKit to solve it.
I could extract a BVP solver if you want or you could change the boundary condition yourself in the functional.
Thank you for the update and giving some insights of where things are going.
This is probably my biggest complaint of other bvp solvers out there. They often become far too specialized on a specific problem or class of problem.
I only have a general understanding of MIRK methods. Do you have any references I could look over in regards to the techniques you plan to use to recover the banded matrices?
I think I will have to spend some time with SciML and Julia before I’d have anything to offer but I’d love to contribute when I do.
Specifically I’m solving for the equations of motion of modified gravity theories. The problem is defined by applying a line element of the form:
ds^2 = -H(r) dt^2 + G(r) d^2_{space}
Where H(r) and G(r) are the unknown metric functions defining how to construct spacetime in a given theory. These are the unknown functions that will appear in the field equations that we need to solve.
These functions typically diverge as they approach the horizon, as they do for a Schwarzschild metric. So our ICs are approximations of a divergence, and the numerical method has to step away from a divergence as well. The two together is what typically makes the problem sensitive to the IC and stiff.
I’ll have to take a look at Bifurcation.jl and ApproxFun.jl. I had looked at Bifurcation briefly but it seemed unrelated to what I was doing but I’ll go back and take another look. These suggestions are appreciated!