Implicit Radau5DAE

[ufechner7]

Please make it possible to use Radau5DAE for implicit problems. See: https://jmodelica.org/assimulo/DAE_Radau5DAE.html
This would make it much easier to port existing code from Python to Julia.

For my problems it has a much better stability and performance than the sundials IDA solver.

In Assimulo both a pure Python version and a Fortran wrapper for Python are available which should make it easy to create a similar solution for Julia.

[ChrisRackauckas]

We have RadauIIA5 and in ODEInterfaceDiffEq we have a wrapper to the Fortran radau5 (and radau). But the original Fortran radau only does linearly implicit, i.e. allows singular mass matrices. It doesn't allow the fully implicit form. So jmodelica has either modified the solver (which, I don't quite see how to change the IRK here), or they must do an automatic conversion of the fully implicit DAE to a linearly-implicit system. @YingboMa have you looked at this?

[ufechner7]

RADAU5 implicit Runge-Kutta method of order 5 (Radau IIA) for problems of the form My'=f(x,y) with possibly singular matrix M; with dense output (collocation solution). Concerning the linear algebra routines the user has the choice to link the program either with DC_DECSOL and DECSOL or with DC_LAPACK and LAPACK and LAPACKC
From: http://www.unige.ch/~hairer/software.html

So it DOES allow a singular matrix M.

[ChrisRackauckas]

Yes sorry, that was a typo. Our RadauIIA5, radau, and radau5 methods all allow DAEs through a singular mass matrix.

[ChrisRackauckas]

https://diffeq.sciml.ai/stable/solvers/dae_solve/ has all of the DAE solvers, mentioning fully implicit and mass matrix form.

[ufechner7]

Any progress on this? What is missing?

[ChrisRackauckas]

It's missing someone to work on it. RadauIIA5 supports DAEs in mass matrix form and that can be used today. A separate form would require a PR.

[ufechner7]

If I would like to implement it, which files would have to be modified?

[ChrisRackauckas]

It would just follow the standards of the OrdinaryDiffEq.jl repository. This is documented in the developer documentation http://devdocs.sciml.ai/latest/. http://devdocs.sciml.ai/latest/contributing/adding_algorithms/ describes how new algorithms are added to OrdinaryDiffEq.jl with an example on explicit RK methods. Some of the fine details of Radau and DAEs are not documented there, but that should get you started on what's important (defining a cache and a perform_step! function). I think the right way to do this might be to add a new DRadauIIA5 that uses the same cache and perform_step! dispatches as RadauIIA5 but hooks into the f calls in the right way with a if prob isa DAEProblem to swap over to the other f call etc.