# Numerical Methods module¶

Miscellaneous numerical methods in Python.

## collect_run_stats¶

collect_run_stats.py: Run an executable a number of times and report.

This is a small application program rather than a library function. Usage:

```\$ \${E3BIN}/cfpylib/nm/collect_run_stats.py executable no_times
```

Adaptive quadrature using Newton-Cotes 5- and 3-point rules.

Example transcript:

```peterj@laval ~/cfcfd3/lib/cfpylib/nm \$ python adapti.py
Estimates of pi/4:  0.785398154922 0.785398173551
errors: 8.47578540686e-09 -1.01539638919e-08
number of function calls: 315 55
Done.
```
`cfpylib.nm.adapti.``rinteg`(f, a, b, tol)

Apply Newton-Cotes 5- and 3-point quadrature rules to the segment [a,b].

Parameters: f – user-supplied function, f(x) b (a,) – range of integration tol – maximum difference between rules, above which the range is split integral of f(x) from a to b.

## least_squares¶

least_squares.py: Fits a generalized polynomial basis to given data.

Example transcript:

```\$python ~/e3bin/cfpylib/nm/least_squares.py
Should be roughly y = x**2, i.e. coefficients [0,0,1].
coeff= [0.06547619047616493, -0.092142857142826, 1.038095238095233]
```
`cfpylib.nm.least_squares.``least_squares`(x, y, fList)

Fit a generalized polynomial to the given data (x, y).

The fitted model is of the form: y(x) = a*fList(x) + a*fList(x) + ... + a[m]*fList[m](x)

Parameters: x – list or array of data for independent variable y – list or array of data for dependent variable fList – list or array of basis functions to use a vector of coefficients a[]

## nelmin¶

nelmin.py: Nelder-Mead simplex minimization of a nonlinear (multivariate) function.

This code has been adapted from the C-coded nelmin.c which was adapted from the Fortran-coded nelmin.f which was, in turn, adapted from the papers:

J.A. Nelder and R. Mead (1965) A simplex method for function minimization. Computer Journal, Volume 7, pp 308-313.

R. O’Neill (1971) Algorithm AS47. Function minimization using a simplex algorithm. Applied Statistics, Volume 20, pp 338-345.

and some examples are in:

D.M. Olsson and L.S. Nelson (1975) The Nelder-Mead Simplex procedure for function minimization. Technometrics, Volume 17 No. 1, pp 45-51.

For a fairly recent and popular incarnation of this minimizer, see the amoeba function in the famous “Numerical Recipes” text. The programming interface is via the minimize() function; see below.

Example transcript:

```\$ python ~/e3bin/cfpylib/nm/nelmin.py
Begin nelmin self-test...
---------------------------------------------------
test 1: simple quadratic with zero at (1,1,...)
x= [1.0000000651818255, 1.000000021516589, 0.9999999925111813]
fx= 4.76771637998e-15
convergence-flag= 0
number-of-fn-evaluations= 300
number-of-restarts= 3
---------------------------------------------------
test 2: Example 3.3 in Olsson and Nelson f(0.811,-0.585)=-67.1
x= [0.8112948625421268, -0.5846355980866065]
fx= -67.1077410608
convergence-flag= 1
number-of-fn-evaluations= 127
number-of-restarts= 0
---------------------------------------------------
test 3: Example 3.5 in Olsson and Nelson, nonlinear least-squares
f(1.801, -1.842, -0.463, -1.205)=0.0009
x= [1.8010249070374442, -1.8417283432511073, -0.46337704853342615, -1.205053720578973]
fx= 0.000908952916125
convergence-flag= 1
number-of-fn-evaluations= 618
number-of-restarts= 0
---------------------------------------------------
Done.
```
`cfpylib.nm.nelmin.``minimize`(f, x, dx=None, tol=1e-06, maxfe=300, n_check=20, delta=0.001, Kreflect=1.0, Kextend=2.0, Kcontract=0.5)

Locate a minimum of the objective function, f.

Parameters: f – user-specified function f(x) x – list of N coordinates dx – list of N increments to apply to x when forming the initial simplex. Their magnitudes determine the size and shape of the initial simplex. tol – the terminating limit for the standard-deviation of the simplex function values. maxfe – maximum number of function evaluations that we will allow n_check – number of steps between convergence checks delta – magnitude of the perturbations for checking a local minimum and for the scale reduction when restarting Kreflect – Kextend – Kcontract – coefficients for locating the new vertex a tuple consisting of: ``` a list of coordinates for the best x location, corresponding to min(f(x)),  the function value at that point,  a flag to indicate if convergence was achieved  the number of function evaluations and  the number of restarts (with scale reduction) ```

## ode¶

ode.py: Integrate a set of first-order ODEs.

This module provides a small number of ODE integration schemes together with a coordination function ode_integrate(). The schemes will integrate a system of (nonstiff) ODEs over a given range of the independent variable and with a specified step. The user is expected to supply a function that accepts the current value of the independent variable and system state and which provides the derivatives of the system as an array. For further details, see the doc string of ode_integrate() and study the sample functions near the end of this file.

Running this module as a Python script gives me the following transcript:

```Start sample integrations...
(1) Constant derivatives with Euler method:
t1= 10.0
y1= [ 10. -20.  30.]
err1= [ 0.  0.  0.]
(2) Second-order linear ODE with Euler method:
t1= 6.28318530718
y1= [ -2.16009598e-04   1.03190177e+00]
err1= [ 0.02031786  0.02030512]
(3) Second-order linear ODE with RKF45 method:
t1= 6.28318530718
y1= [  8.42476994e-14   1.00000000e+00]
err1= [  5.12421830e-11   5.12821769e-11]
Done.
```
`cfpylib.nm.ode.``ode_integrate`(t0, tlast, dt, F, n, y0, method='euler')

Steps the set of ODEs until t reaches tlast.

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 one of a set of more specialised stepping functions that appear below.

Parameters: t0 – is the starting value of the independent variable tlast – the desired finishing value for x dt – 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 method – a string specifying which stepping method to use. “euler”, “rkf45” final values of t, y, and error estimates for y values are returned as a tuple.
`cfpylib.nm.ode.``euler_step`(t0, h, F, n, y0)

Single steps the set of ODEs by the Euler method.

Parameters: 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 final values of t, y, and error estimates for y values are returned as a tuple.
`cfpylib.nm.ode.``rkf45_step`(t0, h, F, n, y0)

Steps the set of ODEs by the Runge-Kutta-Fehlberg method.

Parameters: 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 final values of t, y, and error estimates for y values are returned as a tuple.

## roberts¶

roberts.py: Node distribution and coordinate stretching functions.

These functions should behave the same as the C code functions.

`cfpylib.nm.roberts.``roberts`(eta, alpha, beta)

Computes the stretched coordinate in the range [0.0..1.0] using the boundary-layer-like transformation devised by Roberts.

Parameters: eta – unstretched coordinate, 0 <= eta <= 1.0 beta – stretching factor (more stretching as beta –> 1.0) alpha – location of stretching: alpha = 0.5: clustering of nodes at both extremes of eta alpha = 0.0: nodes will be clustered near eta=1.0

Works for both scalars and arrays.

`cfpylib.nm.roberts.``distribute_points_1`(t1, t2, n, end1, end2, beta)

Generate a set of points nonuniformly distributed from t1 to t2.

Parameters: t1 – parameter value 1 t2 – parameter value 2 n – number of intervals (the number of points is n+1) end1 – (int) cluster flag for end 1: ==0 points are not clustered to end 1 ==1 points ARE clustered to end 1 end2 – (int) cluster flag for end 2: ==0 points are not clustered to end 2 ==1 points ARE clustered to end 2 beta – grid stretching parameter: 1 < beta < +inf : points are clustered The closer to 1, the more the clustering. beta < 1 for no clustering. t[0:n] an array of distributed values.

## secant_method¶

secant_method.py: Function solver, using the secant method.

Example transcript:

```\$ python ~/e3bin/cfpylib/nm/secant_method.py
Begin secant_method test...
one solution at x= 6.28318530718
expected x= 6.28318530718
Done.
```
`cfpylib.nm.secant_method.``solve`(f, x1, x2, tol=1e-09)

Computes x that satisfies f(x) = 0.

Parameters: f – user-defined function x1 – first guess x2 – second guess, presumably close to x1 tol – stopping tolerance on f(x)=0 x for best solution

## sode¶

sode.py: Integrate a set of stiff ODEs.

Example transcript:

```\$ python ~/e3bin/cfpylib/nm/sode.py
Start sample integrations...
(1) Constant derivatives with semi-implicit RK1 method:
t1= 10.0
y1= [ 10. -20.  30.]
err1= [ 0.  0.  0.]
(2) Second-order linear ODE with semi-implicit RK1 method:
t1= 6.28318530718
y1= [ -5.23352415e-05   9.99999999e-01]
err1= [ 0.01999997  0.01998885]
(3) Second-order linear ODE with semi-implicit RK2 method:
t1= 6.28318530718
y1= [ -2.29779805e-10   1.00000000e-00]
err1= [ 0.01999992  0.01998906]
Done.
```
`cfpylib.nm.sode.``ode_integrate`(t0, tlast, dt, F, n, y0, dFdt=None, dFdy=None, method='rk1')

Steps the set of ODEs until x reaches xlast.

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 one of a set of more specialised stepping functions that appear below.

Parameters: t0 – the starting value of the independent variable tlast – the desired finishing value for x dt – 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 dFdt – user-defined function that returns the vector of partial derivatives dFdy – user-defined function that returns the matrix of partial derivatives method – a string specifying which stepping method to use. “rk1”, “rk2” final values of t, y, and error estimates for y values as a tuple.
`cfpylib.nm.sode.``rk1_step`(t0, h, F, n, y0, dFdt, dFdy)

Steps the set of ODEs by the semi-implicit one-stage Runge-Kutta method.

Parameters: t0 – the starting value of the independent variable h – the step size in t 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 dFdt – user-defined function that returns the vector of partial derivatives dFdy – user-defined function that returns the matrix of partial derivatives final values of t, y, and error estimates for y values as a tuple.
`cfpylib.nm.sode.``rk2_step`(t0, h, F, n, y0, dFdt, dFdy)

Steps the set of ODEs by the semi-implicit two-stage Runge-Kutta method.

Parameters: t0 – the starting value of the independent variable h – the step size in t 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 dFdt – user-defined function that returns the vector of partial derivatives dFdy – user-defined function that returns the matrix of partial derivatives final values of t, y, and error estimates for y values as a tuple.

## stats¶

stats.py: Simple statistics for arrays of values.

To replace those in scipy, just in case scipy is not installed.

PJ, 27-Mar-2007

`cfpylib.nm.stats.``mean`(a)

Return the mean value of an array.

`cfpylib.nm.stats.``std`(a)

Return the standard deviation of an array of values.

## zero_finding¶

Used in Brendan’s free-piston engine modelling code.

zero_finding.py

Example transcript:

```\$ python ~/e3bin/cfpylib/nm/zero_finding.py
bisection method
find x for f0(x) < 1.000000e-06 for x(2.20:4.08)
f0(3.141593) = 4.494084e-07

golden section method
find x for min(f1(x)) < 1.000000e-06 for x(0.94:2.20)
f1(1.570796) = -1.000000
```
`cfpylib.nm.zero_finding.``bisection_root`(f, x0, x1, tol, args=(), max_iter=100)

Applies the Bisection-Root method to find f(x) = 0.

`cfpylib.nm.zero_finding.``muller_root`(f, x0, x1, tol, args=(), max_iter=100)

Applies the Bisection-Root method to find f(x) = 0.

`cfpylib.nm.zero_finding.``newton_root`(f, x0, x1, tol, args=())

Applies Newton’s method to find f(x) = 0.

`cfpylib.nm.zero_finding.``golden_section`(f, a, c, tol, args=())

Applies the Golden-Section search to find min(f(x)).

## zero_solvers¶

zero_solvers.py

Example transcript:

```python ~/e3bin/cfpylib/nm/zero_solvers.py
Begin zero_solvers self-test...

Test function 1.
----------------
Example from Gerald and Wheatley, p. 45
Solve f(x) = x^3 + x^2 - 3x -3 = 0 with initial
guesses of x0 = 1 and x1 = 2.
Begin function call secant()...

Iteration    x0             x1              x2       f(x2)
-----------------------------------------------------------------------
1     1.000000       2.000000        1.571429        -1.364431e+00
2     2.000000       1.571429        1.705411        -2.477451e-01
3     1.571429       1.705411        1.735136        2.925540e-02
4     1.705411       1.735136        1.731996        -5.151769e-04
5     1.735136       1.731996        1.732051        -1.039000e-06
6     1.731996       1.732051        1.732051        3.702993e-11
7     1.732051       1.732051        1.732051        1.776357e-15
-----------------------------------------------------------------------
Final result x =  1.73205080757
Gerald and Wheatley report x = 1.732051
Using bisection... x = 1.73205080757

Test function 2.
----------------
Example from Gerald and Wheatley, p.45
Solve f(x) = 3*x + sin(x) - e^x = 0 with initial
guesses of x0 = 0 and x1 = 1.
Begin function call secant()...

Iteration    x0             x1              x2       f(x2)
-----------------------------------------------------------------------
1     1.000000       0.000000        0.470990        2.651588e-01
2     0.000000       0.470990        0.372277        2.953367e-02
3     0.470990       0.372277        0.359904        -1.294813e-03
4     0.372277       0.359904        0.360424        5.530053e-06
5     0.359904       0.360424        0.360422        1.021329e-09
6     0.360424       0.360422        0.360422        -8.881784e-16
-----------------------------------------------------------------------
Final result x =  0.36042170296
Gerald and Wheatley report x = 0.3604217
Using bisection... x = 0.360421702961
```
`cfpylib.nm.zero_solvers.``secant`(f, x0, x1, tol=1e-11, limits=[], max_iterations=1000)

The iterative secant method for zero-finding in one-dimension.

Parameters: f – user-defined function f(x) x0 – first guess x1 – second guess, presumably close to x0 stopping tolerance for f(x)=0 to stop the iterations running forever, just in case... x such that f(x)=0 or the string ‘FAIL’
`cfpylib.nm.zero_solvers.``bisection`(f, bx, ux, tol=1e-06)

The iterative bisection method for zero-finding in one-dimension.

Parameters: f – user-defined function f(x) bx – bottom-limit of bracket ux – upper-limit of bracket stopping tolerance on bracket size x such that f(x)=0