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

adapti

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

Example transcript:

peterj@laval ~/cfcfd3/lib/cfpylib/nm $ python adapti.py 
Begin adapti...
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
Returns:

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[0]*fList[0](x) + a[1]*fList[1](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
Returns:

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
Returns:

a tuple consisting of:

[0] a list of coordinates for the best x location, corresponding to min(f(x)),
[1] the function value at that point,
[2] a flag to indicate if convergence was achieved
[3] the number of function evaluations and
[4] 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”
Returns:

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
Returns:

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
Returns:

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.
Returns:

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
Returns:

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”
Returns:

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
Returns:

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
Returns:

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
Tol:

stopping tolerance for f(x)=0

Max_iterations:

to stop the iterations running forever, just in case...

Returns:

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
Tol:

stopping tolerance on bracket size

Returns:

x such that f(x)=0

Table Of Contents

Previous topic

cfpylib (Library of support modules)

Next topic

Gas Dynamics module

This Page