Numerical methods library

This library provides access to various numerical methods that we have found useful for the construction of our gas-dynamics tools. The import library sits in the gdtk.numeric package.

1. Installing the library

This numerical-methods library for Python3 is part of a larger gas-dynamics toolkit and general getting started notes can be found at https://gdtk.uqcloud.net/docs/getting-started/prerequisites . There, you will see how to get a copy of the source code, and a list of what other software you will need to build and install the tool kit, and a collection of environment variables that need to be set.

To install the library, move to the gas source directory and use the make utility.

cd dgd/src/gas
make install

Even though this part of the package is a pure Python library, the rest of the loadable library, including gas models, will be built and installed with this command. So that the Python interpreter can find the installed library, set your environment variables with something like:

export DGD=$HOME/dgdinst
export PYTHONPATH=${PYTHONPATH}:${DGD}/lib

2. Solve scalar f(x)=0

from gdtk.numeric.zero_solvers import secant, bisect

2.1. Secant method

secant(f, x0, x1, tol=1.0e-11, limits=[], max_iterations=1000, tf=False): Returns x such that f(x)=0.
f: User-supplied function that accepts a scalar value x.
x0: first guess
x1: second guess, presumably close to x0
tol: stopping tolerance for f(x)=0
max_iterations: to stop the iterations running forever, just in case…
tf: boolean flag to turn on printing of intermediate states

2.2. Bisection method

bisection(f, bx, ux, tol=1.0e-6): Returns x such that f(x)=0.
f: User-supplied function that accepts a scalar value x.
bx: bottom-limit of bracket
ux: upper-limit of bracket
tol: stopping tolerance on bracket size

3. Minimize f(x)

from gdtk.numeric.nelmin import minimize

minimize(f, x, dx=None, options={})

Locate a minimum of the objective function, f.

f

user-specified scalar function f(x) of a list of parameters, x

x

list of N coordinates in parameter space.

dx

optional list of N increments to apply to x when forming the initial simplex. These increments determine the size and shape of the initial simplex.

options, a dictionary with entries

tol: (default 1.0e-6) the terminating limit for the standard-deviation of the simplex function values.
P: (default 1) number of points to replace in parallel, each step.
n_workers: (default 1) number of concurrent threads or processes in pool
maxfe: (default 300) maximum number of function evaluations that we will allow
n_check: (default 20) number of steps between convergence checks
delta: (default 0.001) magnitude of the perturbations for checking a local minimum and for the scale reduction when restarting
Kreflect: (default 1.0)
Kextend: (default 2.0)
Kcontract: (default 0.5) coefficients for locating the new vertex

Returns a namedtuple consisting of

x, a list of coordinates for the best x location, corresponding to min(f(x)),
fun, the function value at that point,
success, a flag to indicate if convergence was achieved
nfe, the number of function evaluations and
nrestarts, the number of restarts (with scale reduction)

3.1. Example

from gdtk.numeric.nelmin import minimize

def test_fun(x):
    "Example 3.3 from Olsson and Nelson."
    x1, x2 = x   # rename to match the paper
    if (x1 * x1 + x2 * x2) > 1.0:
        return 1.0e38
    else:
        yp = 53.69 + 7.26 * x1 - 10.33 * x2 + 7.22 * x1 * x1 \
             + 6.43 * x2 * x2 + 11.36 * x1 * x2
        ys = 82.17 - 1.01 * x1 - 8.61 * x2 + 1.40 * x1 * x1 \
             - 8.76 * x2 * x2 - 7.20 * x1 * x2
        return -yp + abs(ys - 87.8)

print("Example 3.3 in Olsson and Nelson f(0.811,-0.585)=-67.1")
result = minimize(test_fun, [0.0, 0.0], [0.5, 0.5], options={'tol':1.0e-4})
print("  x=", result.x)
print("  fx=", result.fun)
print("  convergence-flag=", result.success)
print("  number-of-fn-evaluations=", result.nfe)
print("  number-of-restarts=", result.nrestarts)

4. Integrate ODEs

from gdtk.numeric.ode import ode_integrate, rk45_step

4.1. Integrate to stopping point

ode_integrate(t0, tlast, nstep, f, n, y0): Steps the set of ODEs until independent variable, t, reaches tlast. Returns lists of t, y, and error estimates for y values in a tuple.

This function coordinates the work of integrating a system of first-order differential equations of the form: y'=f(t, y); y(t=t0)=y0 The actual work is done by rkf45_step, a more specialised stepping function, that is described below.

t0: is the starting value of the independent variable
tlast: the desired finishing value for x
nstep: number of steps to take to arrive at tlast
f: a callable function that returns the derivative of y wrt t. The signature of this function is f(t, y, n) where t is a float value, y is an array of float values and n is an integer specifying the number of equations.
n: the number of dependent variables (in y)
y0: an array of starting values for the dependent variables. It is assumed that the y-elements are indexed 0..n-1.

4.2. Single step

rkf45_step(t0, h, f, n, y0): Single-step the set of ODEs by the Runge-Kutta-Fehlberg method. Returns final values of t, y, and error estimates for y values in a tuple.
t0: is the starting value of the independent variable
h: the requested step size
f: a callable function that returns the derivative of y wrt t. The signature of this function is f(t, y, n) where t is a float value, y is a list (or array) or float values and n is an integer specifying the number of equations.
n: the number of dependent variables (in y)
y0: an array of starting values for the dependent variables. It is assumed that the y-elements are indexed 0..n-1.

Numerical methods library

1. Installing the library#

2. Solve scalar f(x)=0#

2.1. Secant method#

2.2. Bisection method#

3. Minimize f(x)#

3.1. Example#

4. Integrate ODEs#

4.1. Integrate to stopping point#

4.2. Single step#