Math. Js: An Advanced Mathematics Library For JavaScript

expression
. 
Math.js automatically compiles these into native JavaScript functions for better performance. 
drdt(r, v, m, phi, gamma) = v sin(gamma) 
dvdt(r, v, m, phi, gamma) = -mu / r^2 sin(gamma) + T / m 
dmdt(r, v, m, phi, gamma) = -T/g0/isp 
dphidt(r, v, m, phi, gamma) = v/r cos(gamma) * rad 
dgammadt(r, v, m, phi, gamma) =
(1/r * (v - mu / (r v)) * cos(gamma)) * rad 
Finally, we group our functions and initial conditions into matrices and pass them to ndsolve: 
result = ndsolve([drdt, dvdt, dmdt, dphidt, dgammadt],
[r0, v0, m0, phi0, gamma0], nSteps, tfinal) 
After the call to ndsolve, we can use the eval function to extract the results in a format that can 
be plotted using a third-party library such as NVD3. We will make use of three additional matrix 
functions and operators: concat, subset [...] and transpose ‘. In the JavaScript “one-liner” below, 
26
January/February 2018
www.computer.org/cise

---

COMPUTING IN SCIENCE & ENGINEERING 
we use the eval function to transpose two subsets of the result matrix, strip the units, and concat-
enate those. We then convert the result to a JavaScript array, and use the map function to bundle 
the results into a format that can be plotted by NVD3: 
var xy_result = parser 
.eval("concat((result[:,4]' - phi0) * r0/rad/km, " +
"(result[:,1]' - r0) / km, 1)' ") 
.toArray() 
.map(function(e) { return {x: e[0], y: e[1]}; }); 
As Figure 2 shows, our initial parameters of pitch-over angle and burn time did not result in a 
stable orbit; in fact, this scenario ends very badly for our brave astronauts. To achieve a stable 
orbit, we will take this example one step further by creating an objective function that ap-
proaches zero as the trajectory approaches the target orbit. We will then pass this objective func-
tion to an optimization routine. 
Figure 2. Lunar module trajectory before optimization. Initial pitch-over angle is 89.2 degrees and 
burn time is 300 seconds. 
We first import the optimization routine, which will not be listed here for the sake of brevity. 
The function signature is optimize(fxn, arr). The routine finds a local minimum of the objective 
function fxn, using the guess values arr. The resulting values are returned in a new array. 
math.import({optimize:optimize}); 
We then enter our target values and define the objective function using the expression parser. In 
the example below, objFxn uses the last point in the provided trajectory y to calculate the objec-
tive function: 
rTarget = r0 + 33 km # 33 km above surface 
gammaTarget = 0 deg # Ensures a circular orbit 
vTarget(r) = sqrt(mu/r) # Ensures a circular orbit 
objFxn(y) = (y[nSteps, 1] - rTarget) ^ 2 / km^2 +
(y[nSteps, 2] - vTarget(y[nSteps, 1])) ^ 2 / (m/s)^2 +
(y[nSteps, 5] - gammaTarget) ^ 2 / deg^2 
Next, we define a wrapper function that accepts a given set of input values x, computes a trajec-
tory, and returns the resulting value of the objective function. Because the optimize function was 
not designed to work with math.js types, we manually add or strip units from guess values as 
necessary: 
27
January/February 2018
www.computer.org/cise

---

THE RISE OF JAVASCRIPT 
ndsolveWrapper(x) = objFxn(ndsolve([drdt, dvdt, dmdt, dphidt,
dgammadt], [r0, v0, m0, phi0, x[1] * deg], nSteps, x[2] * s)) 
We then define our guess values and make the call to optimize: 
xGuess = [89.2, 300] 
xOpt = optimize(ndsolveWrapper, xGuess.toArray()) 
The optimized parameters can be extracted via the get function, and the resulting trajectory is 
shown in Figure 3. 
parser.get("xOpt"); // [89.5039862828086, 422.00063165607446] 
Figure 3. Optimized lunar module trajectory. Initial pitch-over angle is 89.5 degrees and burn time is 
422 seconds. 
Our entire solution involved physical quantities with units. Math.js automatically converts units 
from one system to another as required and checks each operation for dimensional consistency. 
In the event of an invalid operation, such as attempting to evaluate 1 m/s + 2 kg, math.js throws 
an error. These errors might go unnoticed when performing calculations by hand. Students, in 
particular, are taught to take great care when working with units to avoid errors. Despite rigorous 
training and care, errors in unit conversion have been the cause of several high-profile incidents 
in recent history. Math.js greatly reduces the risk of unit-conversion error. 
In solving our problem, we designed ndsolve to have the highest possible generality. By extend-
ing this example into a full application, a developer might, for instance, create an interactive 
online tool that scientists, mathematicians, engineers, or students can use to solve any system of 
first-order ODEs. The developer has the burden of writing and importing the ndsolve and opti-
mize functions, and associated plotting functions—but this will, in turn, leave the user free to 
solve any problem using the familiar mathematical syntax of math.js, without getting bogged 
down in the details of the underlying programming language. 
BENCHMARKS 
JavaScript performance has improved tremendously in recent years, but does it come anywhere 
near optimized math applications such as Matlab, Octave, or Python?
We performed benchmarks comparing several JavaScript libraries, including Math.js, with Oc-
tave.
30
Octave is an open-source math application comparable to Matlab. Octave can be consid-
ered “state of the art” in terms of speed because it is written in C++ and optimized for numerical 
28
January/February 2018
www.computer.org/cise

---

COMPUTING IN SCIENCE & ENGINEERING 
operations. The goal of this benchmark is to get an indication of the order of magnitude differ-
ence in performance, not to come up with an in-depth comparison. 
We chose matrix operations as a benchmark, because these are typically the most computation-
ally demanding and thus have a large impact on overall performance. The benchmark measures 
the duration of four basic matrix operations: adding two matrices (A+A), multiplying two matri-
ces (A*A), transposing a matrix (A’), and calculating the determinant of a matrix (det(A)). The 
particular contents of matrix A are assumed to have no effect on the outcome of the benchmark; 
in other words, the particular algorithms chosen to perform the matrix operations do not assume 
or take advantage of any special properties of matrix A. 
The selected operations were repeated at least 100 times for a 25 × 25 square matrix, and the av-
erage duration is listed in Table 2. 
Table 2. Math.js benchmarks for selected matrix operations compared to other JavaScript libraries. 
Average duration in microseconds is shown, with smaller numbers representing better 
performance.