Add BVDAE solvers

[ErikQQY]

The BVDAE solvers are almost done now, but need further optimization in the style and improving its performance. The BVDAE solvers are similar to MIRK and FIRK methods, and their numerical workflow can be roughly divided into several steps:

(1) Build collocation equations with RK tableau. With Runge-Kutta coefficients at collocation points, build the collocation equations at each subinterval. After the collocation equations are built, we need to use AlmostBlockDiagonals.jl) to solve the built linear systems, which has an almost block diagonal form, we can just factorize(pLU) and substitute to find the solutions.

(2) Construct a nonlinear system. Construct a nonlinear system using collocation equations. In the nonlinear solving process, residual control is deployed to obtain the numerical solution. The sparsity of the nonlinear system is also configured according to what we did in MIRK and FIRK.

(3) Error estimation. While we use defect control for error control in MIRK and FIRK, BVDAE solvers utilize a common error control mechanism, which compares the two solutions obtained from different mesh(solutions from original mesh and halved mesh) and estimates the errors.

(4) Refine mesh according to the error norm and degree of equidistribution. With the error estimation we get, we compute the degree of equidistribution of the approximate integral of function and accordingly decide to refine the mesh or just halve the mesh.

(5) Proceed for next iteration or successful return.

There are still several issues that need to be fixed:

• When solving multi-points BVDAE, we need to specify these fixed points explicitly, I am thinking about adding a problem type like BVDAEProblem like what we did for DAEProblem
• There are not enough test cases for boundary value differential-algebraic equations, so need further investigation for more solid test cases
• Any suggestions would be welcomed

Also implemented an interface for COLDAE in ODEInterface.jl, since the ODEInterface.jl is not actively maintained, just put the implementation in fork: https://github.com/ErikQQY/ODEInterface.jl

Fix: #76

[github-actions]

Benchmark Results

	master	73c4161...	master/73c41613b3fb1e...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	6.27 ± 0.19 ms	6.32 ± 0.19 ms	0.991
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	2.35 ± 0.12 ms	2.38 ± 0.1 ms	0.987
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	0.827 ± 0.066 ms	0.828 ± 0.06 ms	0.998
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	0.869 ± 0.14 ms	0.878 ± 0.16 ms	0.99
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	0.986 ± 0.19 ms	1 ± 0.2 ms	0.983
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	1.91 ± 0.6 ms	1.91 ± 0.61 ms	1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	3.46 ± 0.92 ms	3.51 ± 0.9 ms	0.986
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0551 ± 0.015 s	0.0544 ± 0.016 s	1.01
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.0683 ± 0.017 s	0.066 ± 0.018 s	1.04
Simple Pendulum/IIP/Shooting(Tsit5())	0.242 ± 0.074 ms	0.245 ± 0.074 ms	0.988
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.0386 ± 0.0025 s	0.0389 ± 0.0021 s	0.994
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	11.2 ± 0.71 ms	11.2 ± 0.52 ms	0.998
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	3.22 ± 0.27 ms	3.23 ± 0.27 ms	0.997
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	3.34 ± 0.29 ms	3.34 ± 0.33 ms	0.999
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	3.44 ± 0.31 ms	3.44 ± 0.32 ms	1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.44 ± 2.6 ms	3.31 ± 2.6 ms	1.04
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	6.03 ± 4.2 ms	5.85 ± 4.1 ms	1.03
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0997 ± 0.0042 s	0.0991 ± 0.0053 s	1.01
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.121 ± 0.0093 s	0.12 ± 0.011 s	1.01
Simple Pendulum/OOP/Shooting(Tsit5())	0.594 ± 0.14 ms	0.59 ± 0.099 ms	1.01
time_to_load	15.7 ± 0.12 s	15.7 ± 0.11 s	1

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[ChrisRackauckas]

Why make each one a different type?

[ChrisRackauckas]

This is exactly the radau iia tableau?

[ErikQQY]

The tableaus come from the COLNEW paper: A New Basis Implementation for a Mixed Order Boundary value ODE solver, which I don't think is the radau iia tableau

[ChrisRackauckas]

Okay I dug into https://www.researchgate.net/publication/2435083_Collocation_Software_for_Boundary_Value_Differential-Algebraic_Equations a bit and talked with @YingboMa. It's a Gauss-Legendre quadrature https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature. Instad of calling this Ascher we should just call it a FIRK that is GaussLegendre.

Also, I don't really see the reason to be splitting the FIRK methods by different stage numbers, making that a runtime value could mean that we could adapt the order on the fly. But I think we can fix that in the future.

[ChrisRackauckas]

Are there tests with nonlinearleastsquaresproblem? I didn't see any.

[ErikQQY]

Are there tests with nonlinearleastsquaresproblem? I didn't see any.

Haven't found an appropriate test example for its case, previous work on this PR only ensures the "standard" BVDAE part is working fine, but I suspect COLDAE support over-constraints or under-constraints BVDAE

[ChrisRackauckas]

Trying to understand the general system here. I think it's fine for the Ascher (Gauss-Legendre) to be kept separate from the other FIRK because they have a lot going on with their adaptivity algorithm and such, that's fine. Maybe we keep the name Ascher because of that but in the docstring mention it's Gauss-Legendre but with Ascher's adaptivity algorithm.

But the bigger point, for BVDAEs, is there anything that limits them to the GL methods? What I ask is, any FIRK method should work for the DAE in mass matrix form right? It would be good to test that this is method independent. I'm trying to see where this mass matrix is being used and I only see https://github.com/SciML/BoundaryValueDiffEq.jl/pull/217/files#diff-f206170122c6dd64b293fa8ee49c4b402d313b6ad26dfccc7bb1ad11624f91f9R65 so I'm a bit confused.

[ErikQQY]

Maybe we keep the name Ascher because of that but in the docstring mention it's Gauss-Legendre but with Ascher's adaptivity algorithm.

Good points, will add them then.

But the bigger point, for BVDAEs, is there anything that limits them to the GL methods? What I ask is, any FIRK method should work for the DAE in mass matrix form right? It would be good to test that this is method independent. I'm trying to see where this mass matrix is being used and I only see https://github.com/SciML/BoundaryValueDiffEq.jl/pull/217/files#diff-f206170122c6dd64b293fa8ee49c4b402d313b6ad26dfccc7bb1ad11624f91f9R65 so I'm a bit confused.

Yeah, I can test if the FIRK tableaus can plug into this implementation but I suspected there would be extra work comparing that we can just replace FIRK tableaus with MIRK to get FIRK collocation methods for common BVPs. The current implementation follows the pattern of COLDAE, in which IIUC only needs to specify the differential vars, so the mass matrix here is only used for stating which equations should be the differential or algebraic equations.