Math. Js: An Advanced Mathematics Library For JavaScript

v
is the speed of the object, and 
m
is the 
object’s mass. The “pitch” angle 
γ
is the angle that the object’s velocity vector makes with the 
surface of the moon. 
T
is the thrust in the forward direction, 
0
g
is standard gravity (9.80665 
m/s
2
), 
sp
I
is the specific impulse of the ascent stage engine, and 
μ
is the product of the gravita-
tional constant and the moon’s mass. In this simplified model, we consider motion in just two 
dimensions and exclude relativistic effects. 
During an ascent, propellant is used most efficiently by following a trajectory known as a gravity 
turn. The lunar module lifts off vertically, and then briefly uses the reaction control system to 
pitch over slightly. During the remainder of the ascent, the lunar module accelerates only in the 
forward direction, while gravity causes it to “turn” toward the horizon. The goal is to reach the 
target altitude (
target
r
r
=
) at the correct orbital speed (
/
orbit
target
v v
r
μ
=
=
), while at the same 
moment having the lunar module rotated exactly horizontal (
0
γ
=
), resulting in a circular orbit. 
Several variables can be adjusted to alter the final orbit of the lunar module, including initial 
pitch-over angle and engine burn time. 
Our problem statement, then, is this: find the combination of initial pitch-over angle and engine 
burn time to place the lunar module in a specified orbit. 
The backbone of our tool will be the following JavaScript function that solves a system of ODEs 
using Euler’s method. The function could be readily modified to use a more accurate method, 
such as RK4, but we leave this as an exercise for the reader. 
function ndsolve(f, x0, nsteps, tmax) {
var n = f.size()[0]; // Number of variables 
var x = x0.clone(); // Current values of variables 
var dxdt = []; // Temporary variable to hold time-derivatives 
var result = []; // Contains entire solution 
var dt = math.divide(tmax, nsteps); // Time step 
for(var i=0; i// Compute derivatives 
for(var j=0; jdxdt[j] = f.get([j]).apply(null, x.toArray());
} 
// Euler’s method to compute next time step 
for(var j=0; jx.set([j], math.add(x.get([j]), math.multiply(dxdt[j], dt))); 
} 
result.push(x.clone()); 
} 
return math.matrix(result); 
} 
The function accepts a math.js matrix f of functions representing the time-derivatives of the 
equations of motion. x0 is the set of initial conditions, steps is the number of steps for the Euler 
method, and tmax is the final time of the simulation. The evolution of the dependent variables is 
computed from time 0 to tmax, and the entire solution is returned as a math.js matrix. (Although 
time
is used as the independent variable in this example, this does not necessarily need to be the 
case.) 
The observant reader will note the use of the functions math.add, math.multiply, and math.divide 
in places where one would typically find the more usual operators +, *, and /. This is necessary 
because we have designed the function ndsolve to work with all data types supported by math.js, 
including physical quantities containing units. Likewise, the math.js methods get and set are used 
25
January/February 2018
www.computer.org/cise

---

THE RISE OF JAVASCRIPT 
instead of the array indexing operator [ ], to retrieve and assign, respectively, elements in a 
math.js matrix. 
Our solution will make use of the math.js expression parser. The syntax of the expression parser 
is intended to be more familiar to mathematicians, and is well suited to solve our present prob-
lem. The parser is invoked using the following: 
var parser = math.parser(); 
parser.eval(input); 
where input is a string containing one or more mathematical expressions. To use ndsolve in the 
context of the math.js expression parser, it must first be imported using the following statement: 
math.import({ndsolve:ndsolve}); 
The expression parser accepts any valid mathematical expression, which might contain units, 
complex numbers, or matrices. We will take advantage of ndsolve’s polymorphic nature and de-
fine our parameters and initial conditions using physical units by passing the following parame-
ters and initial conditions into the expression parser: 
G = 6.67408e-11 m^3 kg^-1 s^-2 # Gravitational constant 
mbody = 7.342e22 kg # Mass of Moon 
mu = G * mbody 
nSteps = 200 # Number of steps in Euler’s Method 
tfinal = 300 s # Burn time 
T = 3500 lbf # Engine thrust (fixed) 
g0 = 9.80665 m/s^2 # Standard gravity 
isp = 311 s # Engine isp 
gamma0 = 89.2 deg # Initial pitch-over angle 
r0 = 1738.1 km # Equatorial radius of Moon 
v0 = 1 m/s # Initial velocity 
phi0 = 0 deg # Reference angle 
m0 = 10300 lbm # Mass of Lunar Module and propellent 
We will choose an initial pitch-over angle of 89.2 degrees (0.8 degrees from vertical) and a burn 
time of 300 seconds. The system has a singularity at 
0
0
v
=
, so we choose a small, nonzero ini-
tial velocity to initiate Euler’s method. 
We will represent our system of ODEs by defining several functions, one for each dependent 
variable. The expression parser allows functions to be defined using the syntax f(x) =