We have spent some time now^{[1][2][3]} looking at how we might numerically approximate the integrals of functions but have thus far completely ignored the problem of integrating functions with more than one argument, known as multivariate functions.
When solving integrals of multivariate functions mathematically we typically integrate over each argument in turn, treating the others as constants as we do so. At each step we remove one argument from the problem and so must eventually run out of them, at which point we will have solved it.
For example, we can integrate the square root of the sum of \(x\) and \(y\) with the steps
In listing 1 that helper function is
Program 1 demonstrates that this algorithm is indeed capable of approximating the integral of the square root of \(x\) plus \(y\) that we solved above by comparing its result to our exact solution.
Unfortunately, numerically integrating numerical integrals compounds their approximation errors, which you can see for yourself by setting
In generalising the bivariate integrator to a multivariate integrator we no longer know in advance how many dimensions we're going to be integrating over and so we must construct the nested integrals at run time. The
Finally, the
Program 3 demonstrates that this yields exactly the same result as program 1 for the vector argument version of our test function.
Whilst this certainly demonstrates how we can recursively integrate a multivariate function with a univariate integrator, it can hardly be considered a robust implementation given its hard coded use of
Note that if the
The
Here we're providing a default lower bound of all zeros if
Finally, the
This checks that the
This robust implementation of our recursive numerical integration algorithm is defined in NestedIntegral.js and is demonstrated by program 4.
You can see that this can create a multivariate integrator from any of our univariate integrators by changing
Unfortunately, recursively applying numerical integrators to numerical integrators gets more and more error prone as the number of dimensions increases and we must consequently increase the number of steps to maintain accuracy. This can get extremely costly for more than say three or four dimensions and so we must turn to alternative approaches for such integrals.
[2] The Tip Of The Romberg, www.thusspakeak.com, 2015.
[3] A Jump To The Left And Some Steps To The Right, www.thusspakeak.com, 2016.
When solving integrals of multivariate functions mathematically we typically integrate over each argument in turn, treating the others as constants as we do so. At each step we remove one argument from the problem and so must eventually run out of them, at which point we will have solved it.
For example, we can integrate the square root of the sum of \(x\) and \(y\) with the steps
\[
\begin{align*}
\int_{x_0}^{x_1} \int_{y_0}^{y_1} \sqrt{x + y} \; \mathrm{d}y \; \mathrm{d}x
&=\int_{x_0}^{x_1} \left(\int_{y_0}^{y_1} \left(x + y\right)^\tfrac12 \; \mathrm{d}y\right) \; \mathrm{d}x\\
&=\int_{x_0}^{x_1} \bigg[\tfrac23\left(x+y\right)^\tfrac32\bigg]_{y_0}^{y_1} \; \mathrm{d}x\\
&=\int_{x_0}^{x_1} \tfrac23\left(x+y_1\right)^\tfrac32  \tfrac23\left(x+y_0\right)^\tfrac32 \; \mathrm{d}x\\
&=\bigg[\tfrac4{15}\left(x+y_1\right)^\tfrac52  \tfrac4{15}\left(x+y_0\right)^\tfrac52\bigg]_{x_0}^{x_1}\\
&=\left(\tfrac4{15}\left(x_1+y_1\right)^\tfrac52  \tfrac4{15}\left(x_1+y_0\right)^\tfrac52\right)  \left(\tfrac4{15}\left(x_0+y_1\right)^\tfrac52  \tfrac4{15}\left(x_0+y_0\right)^\tfrac52\right)
\end{align*}
\]
It is consequently extremely tempting to approximate multivariate integrals by recursively applying univariate numerical integration algorithms to each argument in turn.
A Bivariate Integrator
We can illustrate the principle with an integrator for bivariate functions, being those that have two arguments. The trick is to implement a helper function that integrates over their second arguments for a fixed value of their first ones and then to integrate it over their first arguments.In listing 1 that helper function is
fx
which in turn defines another helper function fy
that fixes the value of the first argument of the function f
and then integrates it with our ak.rombergIntegral
integrator. The biIntegral
function then numerically approximates the integral of fx
to yield an approximation of the bivariate integral of f
.Listing 1: A Bivariate Integrator
function biIntegral(f, x0, x1, y0, y1) { function fx(x) { function fy(y) {return f(x, y);} return ak.rombergIntegral(fy)(y0, y1); } return ak.rombergIntegral(fx)(x0, x1); }
Program 1 demonstrates that this algorithm is indeed capable of approximating the integral of the square root of \(x\) plus \(y\) that we solved above by comparing its result to our exact solution.
Program 1: Using biIntegral



Unfortunately, numerically integrating numerical integrals compounds their approximation errors, which you can see for yourself by setting
x0
and y0
to zero in program 1 and noting that the integral fails to converge. At the cost of increased approximation error we can fix this by using a larger convergence threshold, as demonstrated by program 2.
Program 2: Fixing The Convergence Failure



A Multivariate Integrator
Rather than defining multivariate functions as having long lists of arguments, it is convenient to define them as taking a single vector argument, which we can represent with ourak.vector
type. Listing 2 shows how we might do so for the function that we've been using as an example so far.Listing 2: A Vector Argument Multivariate Function
function f(x) { return Math.sqrt(x.at(0)+x.at(1)); }
In generalising the bivariate integrator to a multivariate integrator we no longer know in advance how many dimensions we're going to be integrating over and so we must construct the nested integrals at run time. The
nested
function in listing 3 does so recursively, using an array x
as a buffer for the elements of the ak.vector
argument of the multivariate function f
which is constructed in the function returned in the last step. Note that the functions created at each step populate their elements of that array with their arguments and then integrate over the function created at the next step.Finally, the
multiIntegral
function integrates the outermost nested function over the last dimension.Listing 3: A Multivariate Integrator
function nested(f, x0, x1, x, n) { if(n===0) { return function(xn) { x[0] = xn; return f(ak.vector(x)); }; } var g = nested(f, x0, x1, x, n1); return function(xn) { x[n] = xn; return ak.rombergIntegral(g)(x0.at(n1), x1.at(n1)); }; } function multiIntegral(f, x0, x1) { var n = x0.dims(); var x = new Array(n); var g = nested(f, x0, x1, x, n1); return ak.rombergIntegral(g)(x0.at(n1), x1.at(n1)); }
Program 3 demonstrates that this yields exactly the same result as program 1 for the vector argument version of our test function.
Program 3: Using multiIntegral



Whilst this certainly demonstrates how we can recursively integrate a multivariate function with a univariate integrator, it can hardly be considered a robust implementation given its hard coded use of
ak.rombergIntegral
and its complete lack of error checking.
ak.nestedIntegral
Whilst we go about bringing the implementation of the algorithm up to scratch, we might as well put some thought into optimising it. For one thing, we can convert the vector bounds into arrays to slightly reduce the cost of accessing their elements at each step of the recursion. For another, we can move the construction of the integrators outside of the recursively defined functions that use them. Thenested
function defined in listing 4 exploits both of these optimisations whilst allowing the integrator
, its convergence threshold
and maximum number of steps
to be specified rather than hard coded. Apart from these minor improvements it has exactly the same structure as our prototypical implementation.Listing 4: An Improved Nested Integral
function nested(integrator, f, threshold, steps, x0, x1, x, n) { var g, integral; if(n===0) { return function(xn) { x[0] = xn; return f(ak.vector(x)); }; } g = nested(integrator, f, threshold, steps, x0, x1, x, n1); integral = integrator(g, threshold, steps); x0 = x0[n1]; x1 = x1[n1]; return function(xn) { x[n] = xn; return integral(x0, x1); }; }
Note that if the
integrator
is ak.trapeziumIntegral
, the threshold
will be interpreted as the trapezium width and the number of steps
will be ignored, so the nested
function is compatible with it as well as our other, more sophisticated, numerical integrators.The
multiIntegral
function given in listing 5 adds error checking to the arguments. In particular, that the function to be integrated is indeed a function, that the number of dimensions is greater than zero, that the number of steps is a number, that the bounds have the same number of dimensions, and that the integration algorithm produces a function.Listing 5: An Improved Multivariate Integral
function multiIntegral(integrator, f, threshold, steps, x0, x1) { var n = x0.dims(); var x, g, integral; if(ak.nativeType(f)!==ak.FUNCTION_T) { throw new Error('invalid function in ak.nestedIntegral'); } if(n===0) throw new Error('invalid lower bound in ak.nestedIntegral'); steps = ak.nativeType(steps)===ak.UNDEFINED_T ? 18 : ak.floor(Math.abs(steps)); if(isNaN(steps)) throw new Error('nonnumeric steps passed to ak.nestedIntegral'); steps = ak.ceil(steps/n); if(ak.nativeType(x1)===ak.UNDEFINED_T) { x1 = x0; x0 = ak.vector(n, 0); } else if(ak.type(x1)!==ak.VECTOR_T  x1.dims()!==n) { throw new Error('invalid upper bound in ak.nestedIntegral'); } x0 = x0.toArray(); x1 = x1.toArray(); x = new Array(n); g = nested(integrator, f, threshold, steps, x0, x1, x, n1); integral = integrator(g, threshold, steps); if(ak.nativeType(integral)!==ak.FUNCTION_T) { throw new Error('invalid integrator in ak.nestedIntegral'); } return integral(x0[n1], x1[n1]); }
Here we're providing a default lower bound of all zeros if
x1
is undefined and a default maximum number of steps
if it isn't specified. Furthermore, if each step of the univariate integrator
involves a subdivision into \(k\) intervals, as they typically will, then in one dimension we have for \(n_s\) steps
\[
n_f = k^{n_s}
\]
function evaluations, whereas in \(d\) dimensions we have
\[
n_f = \left(k^{n_s}\right)^d = k^{n_s \times d}
\]
and so we divide the maximum number of steps by the number of dimensions so that we have the same maximum number of function calls for any number of dimensions.Finally, the
ak.nestedIntegral
function given in listing 6 creates a function that uses this to perform the actual integration of multivariate functions.Listing 6: ak.nestedIntegral
function uniIntegral(integrator, f, threshold, steps, x0, x1) { var integral = integrator(f, threshold, steps); if(ak.nativeType(integral)!==ak.FUNCTION_T) { throw new Error('invalid integrator in ak.nestedIntegral'); } return integral(x0, x1); } function nestedIntegral(integrator, f, threshold, steps) { return function(x0, x1) { return ak.type(x0)===ak.VECTOR_T ? multiIntegral(integrator,f,threshold,steps,x0,x1) : uniIntegral(integrator,f,threshold,steps,x0,x1); }; } ak.nestedIntegral = function(integrator) { if(ak.nativeType(integrator)!==ak.FUNCTION_T) { throw new Error('invalid integrator in ak.nestedIntegral'); } return function(f, threshold, steps) { return nestedIntegral(integrator, f, threshold, steps); }; };
This checks that the
integrator
is itself a function and the function that it creates with nestedIntegral
is an integrator that switches between the multiIntegral
and uniIntegral
helper functions depending upon whether the lower bound of integration is a vector or not, with the latter simply integrating the function with the univariate integrator
.This robust implementation of our recursive numerical integration algorithm is defined in NestedIntegral.js and is demonstrated by program 4.
Program 4: Using ak.nestedIntegral



You can see that this can create a multivariate integrator from any of our univariate integrators by changing
uni
from ak.rombergIntegral
to ak.simpsonIntegral
, for example, and running the program again.Unfortunately, recursively applying numerical integrators to numerical integrators gets more and more error prone as the number of dimensions increases and we must consequently increase the number of steps to maintain accuracy. This can get extremely costly for more than say three or four dimensions and so we must turn to alternative approaches for such integrals.
References
[1] Trapezium Artistry, www.thusspakeak.com, 2015.[2] The Tip Of The Romberg, www.thusspakeak.com, 2015.
[3] A Jump To The Left And Some Steps To The Right, www.thusspakeak.com, 2016.
Leave a comment