Developer Guide
When you write developer guide or docs in this project, please be concise.
Mathematics
All mathematics below is written without assuming a uniform grid for the functions, and the current library code also works for non-uniform grids.
Trapezoid Rule for Indefinite Integral
This is used in GridFunction.integrate(). Given grid function \(\{x_j\}_{j=0}^n, \{y_j\}_{j=0}^n\), we derive the grid function of its integral, denoted by \(\{x_j\}_{j=0}^n, \{\tilde y_j\}_{j=0}^n\). To begin with, clearly we will have \(\tilde y_0 = 0\).
Note that the area of the first small trapezoid is \((y_1 + y_0)(x_1 - x_0)/2\), and the area of the second small trapezoid is \((y_2 + y_1)(x_2 - x_1)/2\), and so on. Thus we have a list of small trapezoid areas
\begin{align*} 0, \quad\frac{1}{2}(y_1 + y_0)(x_1 - x_0), \quad\frac{1}{2}(y_2 + y_1)(x_2 - x_1), \quad\ldots, \quad\frac{1}{2}(y_n + y_{n-1})(x_n - x_{n-1}), \end{align*}
the cumulative sum of which is the desired \(\{\tilde y_j\}_{j=0}^n\).
Differentation Operators
Below is from p.57 of Interest Rate Modeling by Leif Andersen and Vladimir Piterbarg.
Define
\begin{align*} \Delta^+_j = x_{j+1} - x_j, \qquad\Delta^-_j = x_j - x_{j-1}. \end{align*}
The forward and backward difference operators are, respectively,
\begin{align*} \delta^+_x f(x_j) = \frac{f(x_{j+1}) - f(x_j)}{\Delta^+_j}, \qquad \delta^-_x f(x_j) = \frac{f(x_j) - f(x_{j-1})}{\Delta^-_j}. \end{align*}
As an approximation to \(f'(x_j)\), they have errors \(O(\Delta^+_j)\) and \(O(\Delta^-_j)\), respectively. A second-order approximation of \(f'(x_j)\) is given by
\begin{align*} \delta_x f(x_j) = \frac{\Delta_j^-}{\Delta_j^+ + \Delta_j^-} \delta^+_x f(x_j) + \frac{\Delta_j^+}{\Delta_j^+ + \Delta_j^-} \delta^-_x f(x_j). \end{align*}
It is second-order accurate in the sense that reducing both \(\Delta_j^+\) and \(\Delta_j^-\) by a factor of \(k\) will reduce the error by a factor of \(k^2\). An approximation of \(f''(x_j)\) is
\begin{align*} \delta_{xx} f(x_j) = \frac{\delta^+_x f(x_j) - \delta^-_x f(x_j)}{\frac{1}{2}(\Delta_j^+ + \Delta_j^-)}, \end{align*}
which is only first-order accurate, unless \(\Delta_j^+ = \Delta_j^-\). Despite this, the global discretization error will typically remain second-order even for a nonuniform grid.
Derivatives at End Points
Currently a GridFunction’s derivative is computed by the central difference at the interior points and one-sided difference at the end points. As such, it is second-order accurate in the interior but first-order accurate at the end points. There are at least two ways to potentially achieve second-order accuracy at end points:
Take 3 points at the end and find the derivative of the interpolating parabola at the end point, and
Extrapolate from derivative at the interior points by the central difference.
Here we illustrate the first method. Given the function values at the first 3 grid points, namely the values \(\{x_j\}_{i=0}^2\) and \(\{y_j\}_{i=0}^2\), interpolating quadratic function is \begin{align*} y_0 \frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} + y_1 \frac{(x - x_0)(x - x_2)}{(x_1 - x_0)(x_1 - x_2)} + y_2 \frac{(x - x_0)(x - x_1)}{(x_2 - x_0)(x_2 - x_1)}, \end{align*} whose first order derivative is \begin{align*} y_0 \frac{(x - x_1) + (x - x_2)}{(x_0 - x_1)(x_0 - x_2)} + y_1 \frac{(x - x_0) + (x - x_2)}{(x_1 - x_0)(x_1 - x_2)} + y_2 \frac{(x - x_0) + (x - x_1)}{(x_2 - x_0)(x_2 - x_1)}. \end{align*} Plug in \(x_0\) to get \begin{align*} y_0 \frac{(x_0 - x_1) + (x_0 - x_2)}{(x_0 - x_1)(x_0 - x_2)} + y_1 \frac{x_0 - x_2}{(x_1 - x_0)(x_1 - x_2)} + y_2 \frac{x_0 - x_1}{(x_2 - x_0)(x_2 - x_1)}. \end{align*}
Notes on Design
The Grid Class
In early versions GridFunctions simply store the \(x\) and \(y\) values in two Numpy arrays x and y, but now the \(x\) values are in a Grid object for
Encapsulation: Grid logic is separated from function values for clarity and maintainability
Performance: Grid differences are cached, avoiding repeated calculations in integral and differentition
Extensibility: Future grid features can be added without changing function code
Cleaner API:
GridFunction.xgives the grid array,GridFunction.gridgives the grid object—no moreself.x.xConsistency: All grid-based operations use a single grid instance, reducing errors
Extend the Grid class for grid logic; use GridFunction.x for the array and GridFunction.grid for advanced features.
Miscellaneous
Dynamic Math Function Dispatch
All univariate Numpy math functions (sin, cos, tan, exp, log, sqrt, etc.) are supported automatically, as GFuncPy uses Python’s module-level __getattr__ to catch calls like gfuncpy.sin. If the function exists in NumPy, it creates a wrapper that applies it to the grid values. No manual function definitions needed.
If you define a custom function in GFuncPy with the same name as a NumPy function (e.g., maximum, minimum), you must explicitly import it in __init__.py. This prevents the dynamic dispatch mechanism from interfering and ensures your function is used instead of NumPy’s version.
Usage
Call any NumPy math function from gfuncpy on a GridFunction object, e.g. gfuncpy.sin(x).
max and min Naming Issue
The functions max and min in GFuncPy overwrite the names of Python’s built-in max and min. Use these with caution, especially in code that also relies on the built-in versions. To avoid confusion or bugs, call them via the module (e.g., gfuncpy.max) if you want to use them with the built-in max and min.
Release Management Tip
Before creating a new git tag for this project, always update the release number in both docs/source/conf.py and pyproject.toml to match the new version.
PyData Theme Issue
The pydata-sphinx-theme has an issue that equations will not be centered, which is claimed to have been fixed here, but I tried the latest version (0.16.1) and no luck, so switching back to sphinx-rtd-theme for now.
To-Do List
Docstring and sphinx autodoc
Integration logic not documented (refactor needed?)
Derivative is not even mentioned in the usage page