processmodeling.org

 Mathematics

cos.cpp - Compute the cosine of x to within tolerance eps
asin.cpp - asin(arg) and acos(arg) return the arcsin, arccos, respectively of their arguments. Arctan is called after appropriate range reduction.
atan.cpp - atan returns the value of the arctangent of its argument in the range [-pi/2,pi/2]. atan2 returns the arctangent of arg1/arg2 in the range [-pi,pi].
tan.cpp - Floating point tangent. A series is used after range reduction.
erf.cpp - Floating point error function erf(x) returns the error function of its argument erfc(x) returns 1.0-erf(x) erf(x) is defined by \${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt\$ the entry for erfc is provided because of the extreme loss of relative accuracy if erf(x) is called for large x and the result subtracted from 1. (e.g. for x= 10, 12 places are lost).
exp.cpp , exp2.cpp - exp returns the exponential function of its floating-point argument.
floor.cpp - floor and ceil--greatest integer <= arg (resp least >=)
frexp.cpp - the call x = frexp(arg, exp); must return a double fp quantity x which is 1.0 and the corresponding binary exponent "exp". such that arg = x*2^exp
gamma.cpp - gamma(x) computes the log of the absolute value of the gamma function. The sign of the gamma function is returned in the external quantity signgam.
hypot.cpp - sqrt(a^2 + b^2)
tanh.cpp tanh(arg) computes the hyperbolic tangent of its floating point argument. sinh and cosh are called except for large arguments, which would cause overflow improperly.
sqrt.cpp - Returns the square root of its floating point argument. Newton's method.
sinh.cpp - sinh(arg) returns the hyperbolic sine of its floating-point argument. The exponential function is called for arguments greater in magnitude than 0.5. A series is used for arguments smaller in magnitude than 0.5.
sin.cpp - floating point sin/cos.
pow.cpp - computes a^b. Uses log and exp.
log.cpp - Returns the natural logarithm of its floating point argument.
j0.cpp - Floating point Bessel's function of the first and second kinds of order zero. j0(x) returns the value of J0(x) for all real values of x. y0(x) returns the value of Y0(x) for positive real values of x.
j1.cpp - Floating point Bessel's function of the first and second kinds of order one j1(x) returns the value of J1(x) for all real values of x.
jn.cpp Floating point Bessel's function of the first and second kinds and of integer order. Returns the value of Jn(x) for all integer values of n and all real values of x.
bessel fun.cpp , butil.h - Returns the Bessel functions All of them. Pick up our choice.
runge.cpp , butil.h - Runge-Kutta solution for harmonic oscillator. In this case, d^y2/dx^2 = - y, can be written as dy_1/dx = y_2 and dy_2 /dx = - y_1, with y_1 = y, and y_2 = dy/dx.

runge kutta.cpp - This program uses a forth order Runge-Kutta method to solve a second order differential equation. The equation is for a damped harmonic oscilltor. When the program starts the equation is set up and solved. The solution curve is stored so that it can be plotted by the paint routine. The damped harmonic oscillator equation can be written as:
m * x'' + b * x' + k * x = 0

where m is the mass of the oscillator, b is the damping coefficient and k is the spring constant. We are give the initial position and velocity of the mass and we wish to find x(t) that satisfies these conditions. We expect x(t) to be oscillatory so we want to use enough t values to cover several cycles of the oscillation. The angular frequency w of the oscillator is given by:

w = sqrt((k / m) - (b / 2 m)^2)

(We are assuming that k, m and b are greater than zero and that 4 * k * m > b^2.) Thus we take tmin = 0 and tmax = n * w, for some positive integer n > 1. To use the Runge-Kutta method we reduce this second order equation to a pair of first order equations. We do this by introducing a new variable v = x'. (In physical terms this is making the velocity an explicite term.)

x' = v
v' = -((b / m) * v + (k / m) * x)

Given a system of equations:

x' = f(t, x, v)
v' = F(t, x, v)

and the initial conditions x0 and v0 at t0, the Runge-Kutta method of finding x1 and v1 at t0 + dt is done by finding the following four pairs of numbers:

k1 = dt * f(t, x, v)
l1 = dt * F(t, x, v)
k2 = dt * f(t + dt / 2, x + k1 / 2, v + l1 / 2)
l2 = dt * F(t + dt / 2, x + k1 / 2, v + l1 / 2)
k3 = dt * f(t + dt / 2, x + k2 / 2, v + l2 / 2)
l3 = dt * F(t + dt / 2, x + k2 / 2, v + l2 / 2)
k4 = dt * f(t + dt, x + k3, x + l3)
l4 = dt * F(t + dt, x + k3, x + l3)

These items are then used in:
x1 = x0 + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
v1 = v0 + l1 / 6 + l2 / 3 + l3 / 3 + l4 / 6

These calculations are then repeated with x1 and v1 replacing x0 and y0 to find x2 and v2 at t2 = t1 + dt. This is repeated until all the desired values are found. For what we are doing here:

f(t, x, v) = v and F(t, x, v) = -(b * v + k * x) / m

The k's and l's become:
k1 = dt * v0
l1 = -dt * (b * v0 + k * x0) / m
k2 = dt * (v0 + l1 / 2)
l2 = -dt * (b * (v0 + l1 / 2) + k * (x0 + k1 / 2)) / m
k3 = dt * (v0 + l2 / 2)
l3 = -dt * (b * (v0 + l2 / 2) + k * (x0 + k2 / 2)) / m
k4 = dt * (v0 + l3)
l4 = -dt * (b * (v0 + l3) + k * (x0 + k3)) / m

These are then used to compute:
x1 = x0 + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
v1 = v0 + l1 / 6 + l2 / 3 + l3 / 3 + l4 / 6

rk45.cpp - Runge-Kutta 45

Example:
Enter value for x:1; Enter value for y:2; Enter value for Vx:3; Enter value for Vy:4

 N X Y Vx Vy 0 x = 1 y = 2 Vx = 3 Vy = 4 1 x = -2147483648 y = 402 Vx = -1073741821 Vy = 4 2 x = -2147483648 y = 802 Vx = -2147483648 Vy = 4 3 x = -2147483648 y = -2147483648 Vx = -2147483648 Vy = 177585892 4 x = -2147483648 y = -1568763632 Vx = -2147483648 Vy = 177585892 5 x = -2147483648 y = -2147483648 Vx = -2147483648 Vy = -811648420 6 x = -2147483648 y = -1707947024 Vx = -2147483648 Vy = -811648420 7 x =-2147483648 y =-2147483648 Vx =-2147483648 Vy =-1287506838 8 x = -2147483648 y = -2049148568 Vx = -2147483648 Vy = -1287506838 9 x = -2147483648 y = -2147483648 Vx = -2147483648 Vy = 255514688 10 x = -2147483648 y = 1929148672 Vx = -2147483648 Vy = 255514688

airy.cpp , butil.h - Computes Airy fnctions Ai(x), Bi(x), and their derivatives, for any real x.
dfsafe.cpp , butil.h - Returns the derivative of a function func at a point x by Ridders' method of polynomial extrapolation. The value h is input as an estimated initial stepsize; it need not to be small, but rather should be an increment in x over which func changes "substantially". An estimate of the error in the derivative is returned as err.
error fun.cpp - Returns the error function erf(x)=2*(int_0^x e^{-t^2}dt)/sqrt(pi)
exp integral.cpp - Returns the exponential integral functions E_n(x) and E_i(x). E_n(x) = int_1^infinity e^{-x*t}/t^n dt, for x > 0.
expdev.cpp - Returns an exponentially distributed, positive, random deviate of unit mean, using random1 as the source of uniform deviates.
fft.cpp - Fast Fourier transform Replaces data[1..2*n] by its discrete Fourier transform, if isign is input as 1; or replaces data[1..2*n] by nn times its inverse discrete Fourier transform, if isign is input as -1. data is a complex array of length nn, or equivalently, a real array of length nn, or, equivalently, a real array of length 2*nn. nn MUST be an integer power of 2 (this is not checked for!). If tehy don't, just pad your data with zeros up to the next power of two. The extended version uses inline method developed by Todd Veldhuizen and compares to a similar program from "Numerical Recipes". four1() implements a Cooley-Tukey FFT on N points, where N is a power of 2. "Numerical Recipes in C", 2nd Edition (Teukolsky,Press,Vetterling,Flannery). It serves as the basis for the inline version found below. All of the variable names have been preserved in the inline version for comparison purposes. Also, the line doubles are referred to in the corresponding template classes.
find root.cpp , butil.h - Finds the root of a function, f(x)=0, in an interval x1, x2, using Brent's method. The root, returned as zbrent, will be refined until is accuracy is tol. This is a robust program, which should be used for general purpose.
incomplete beta.cpp Returns the incomplete beta function, I_x(a,b), and the binomial probability distribution P. I_x(a,b)=(int_0^x t^(a-1)*(1-t)^(b-1)dt)/B(a,b), where B(a,b) is the betafunction B(a,b)=int_0^1 t^(a-1)*(1-t)^(b-1)dt = Gamma(a) * Gamma (b) / Gamma(a+b) The Probability distribution function is given by P=Sum_(j=k)^n (n j)p^j(1-p)^(n-j)=I(k,n-k+1)
incomplete gamma.cpp - Returns the imcomplete gamma functions P(a,x) and Q(a,x).
random.cpp Tests use of random double generators routine random. Generates random doubles between 0.0 and 1.0. Call with idum a negative double to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX approximates the largest floating value that is less than 1. The period of this routine is ~ 10^8.
fcomplex.cpp , fcomplex.h - Operations with complex doubles. For a and b complex, and x,y real: Absolute value of a: a.Cabs(); Conjugate of a: a.Conjg(); Square root of a: a.Csqrt(); Real part of a: a.r() Imaginary part of a: a.i(); Transforming x,y pair into complex: a=fcomplex(x,y); Adding a+b: a.Cadd(b) Subtracting a-b: a.Csub(b); Multiplying a*b: a.Cmul(b); Dividing a/b: a.Cdiv(b); multiplying by a real double a*x: a.RCmul(a)
binomial coeff.cpp - Returns the binomial coefficient (n k), 0<=k<=n, as a floating point double.
gamma function.cpp - Returns the value of Gamma(z)=int_0^infinity t^{z-1} e^{-t} dt for z>0
screen plot.cpp - For interactive terminal use. Produce a crude graph of the function fx over the prompted for interval x1, x2. Querey for another plot until user sginals satisfaction.
statistics.cpp Given an array of data in a file, this routine returns its mean, average deviation adev, standard deviation sdev, skewness skew, and kurtosis, kurt.
gasdedv.cpp - Returns a normally distributed deviate with zero mean and unit variance, using random1 as the source of uniform deviates.
legendre.cpp - Computes Legendre Polynomials P_l^m(x), where m and l are integers satisfying 0<= m <= l, while x lies in the range -1 <= x <= l.
deriv.cpp butil.h - Returns the derivative of a tabulated function func (array func), of equally spaced points with stepsize h using a five point differentation formula.
deriv.cpp , butil.h - Returns the derivative of a tabulated function.
interp.cpp , butil.h - Test interpolation routines. Enter number of points, and their values xa's and ya's. Calculate the interpolation at a desired Given arrays xa[1..n] and ya[1..n], and given a value of x, this routine returns a value y, and an error estimate dy. If P(x) is the polynomial of degree N-1 such that P(xa_i) = ya_i, i=1,..n, then the returned value y = P(x).
minimum search.cpp , butil.h - Finds the minimum of a function in an interval ax, bx, using the. Golden search method.

 Other
 ssort.cpp - A simple sort program. A real program should use C++ standard library algorithm 'sort'. test_sort.cpp , sort.h , butil.h - This overloaded class (a) organizes a matrix in order of increasing values - uses the quicksort algorithm, (b) constructs an index table. test_quick_sort.cpp , butil.h - Tests use of sorting a list of integers with the quick sort algorithm and the recursion method address book.cpp - Builds an address book from input screen. This program prints a calendar for a year. The user enters the year, after which the year and calendar are printed to the file calendar.txt. treedemo - Demonstrates binary search tree. stack - Demonstrates recursive high level procedures on stack.