Radau with mass matrix

[griff10000]

Did you manage to get your Radau version accepted into scipy ?

[laurent90git]

It's been sitting idle a long time without any feedback from the community. I have eventually got a review a month ago or so. I still need to do some changes but it should be accepted !
Have you used it some more ?

[griff10000]

Hi, I have not been active in this area for some while, but recently mentioned your work to someone interested in python integrators. I looked at scipy and the Radau function still does not seem to implement the mass matrix method. This prompted me to enquire as to the status of your implementation. Kind regards, Graham W Griffiths www.pdecomp.net <http://www.pdecomp.net> From: laurent90git ***@***.***> Sent: 17 June 2022 22:38 To: laurent90git/DAE-Scipy ***@***.***> Cc: griff10000 ***@***.***>; Author ***@***.***> Subject: Re: [laurent90git/DAE-Scipy] Radau with mass matrix (Issue #2) It's been sitting idle a long time without any feedback from the community. I have eventually got a review a month ago or so. I still need to do some changes but it should be accepted ! Have you used it some more ? — Reply to this email directly, view it on GitHub <#2 (comment)> , or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALXDUOJH7IZHZRTOGWRLMUDVPTV2FANCNFSM5ZDIRRJA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[dullemond]

Hi Laurent,

I love your version of SciPy's Radau, thanks for providing it! I need it for a 1D time-dependent PDE project I am doing, in which I would like to impose boundary conditions and possibly internal conditions within the domain. If I am not mistaken, the only way to do this in a stable and reliable way is to replace the boundary equations with algebraic equations representing the boundary conditions. That is where the mass matrix of your version of Radau comes in handy.

While digging through your code, I noticed some small things:

Where you compute the error_norm some of your computations are overwritten by a later statement, and the scale is not yet computed:

    ZE = Z.T.dot(E) / h
    if self.mass_matrix is None:
        error = self.solve_lu(LU_real, f + ZE)
        error_norm = norm(error / scale)
    else: # see Hairer II, chapter IV.8, page 127
        error = self.solve_lu(LU_real, f + self.mass_matrix.dot(ZE))
        if self.index_algebraic_vars is not None:
            error[self.index_algebraic_vars] = 0. # ideally error*(h**index) 
            error_norm = np.linalg.norm( error/scale ) / (n - self.nvars_algebraic)** 0.5
            # we exclude the algebraic components, as they would otherwise artificially lower the error norm
            # error_norm = norm(error / scale)
        else:
            error_norm = norm(error / scale)

    scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol
    error_norm = norm(error / scale)

Shouldn't the scale = ... line be moved to the top, and the error_norm = norm(...) line at the end be removed?

The other issue is the detection of algebraic equations with a sparse matrix:

    if issparse(mass):
        self.mass_matrix = csc_matrix(mass)
        self.index_algebraic_vars = np.where(np.all(self.mass_matrix.toarray()==0,axis=1))[0] # TODO: avoid the toarray()

I think it can be done in sparse-matrix-style with:

    if issparse(mass):
        self.mass_matrix = csc_matrix(mass)
        mm = csr_matrix(self.mass_matrix)
        mm.eliminate_zeros()
        self.index_algebraic_vars = np.where(np.diff(mm.indptr)==0)[0]

Finally, wouldn't it be good to add after the line

    self.nvars_algebraic = self.index_algebraic_vars.size

the line:

    if self.nvars_algebraic==0: self.index_algebraic_vars = None

?

I very much hope that your RadauDAE will soon be incorporated as a new version of Radau in SciPy. Perhaps it would help to provide a minimally-modified version of Radau in which only the mass matrix is included? That might be easier for the community to gain confidence in. Although I do, in fact, like the other features you implemented into RadauDAE as well.

Best wishes,
Kees Dullemond