For Science and Engineering

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)	2^x
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[1] + 2*x[2] - 4, 2) + math.pow(2*x[1] + x[2] - 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[1], x[2], 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:
ODE System

* 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[1] = y[2]
   ydot[2] = (-c * y[2]) + (-k * y[1])
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[1]
   dy[i] = yy[2]
   t0 = t1
   t1 = t1 + dt
   y0[1] = yy[1]
   y0[2] = yy[2]
end

-- plot solution y and y'
p = plot.new(800, 600)
plot.set(p, "title", "Damped Oscillator")
plot.add(p, tm, y)
plot.add(p, tm, dy)
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 ax + by
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.