roottools

This module specialises class Submodel for root-finding problems.

Module roottools provides Python interfaces only. See the Cython extension module rootutils for the actual implementations of the mathematical algorithms.

Module roottools implements the following members:

• Pegasus Root-finder based on the Pegasus method.

---

class hydpy.auxs.roottools.Pegasus(model: Model)

Bases: Submodel

Root-finder based on the Pegasus method.

The Pegasus method is a root-finding algorithm which sequentially decreases its search radius (like the simple bisection algorithm) and shows superlinear convergence properties (like the Newton-Raphson algorithm). Our implementation adds a simple “interval widening” strategy that comes into play when the given initial interval does not contain the root. This widening should be left for emergencies as it might be not always efficient. So please try to provide to initial interval estimates that rather overestimate than underestimate the interval width. Additionally, make sure you apply Pegasus on monotone functions only.

CYTHONBASECLASS

alias of PegasusBase

PYTHONCLASS

alias of PegasusPython

find_x(x0: float, x1: float, xmin: float, xmax: float, xtol: float, ytol: float, itermax: int) → float

Find the relevant root within the interval \(x0 \leq x \leq x1\) with an accuracy meeting at least one of the absolute tolerance values xtol and ytol (or the accuracy achieved by performing the maximum number of iteration steps, defined by argument itermax).

The arguments xmin and xmax define the smallest and largest allowed x value, respectively. Method find_x() never evaluates the underlying function for any x value outside these bounds, even if its “interval widening” strategy suggests so. Instead, it returns the relevant xmin or xmax value.

In the following, we explain the details of our Pegasus implementation by using method Return_TempSSurface_V1.

Method Return_TempSSurface_V1 (used by application model lland_knauf) implements the Pegasus iteration to determine the surface temperature of the snow layer with the help of the Pegasus subclass PegasusTempSSurface. For the correct surface temperature, the net energy gain of the surface (defined by method Return_EnergyGainSnowSurface_V1) must be zero. Iteration is necessary, as some terms of the net energy gain (for example, the emitted longwave radiation) depend on the surface temperature in a non-linear manner.

As a starting point, we use the setting provided by the documentation on method Return_TempSSurface_V1 but work in pure Python mode for flexibility (more specifically, to have direct access to method find_x) and focus on a single hydrological response unit for simplicity:

>>> from hydpy import pub
>>> with pub.options.usecython(False):
...     from hydpy.models.lland import *
...     simulationstep("1d")
...     parameterstep("1d")
>>> nhru(1)
>>> lnk(ACKER)
>>> turb0(2.0)
>>> turb1(2.0)
>>> ktschnee(5.0)
>>> inputs.relativehumidity = 60.0
>>> states.waes.values = 1.0
>>> fluxes.tkor = -3.0
>>> fluxes.reducedwindspeed2m = 3.0
>>> fluxes.actualvapourpressure = 2.9
>>> fluxes.netshortwaveradiationsnow = 10.0
>>> aides.temps = -2.0
>>> aides.rlatm = 200.0

Method Return_TempSSurface_V1 finds the following surface temperature value:

>>> model.idx_hru = 0
>>> from hydpy import round_
>>> round_(model.return_tempssurface_v1(0))
-8.307868

We confirm that the net energy gain corresponding to this surface temperature is approximately zero to validate this result:

>>> round_(model.return_energygainsnowsurface_v1(-8.3078682))
0.0

Method apply_method0 always calls the appropriate target method, which is, for PegasusTempSSurface, method Return_EnergyGainSnowSurface_V1:

>>> round_(model.pegasustempssurface.apply_method0(-8.3078682))
0.0

To check for some exceptional cases, we call the find_x method of class PegasusTempSSurface directly. First, we pass the same arguments as in method Return_TempSSurface_V1 and thus get the same result:

>>> round_(model.pegasustempssurface.find_x(
...     -50.0, 5.0, -100.0, 100.0, 0.0, 1e-8, 10))
-8.307868

Through decreasing the maximum number of iterations, the result becomes less accurate:

>>> round_(model.pegasustempssurface.find_x(
...     -50.0, 5.0, -100.0, 100.0, 0.0, 1e-8, 3))
-8.312628
>>> round_(model.return_energygainsnowsurface_v1(-8.312628))
0.09905

Use the ytol argument to control the achieved accuracy more explicitly:

>>> round_(model.pegasustempssurface.find_x(
...     -50.0, 5.0, -100.0, 100.0, 0.0, 0.1, 10))
-8.312628
>>> round_(model.return_energygainsnowsurface_v1(-8.312628))
0.09905

Pass an xtol value larger than zero to stop iteration as soon as the search interval is small enough:

>>> round_(model.pegasustempssurface.find_x(
...     -50.0, 5.0, -100.0, 100.0, 0.1, 1e-8, 10))
-8.30786
>>> round_(model.return_energygainsnowsurface_v1(-8.30786))
-0.000171

x0 does not need to be smaller than x1, and xmin does not need to be smaller than xmax by necessity (the algorithm swaps values when necessary):

>>> round_(model.pegasustempssurface.find_x(
...     50.0, -50.0, 100.0, -100.0, 0.0, 1e-8, 10))
-8.307868

An ill-defined initial interval might decrease efficiency but should usually not affect the result:

>>> round_(model.pegasustempssurface.find_x(
...     0.0, 5.0, -100.0, 100.0, 0.0, 1e-8, 10))
-8.307868
>>> round_(model.pegasustempssurface.find_x(
...     -50.0, -20.0, -100.0, 100.0, 0.0, 1e-8, 10))
-8.307868

Defining proper values for arguments xmin and xmax provides a way to prevent possible program crashes. As an example, we do not care about physics and set the snow layer’s temperature to -300 °C. Theoretically, the surface temperature must also be (unrealistically) low to prevent a considerable energy gain of the surface due to an extreme temperature gradient within the snow layer. However, -100 °C is the lowest value allowed, so method PegasusTempSSurface returns this boundary-value. It is up to the model developer to handle such misleading responses (if possible):

>>> aides.temps = -400
>>> round_(model.pegasustempssurface.find_x(
...     -50.0, 5.0, -100.0, 100.0, 0.0, 1e-8, 10))
-100.0
>>> round_(model.return_energygainsnowsurface_v1(-100.0))
-524.133374

To demonstrate the upper boundary’s functionality, we set the snow layer temperature to 4’000 °C (also not the most plausible value, of course):

>>> aides.temps = 4000.0
>>> round_(model.pegasustempssurface.find_x(
...     -50.0, 5.0, -100.0, 100.0, 0.0, 1e-8, 10))
100.0
>>> round_(model.return_energygainsnowsurface_v1(100.0))
3561.059134

Finally, we evaluate some additional snow layer temperatures to show that everything works consistently:

>>> from hydpy import print_vector
>>> for temps in [-500, -400, -300, -4, -2, 2, 2000, 3000, 4000]:
...    aides.temps = temps
...    tempssurface = model.pegasustempssurface.find_x(
...        -50.0, 5.0, -100.0, 100.0, 0.0, 1e-8, 10)
...    energygain = model.return_energygainsnowsurface_v1(tempssurface)
...    print_vector([temps, tempssurface, energygain])
-500, -100.0, -1024.133374
-400, -100.0, -524.133374
-300, -100.0, -24.133374
-4, -8.790047, 0.0
-2, -8.307868, 0.0
2, -7.353441, 0.0
2000, 85.175508, 0.0
3000, 97.204856, 0.0
4000, 100.0, 3561.059134

apply_method0(value: float) → float

Apply the model method relevant for root-finding.

name: ClassVar[str] = 'pegasus'