Constrained root finding with scipy
I'll go ahead and post the solution I've implemented. I wasn't able to solve the above problem for any root-finding method that numerically computes derivatives. However, I realized a much simpler solution is to use the bisection method, and then return an error of +/- infinity, as appropriate, in those regions where the error is poorly defined. This is pretty slow but seems to work every time so far, and speed is not critical in my application.
I'm trying to root find a scalar function using scipy.optimize.root_scalar. For my problem, the error function becomes poorly defined on some areas of the domain. Unfortunately, I don't know what these regions are beforehand (that's part of the purpose of this root-finding). I am able to provide a bracketing interval, but in general there will be some poorly-defined regions within that bracket.
The behavior I'd ideally like is to return NaN for the error function whenever these regions are encountered, and to have some root finding algorithm "learn" to avoid them. This doesn't work with root_scalar though, so I've been returning a fixed large error of the correct sign when these regions are encountered. Any suggestions for a less-hacky solution would be appreciated.