Programming Environment for Numerical Computing

Math Functions

With Comet, the math functions are mapped to global functions (e.g.: user can use cos or math.cos).

Example:

Copy the following script in the editor...

-- Math
cls()
a = math.cos(math.pi/4)
print("a = ", a)

... and click Run (or press F12 on Windows and Linux)

Below is a summary of the Lua math functions:

 abs acos asin atan atan2 ceil cos cosh deg exp floor fmod frexp huge ldexp log log10 max min modf pi pow rad random randomseed sin sinh sqrt tanh tan

And a summary of the Comet global math functions:

 abs acos asin atan atan2 ceil cos cosh deg exp floor fmod frexp huge ldexp log log10 max min modf pi pow rad random randomseed sin sinh sqrt tanh tan

NB: Lua gives the Napierian logarithm the name log and the decimal logarithm is named log10, as in C language.
The Comet Math Console uses ln for Napierian logarithm and log for decimal logarithm.

Extended Math Functions

Special math functions are added to Comet (math namespace), including and extending the Lua math functions. Summary of math functions and constants:

Functions:

 math.exp2(x) 2x math.logb(x) log base 2 of x math.cbrt(x) cubic root math.hypot(x,y) sqrt(x²+y²) math.erf(x) error function math.erfc(x) complementary error function math.lgamma(x) ln(gamma(x)) math.tgamma(x) gamma(x) math.round(x) nearest integer, rounding math.isinf(x) number is infinite ? math.isnan(x) not a number ? math.isnormal(x) number is normal ? math.asinh(x) hyperbolic arc sine math.acosh(x) hyperbolic arc cosine math.atanh(x) hyperbolic arc tangent math.gauss(x,b,c) G(x) = exp(-(x - b)² / 2c²) math.lorentz(x,b,c) L(x) = (c / ((x - b)² + c²))

Constants (Universal constants in international units (SI)):

 math.q Electron charge (in C) math.me Electron mass (in kg) math.kb Boltzmann constant (J/K) math.h Planck constant (Js) math.c Speed of Light in vacuum (m/s) math.na Avogadro constant (1/mole)

Mathematical Optimization

Comet integrates a mathematical optimization module including minimization of real function of n variables.
The Comet minimization function uses the Hooke and Jeeves algorithm which do not require the jacobian to be evaluated.

iters = optim.minimize(func, pars, maxiters, tol, rho)

func: name of the function of n variables to be minimized.
func defined as func(x) where:
x: vector containing the variables
pars: table containing the initial values (will be updated to the calculated values)
maxiters: maximum number of iterations (optional)
tol: the algorithm tolerance (optional)
rho: the algorithm parameter, between 0 and 1 (optional). Decrease rho to improve speed and increase it for better convergence.

xr = math.root(func,a,b,iters,tol)

Find the root of function in the [a,b] interval, with the given maximum number of iterations (optional) and tolerance (optional)

Examples:

-- Optimization
function booth(x)
return (math.pow(x + 2*x - 4, 2) + math.pow(2*x + x - 5, 2))
end

x = { -5, 5 }
iters = optim.minimize("booth", x, 100, 1e-6)
cls()
io.write(string.format("\n x = [%g %g] iters = %d \n", x, x, iters))

function func(x)
return x^2 - x - 1
end

-- expected: x ~ -0.618034
x = optim.root("func", -5, 5)

print(x, func(x))

Ordinary Differential Equations

With the Comet ODE solver can integrate differential system given:
the system functions, the initial values and the independent variable value (x, t, ...).
It is not necessary to provide the Jacobian which is approximated by the solver.

The ODE solver:

 y = ode.solve(func, y0, t0, t1, tol) Integrates ODE system, where: * func: name of the ODE system function. func defined as func(t, y, ydot) where: ydot vector updated with respect to y and t: * y0 table containing the initial values * t0 value for y0 * t1 value to integrate for * tol the solver tolerance (optional)

Example:

Copy the following script in the editor and click Run (or press F12 on Windows and Linux)

-- Damped Oscillator: y'' + c y' + k y = 0

local c = 0.25
local k = 1

-- Oscillator function
function func(t, y, ydot)
ydot = y
ydot = (-c * y) + (-k * y)
end

-- Solve with 0.1 seconds as interval
y0 = {2, 0}
t0 = 0
dt = 0.1
t1 = t0 + dt
tol = 1e-3
tm = {}
y = {}
dy = {}
for i = 1,100,1 do
yy = ode.solve("func", y0, t0, t1, tol)
tm[i] = t1
y[i] = yy
dy[i] = yy
t0 = t1
t1 = t1 + dt
y0 = yy
y0 = yy
end

-- plot solution y and y'
p = plot.new(800, 600)
plot.set(p, "title", "Damped Oscillator")
plot.set(p, "xlabel", "time (s)")
plot.set(p, "ylabel", "amplitude")
plot.set(p, 1, "legend", "y(t)")
plot.set(p, 1, "style", "-")
plot.set(p, 1, "color", "red")
plot.set(p, 2, "legend", "y'(t)")
plot.set(p, 2, "style", "-")
plot.set(p, 2, "color", "blue")
plot.update(p)

Numerical Array Functions

Functions to handle numerical arrays that can be called through the linalg namespace.

Below is a summary of the linalg functions:

 x = linalg.array(n,a) create array of size n and initialize its values to a x = linalg.zeros(n) create n-size array and initialize its values to zero x = linalg.ones(n) create n-size array and initialize its values to one linalg.swap(x,y) swap two arrays linalg.copy(x,y) copy array x to y s = linalg.dot(x,y) scalar product z = linalg.add(x,y,a,b) return a*x + b*y d = linalg.mult(a,b,m,n) multiply a matrix (with m columns) with a matrix (with n columns) y = linalg.get(x,is,ie) return x subarray from index is to ie z = linalg.cat(x,y) concatenate x and y arrays B = linalg.transpose(A,m) transpose a matrix with m columns s = linalg.format(A,m) format an array with m columns and return a string s = linalg.print(A,m) print an array with m columns x = linalg.rand(n,rmin,rmax) return an array of random values between rmin and rmax i = linalg.imin(x) return index of the min value in x i = linalg.imax(x) return index of the max value in x s = linalg.sum(x) return the array sum d = linalg.norm(x) return the array norm-2 iA = linalg.inv(A) calculate the inverse of a square matrix d = linalg.det(A) calculate the determinant of a square matrix x = linalg.solve(A,b,mt,x0,tol,iters) solve the linear system Ax = b with: mt = "sparse" if A is a sparse matrix; x0 the initial guess; tol a given tolerance; iters the maximum number of iterations. x = linalg.tridiag(al,ad,au,b) solve the tridiagonal linear system Ax = b with: Al the lower diagonal; Ad the main diagonal and Au the upper diagonal.