# 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.