Dear Laurent, Thanks very much for the prompt response and detailed explanation. I was not able to access the specific pages of the book, but your explanation has clarified the problem. I also noticed the issue about the *dy_norm == 0*, that you pointed out. I also had another observation. Isn't it more stable to use a safety factor while incrementing the solution?, i.e. *y += safety * dy*, instead of *y += dy.* However, I don't have a detailed understanding of the Newton method, so I could be wrong about that. I am working on the TRBDF2 method analysed by L.F. Shampine in [1]. It is available in the MATLAB ODE suite as ode23tb. It is considered to be the "optimal" second order method due to its following properties: 1.) It is first-same-as-last, so only two implicit stages need to be computed. 2.) It is diagonally implicit, so the two stages have the same iteration matrix. 3.) It is strongly s-stable. 4.) The two stages also provide for a cubic hermite interpolant for dense output. My current state of progress is that I have managed to implement the method and it's working on simple tests. For the Newton solver, I referred to the bdf code but I'll also work on improving my own understanding. The next step is figuring out how the dense output is implemented and including it in my own code. I also need to figure out how to write the unit tests. Once this is done, I will create a pull request, hopefully sometime in the next couple of weeks. I would also like to mention here that I would be very interested in contributing more solvers to the integration suite which I worked on in my PhD. Specifically, I think the following two would be valuable additions: 1.) SDIRK methods, i.e. higher order cousins of the TRBDF2, which will all have the advantage of having the same iteration matrix. 2.) Rosenbrock methods which are linearised versions of implicit RK methods, so that only linear systems need to be solved to compute the method. Kind regards, Akshay Ranade [1] Hosea, M. E., and L. F. Shampine. "Analysis and implementation of TR-BDF2." *Applied Numerical Mathematics* 20.1-2 (1996): 21-37. link: https://www.sciencedirect.com/science/article/abs/pii/0168927495001158 On Fri, 31 Dec 2021 at 11:29, FRANCOIS Laurent < francois.laurent.90@gmail.com> wrote:

Hi,

it is the same as described by Hairer and Wanner for the Radau code (also available in Scipy) in their book "Solving ordinary differential equations II: stiff and differential-algebraic equations" (page 131), an extract of which is readable here:

https://books.google.fr/books?id=NGvrCAAAQBAJ&pg=PR8&lpg=PR8&dq=newton+stopping+criteria+hairer+wanner+radau&source=bl&ots=clGOqQfdJJ&sig=ACfU3U03pe0cePHk8k1EpVfnekAFuf8WdQ&hl=fr&sa=X&ved=2ahUKEwjej-SV8Y31AhXSxYUKHfdFDuAQ6AF6BAgPEAM#v=onepage&q=newton%20stopping%20criteria%20hairer%20wanner%20radau&f=false

Basically, the convergence rate of the Newton method can be estimated after one iteration by comparing the successive norms of the Newton increment. Assuming linear convergence (actually the quasi-Newton converges superlinearly), we can estimate how many iterations are required to meet the convergence criterion (norm of the increment < tolerance). If this is predicted to occur in a number of iterations larger than the maximum allowed one, the Newton loop is stopped. That corresponds to lines 55-57 in the "bdf.py" file (solve_bdf_system function):

* if (rate is not None and (rate >= 1 or rate ** (NEWTON_MAXITER - k) / (1 - rate) * dy_norm > tol)): break*

One cause can be that the Jacobian of the ODE function has not been updated since a few steps and may have become quite a wrong approximation of the true Jacobian. Hence, the Jacobian is updated and the Newton loop is restarted. If convergence still fails, the time step is lowered, as this usually improves convergence.

Lines 62-65 of the same file checks whether convergence is obtained:

* if (dy_norm == 0 or rate is not None and rate / (1 - rate) * dy_norm < tol): converged = True break*

i.e. either if the norm has become zero (perfect convergence, which should be impossible in practice), or if it is predicted that the next iteration will lead to an increment below the tolerance. Hence, it is possible that the Newton loop is stopped before the tolerance is actually reached. However, this tolerance is sufficiently low in practice any way, such that this does not cause any issue (in particular it does not pollute the error estimate).

Actually, taking a look at that has got me wondering whether the test *dy_norm==0* should not be replaced simply by *dy_norm<tol*. Even for linear systems, *dy_norm *should never be able to reach 0 because of machine accuracy concerns. It is therefore possible that only the second criterion ( *rate / (1 - rate) * dy_norm < tol)* ) is actually triggered. But this can only occur with a second linear LU solve per step. Thus for linear systems of equations, it seems replacing (*dynorm==0*) by ( *dy_norm<tol*) could lower the number of LU solves by a factor of two !

I hope this gives some insights :)

Out of interest, which IVP solver are you implementing ?

Regards, Laurent

<http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail> Garanti sans virus. www.avg.com <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail> <#m_-4859047978151645998_DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>

Le dim. 26 déc. 2021 à 11:02, Akshay Ranade <akshay.ranade.mate@gmail.com> a écrit :

...

Hi All,

I am working on implementing a new ivp solver. I was wondering if someone could provide some reference material to understand the Newton solver's stopping criteria in scipy/integrate/_ivp/bdf.py. Many thanks.

Kind regards, Akshay Ranade

_______________________________________________ SciPy-Dev mailing list -- scipy-dev@python.org To unsubscribe send an email to scipy-dev-leave@python.org https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: francois.laurent.90@gmail.com

_______________________________________________ SciPy-Dev mailing list -- scipy-dev@python.org To unsubscribe send an email to scipy-dev-leave@python.org https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: akshay.ranade.mate@gmail.com