Chapter 2. Ax=b

Linear equation solution by Gauss-Jordan elimination... a(1:n,1:n) is an input matrix stored in an array of physical dimensions np by np. b(1:n,1:m) is an input matrix containing the m right-hand side vectors, stored in an array of physical dimensions np by mp. On output, a(1:n,1:n) is replaced by its matrix inverse, and b(1:n,1:m) is replaced by the corresponding set of solution vectors. (My notes: Can handle up to dimension=50. NR bug: mp is never used; thus, it can be anything -- simply set "mp=1".)

Chapter 3. Interpolation and Extrapolation

Given arrays xa and ya, each of length n, and given a value x, this routine returns a value y, and an error estimate dy.

Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., y_i=f(x_i), with x₁<x₂<...<x_n, and given values yp1 and ypn for the first derivative of the interpolating function at points 1 and n, respectively, this routine returns an array y2(1:n) of length n which contains the second derivatives of the interpolating function at the tabulated points x_i. If yp1 and/or ypn are equal to 1X10³⁰ or larger, the routine is signaled to set the corresponding boundary conditions for a natural spline, with zero second derivative on that boundary.

Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a function (with x₁<x₂<...<x_n), and given the array y2a(1:n), which is the output from spline above, and given a value of x, this routine returns a cubic-spline interpolated value y.

Chapter 4. Integrals

This routine computes the nth stage of refinement of an extended trapezoidal rule. func is input as the name of the function to be integrated between limits a and b, also input. When called with n=1, the routine returns as s the crudest estimate of ò_a^b f(x)dx. Subsequent calls with n=2,3,... (in that sequential order) will improve the accuracy of s by adding 2^(n-2) additional interior points. s should not be modified between sequential calls. (My note, be sure to declare "external func" in the calling program. Do not simply call with, say, n=5; you MUST loop from n=1 to n=5.)

Returns as s the integral of the function func from a to b. ... Integration is performed by the trapezoidal rule.

Romberg integration on an open interval. Returns as ss the integral of the function func from a to b, using any specified integrating subroutine choose and Romberg's method. Normally choose will be an open formula, not evaluating the function at the endpoints. It is assumed that choose triples the number of steps on each call, and that its error series contains only even powers of the number of steps. The routines midpnt, midinf, midsql, midsqu, are possible choices for choose. The parameters have the same meaning as in qromb.

This routine is an exact replacement for midpnt, except that bb is assumed to be infinite (value passed not actually used). It is assumed that the function funk decreases exponentially rapidly at infinity.

Chapter 5. Evaluation of Functions

Given the coefficients of a polynomial of degree nc-1 as an array c(1:nc) with c(1) being the constant term, and given a value x, and given a value nd>1, this routine returns the polynomial evaluated at x ad pd(1) and nd-1 derivatives as pd(2:nd)

Chapter 9. f(x)=0

Using bisection, find the root of a function func known to lie between x1 and x2. The root, returned s rtbis, will be refined until its accuracy is xacc.

Find all the roots of a polynomial with real coefficients, S_(i=1)^(m+1) a(i)*x^(i-1), given the degree m and the coefficients a(1:m+1). The method is to construct an upper Hessenberg matrix whose eigenvalues are the desired roots, and then use the routines balanc and hqr. The real and imaginary parts of th roots are returned in rtr(1:m) and rti(1:m), respectively.

Given the degree m and the complex coefficients a(1:m+1) of the polynomial S_(i=1)^(m+1) a(i)*x^(i-1), this routine ... finds all m complex roots. The logical variable polish should be input as .true. if polishing is desired, .false. if the roots will be subsequently polished by other means.

Given an initial guess x for a root in n dimensions, take ntrial Newton-Raphson steps to improve the root. Stop if the root converges in either summed absolute variable increments tolx or summed absolute function values tolf. (My notes: The vector function f(x)=0 must be specified in subroutine usrfun(x, dfdx, f); the user does not get to choose a name. NR should have make the function name one of the arguments to MNEWT. Furthermore, the leading dimension of dfdx MUST be 15, which should have been a variable if the subroutine was written professionally.)

Chapter 10. Minimization

Given a function func, and given distinct initial points ax and bx, this routine searches in the downhill direction (defined by the function as evaluated at the initial points) and returns new points ax, bx, cx that bracket a minimum of the function. Also returned are the function values at the three points, fa, fb, and fc.

Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is between ax, and cx, and f(bx) is less than both f(ax) and f(cx)), this routine performs a golden section search for the minimum, isolating it to a fractional precision of about tol. The abscissa of the minimum is returned as xmin, and the minimum function value is returned as golden, the returned function value.

Minimization of a function func of n variables. (func is not an argument, it is a fixed function name.) Input consists of an initial starting point p(1:n); an initial matrix xi(1:n,1:n) with physical dimensions np by np, and whose columns contain the initial set of directions (usually, the n unit vectors); and ftol, the fractional tolerance in the function value such that failure to decrease by more than this amount on one iteration signals doneness. On output, p is set to the best point found, xi is the then-current direction set, fret is the returned function value at p, and iter is the number of iterations taken. The routine linmin is used.

Given a starting point p that is a vector of length n, Fletcher-Reeves-Polak-Ribiere minimization is performed on a function func, using its gradient as calculated by a routine dfunc. The convergence tolerance on the function value is input as ftol. Returned quantities are p (the location of the minimum), iter (the number of iterations that were performed), and fret (the minimum value of the function). The routine linmin is called to perform line minimizations. (Note: func(p) is not an argument, it is a fixed function name. dfunc(p,xi) is not an argument.)

Chapter 15. Regression

Given a set of data points x(1:ndat), y(1:ndat) with individual standard deviations sig(1:ndat), use chi² minimization to fit for some or all of the coefficients a(1:ma) of a function that depends linearly on a, y=sum a_i*afunc_i(x). The input array ia(1:ma) indicates by nonzero entries those components of a that should be fitted for, and by zero entries those components that should be held fixed at their input values. The program returns values for a(1:m), chi²=chisq, and the covariance matrix covar(1:ma,1:ma). (Parameters held fixed will return zero covariances.) npc is the physical dimension of covar(npc,npc) in the calling routine. The user supplies a subroutine funcs(x,afunc,ma) that returns the ma basis functions evaluated at x in the array afunc.

Levenberg_Marquardt method, attempting to reduce the value X² of a fit between a set of data points x(1:ndata), y(1:ndata) with individual standard deviations sig(1:ndata), and a nonlinear function dependent on ma coefficients a(1:ma). The input array ia(1:ma) indicates by nonzero entries those components of a that should be fitted for, and by zero entries those components that should be held fixed at their input values. The program returns current best-fit values for the parameters a(1:ma), and X²=chisq. The arrays covar(1:nca,1:nca), alpha(1:nca,1:nca) with physical dimension nca(>= the number of fitted parameters) are used as working space during most iterations. Supply a subroutine funcs(x,a,yfit,dyda,ma) that evaluates the fitting function yfit, and its derivatives dyda with respect to the fitting parameters a at x. On the first call provide an initial guess for the parameters a, and set alamda<0 for initialization (which then sets alamda=0.001). If a step succeeds, chisq becomes smaller and alamda decreases by a factor of 10. If a step fails, alamda grows by a factor of 10. You must call this routine repeatedly until convergence is achieved. Then make one final call with alamda=0, so that covar(1:ma,1:ma) returns the covariance matrix, and alpha the curvature matrix. (Parameters held fixed will return zero covariances.)

Chapter 16. ODE

Given values for the variables y(1:n) and their derivatives dydx(1:n) known at x, use the fourth-order Runge_Kutta method to advance the solution over an interval h and return the incremented variables as yout(1:n), which need not be a distinct array from y. The user supplies the subroutine derivs(x,y,dydx), which returns derivatives dydx at x. (Note: rk4 makes only three calls to derivs; thus, you must make the first call to derivs yourself.)

Return to Prof. Nam Sun Wang's Home Page
Return to Computer Methods in Chemical Engineering (ENCH250)

Computer Methods in Chemical Engineering -- Selected Routines from Numerical Recipes
Forward comments to:

Nam Sun Wang Department of Chemical & Biomolecular Engineering University of Maryland College Park, MD 20742-2111 301-405-1910 (voice) 301-314-9126 (FAX) e-mail: nsw@umd.edu ©1996-2010 by Nam Sun Wang