Radau - missing option "mass"

[griff10000]

I am trying to use scipy.integrate.solve_ivp for a system with a mass matrix 'M', but setting method='Radau' and using **options as options={'mass':M} throws the warning:

UserWarning: The following arguments have no effect for a chosen solver: mass.

Q1: is the mass matrix option implemented in Radau?
Q2: If the mass matrix option is not implemented, is it likely to be added to Radau or any other solver in the near future?

Thanks for considering these questions.

[nmayorov]

A1: No, mass matrix option is not implemented.

A2: If someone takes the initiative to implement that.

When we considered implicit ODEs in solve_ivp we had decided, that if we try to do that we won't be able to come up with any final design and implementation, because it was already getting quite convoluted. So we decided to go with a classic explicit initial value problem formulation.

Now as solve_ivp has settled a bit, we can try to enhance Radau by adding the mass matrix. As far as I remember the mass matrix fits quite naturally into this method, so it should be doable. I'll probably take a look at how to do that.

[griff10000]

@nmayorov Many thanks for your reply.
The Fortran version of Radau5 does include a mass matrix capability so the algorithm is available.
Including a 'mass matrix' option to the scipy version of Radau would provide it with the capability to solve a class of index-1 DAE problems, which would be great.
At present, there appear to be very few python packages that offer solution capabilities for DAE problems, and those that do seem to be rather specialized.

[griff10000]

The system closed this issue and attributed the closure to me - some how I mistakenly hit the wrong button :-(.

[griff10000]

A Matlab version of RadIIA is also available at: https://uk.mathworks.com/matlabcentral/fileexchange/56162-radau-iia

[laurent90git]

Hi,
I've also been interested in this enhancement and I am currently working on it. I only consider constant mass matrices. It involves a few changes on the Newton loop, the jacobian formulation and the error estimate. So far I get good results on a simple linear vector ODEs.
I am still to test it on nonlinear ODEs, as well as index-1 DAEs.
I'll report with a more detailed explanation when I'm satisfied with the implementation.

[griff10000]

Hi, Thanks for the update regarding Radau and the mass matrix addition. Glad that you are making progress. Please continue to keep me informed. Kind regards, Graham W Griffiths

[laurent90git]

Here is a summary of my implementation so far.

I have added the mass matrix as an additional optional field "mass" in the Radau solver constructor.
If it is set to None, then the previous implementation is used. Otherwise it may be a sparse or full matrix. If it is sparse, it is converted to the CSC format.

The mass matrix appears in the Newton loop (function "solve_collocation_system"), as well as in the Jacobian of the corresponding residual function, and in the integration error estimate.

I have also added the field "solver" to the OdeResult object (in a hacky way...), so that we can access more detailed statistics after the solution has been computed (e.g. number of Newton iterations...).

To verify the implementation, I have used 2 test cases.

ODE test case

The first one is a simple vector ODE on the state vector y=(y0,y1,y2,...,yn), which reads:
    d(y[k])/dt = (k+1)*y       (1)
  <=>   dy/dt = diag([1,2,...,n+1])*y
This simple ODE (no mass matrix) can also be reformulated as a:
    M * dy/dt = y       (2)
with M=diag([1, 1/2, ..., 1/(n+1)])
I solve (1) with Radau and mass=None, and solve (2) with mass=M. The theoretical solution is also available and reads:
    y[k] (t) = y0[k] * exp( (k+1) * t )

I compared the three solutions, using tight tolerances on Radau (atol=rtol=1e-10), and they are all close up to the relative precision rtol. Radau uses the same number of steps in both cases. Thus the mass matrix implementation is coherent in the case where this matrix is non-singular.

DAE test case

The second test is based on the stiff nonlinear Robertson ODE, which reads:
dy[0]/dt = -0.04 * y[0] + 1e4 * y[1] * y[2]
dy[1]/dt = 0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * (y[1]**2)
dy[2]/dt = 3e7 * (y[1]**2)
with initial conditions y(t0)=[1,0,0]
By summing the three time derivatives, we see that the system has the following invariant:
y[0] + y[1] + y[2] = 1

This allows the original ODE (no mass matrix) to be reformulated as an index-1 semi-explicit DAE (also called Hessenberg index-1 DAE):
     dy[0]/dt = -0.04 * y[0] + 1e4 * y[1] * y[2]
     dy[1]/dt = 0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * (y[1]**2)
       0 = y[0] + y[1] + y[2] - 1
where y[2] has become a differential variable.

This DAE system is implemented by introducing a singular mass matrix M such that:
     M * Y = f
with M=[[1,0,0],
      [0,1,0],
      [0,0,0]]
and f = [ -0.04 * y[0] + 1e4 * y[1] * y[2],
      0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * (y[1]**2),
      y[0] + y[1] + y[2] - 1 ]

We can then verify that the DAE is solved correctly by comparing its solution to the one of the original ODE system.

Both equations have been solved with Radau (with mass set accordingly). Up to t~0.1, both are exactly identical in terms of solution values, number of LU decompositions, number of Newton iterations.
If the final time is pushed to 5e6, the DAE formulation struggles more to converge. I suspect it is due to the fact that the step size h becomes huge, which makes the Newton Jacobian of the system slightly worse in terms of condition number. I am still investigating this issue. This might however not be a problem for most users.

Here is a plot of the evolution of solution components throughout the integration, as well as the time step used and the cumulated number of step, with the maximum number of Newton iterations for Radau set to the default value of 6 (via the constant variable NEWTON_MAXITER):

The statistics are the following:

	time steps	function evaluations	jacobian evaluations	LU decompositions	Forward-backward substitutions
ODE	271	2029	47	256	1442
DAE	645	11616	612	2284	7958

If we set NEWTON_MAXITER to 100, the DAE solution is able to use time steps equivalent to those of the ODE solution, but of course the number of LU decompositions and forward-backward substitutions explodes.

I will continue to do some tweaks on this and will write these two test cases as unit tests for Radau. I might another test case with a singular mass matrix, without having rows full of zeros. Then I'll commit this on Scipy's fork and suggest a merge.
Note that a similar implementation may probably be performed for the BDF solver.

Regards,

[griff10000]

Hi,
@laurent90git Well done! It Looks like you are making very good progress.

[laurent90git]

Hi,

here's an update. I've made a Git with my modifications to Scipy's Radau and a few test cases:
https://github.com/laurent90git/DAE-Scipy

I have added various printouts in the code to analyse the evolution of the Newton solution, the error estimate and other parameters.

Various test problems are given, which attest that the new version of Radau is able to solve DAE up to index 3. These tests include:

• a system linear ODEs with a mass matrix
• an index-1 DAE with singular (but no nil rows) mass matrix (transistor amplifier)
• a stiff DAE of index 1 (Roberston chemical system)
• a higher-index problem (pendulum, different forumlations with index 0 to 3).

Maybe I'll add a discretised PDE next. I'll also verify the order of convergence on some of the test cases.

Following my read of "The numerical solution of differential-algebraic systems by Runge-Kutta methods" by Hairer, Lubich, Roche, I have decided to force to 0 the error estimate for the algebraic variables. These variables are at the moment identified as those that correspond to nil rows in the mass matrix. The transistor amplifier does not allow for such an identification, however it works fine without this trick. The authors of the previously mentioned book recommend multiplying the errors on the algebraic variables by the time step h (index-1), or h^2 (index 2).

Is the Scipy team willing to integrate this modification of the solver ? :)

Regards