Revisiting Boundary Value Problems

[avik-pal]

With the recent developments around BoundaryValueDiffEq.jl and adjoint sensitivity computation for BVPs, there seem to be some clear shortcomings in how these problems and functions are defined.

Proposals

1. TwoPointBVPFunction: Should take 2 residuals, residuala to store the residual dependent on ua and similarly for ub

   • One thing that isn't very clear to me rn, is how do we initialize the separate residuals. Do we ask the user for len(residuala) & len(residualb)? If so, then we can simply slice the residual at a specific point and still maintain the current API.
   • Why do this? TwoPointBVPs have a specific residual structure if we know the exact bc dependence on ua and ub. So we won't have to do sparsity detection and can manually reorder the Jacobian Matrix to provide a fast path compared to standard BVPs
   • To smoothen the transition, should we check that if the bc has no 5 argument function, we construct a regular BVProblem?

2. Introduce a BVPFunction. This allows us to provide things like jac_prototype and even jac and similar stuff.

   • This is of particular interest if we want to support overconstrained and underconstrained BVPs. We can always use ODEFunction but it might be better to have a specific function for this.

3. Allow bc to be out of place (if we want to use AD through it peacefully). I see two options here:

   • If f is IIP then bc has to be IIP and similar for OOP f.
   • Allow f and bc to be independent, which means we will have to do more isinplace checks but also provides a smoother end-user experience.

Implemented Features

• Updated the TwoPointBVPFunction

cc @ChrisRackauckas @ErikQQY @EthanDecleyn

[codecov[bot]]

Codecov Report

Merging #477 (593e7f0) into master (eb92995) will decrease coverage by 57.66%. The diff coverage is 0.00%.

@@            Coverage Diff             @@
##           master    #477       +/-   ##
==========================================
- Coverage   57.65%   0.00%   -57.66%     
==========================================
  Files          50      50               
  Lines        3710    3635       -75     
==========================================
- Hits         2139       0     -2139     
- Misses       1571    3635     +2064     

Files Changed	Coverage Δ
src/problems/bvp_problems.jl	0.00% <0.00%> (-44.00%)	:arrow_down:
src/scimlfunctions.jl	0.00% <0.00%> (-61.34%)	:arrow_down:

... and 45 files with indirect coverage changes

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[ChrisRackauckas]

If f is IIP then bc has to be IIP and similar for OOP f. Allow f and bc to be independent, which means we will have to do more isinplace checks but also provides a smoother end-user experience.

The first. We require matching in all other instances (for example, SDEs).

This is of particular interest if we want to support overconstrained and underconstrained BVPs. We can always use ODEFunction but it might be better to have a specific function for this.

It's just been a waiting for a BVPFunction. Yes we should do this. The reason why it was waiting though was to confirm whether we wanted two different jac and jac_prototypes, i.e. one also for bcjac and bcjac_prototype separate from the other one. I think we want jac to be the jac of f and bcjac as the jac of bc, and then the solver will make do with that information. W and W_prototype could be reserved for potentially the "global Jac" i.e. the jacobian of the nonlinear solve, which we would probably not expect most people to touch.

TwoPointBVPFunction: Should take 2 residuals, residuala to store the residual dependent on ua and similarly for ub

Well justified.

One thing that isn't very clear to me rn, is how do we initialize the separate residuals. Do we ask the user for len(residuala) & len(residualb)? If so, then we can simply slice the residual at a specific point and still maintain the current API.

Length assumes Array. We should ask for prototypes.

To smoothen the transition, should we check that if the bc has no 5 argument function, we construct a regular BVProblem?

Yeah having a transition path like that is good. That will keep old codes working, if unoptimized (but it wasn't well-optimized before anyways so whatever).

[ErikQQY]

I want to add a few more things here:

1. I implemented a draft of BVPFunction before: Add BVPFunction, I can add more features if we need.

2. Recently I am looking at how can we supply initial guess to the construction of BVProblem and TwoPointBVProblem(e.g. we pass some solution from ODE solving as the initial guess to our BVP problem, and use these initial guess as our initial state to avoid extensive computing from scratch), and I found the current ways of passing initial guess to BVP problem constructing is difficult to use. To be specific, let's see here:

https://github.com/SciML/SciMLBase.jl/blob/cce1ce0bddbbd0a1b6f313e2bfe13bbb9598061a/src/problems/bvp_problems.jl#L117-L122

https://github.com/SciML/SciMLBase.jl/blob/cce1ce0bddbbd0a1b6f313e2bfe13bbb9598061a/src/problems/bvp_problems.jl#L153-L160

I think maybe we also need another property in BVProblem and TwoPointBVProblem to store the initial guess?

[ChrisRackauckas]

I think this looks good. Failing test though? I want this to be merged at the same time as https://github.com/SciML/SciMLBase.jl/pull/497

[avik-pal]

Give me a few hours, I am finishing up the ODEInterface integration to confirm we have everything we need. The tests here also need updating due to the breaking change

[avik-pal]

If this tests, pass we can bump it to 2.0 and then merge. Testing with 2.0 doesn't work due to version conflicts while installing StochasticDiffEq and others