|
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)
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)
x1 = x0 + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
v1 = v0 + l1 / 6 + l2 / 3 + l3 / 3 + l4 / 6
f(t, x, v) = v and F(t, x, v) = -(b * v + k * x) / m
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
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.
|
|
|
|
