How to optimise the computation of multiple integral for a specific integrand?

I am working on the problem of numerically integrating a function for , , .

To take the specific example of , the ambition would be to compute numerically

where denotes the (hyper)cube described by .

It is very much the case that the packaged functions within Matlab R2022a: integral, integral2 and integral3 do achieve this objective but because of the specific nature of q it would seem reasonable to think that execution times could be reduced. Note that I may be evaluated as an iterated integral as implemented in the Matlab functions.

In this application q takes as one of its arguments a vector of handles to the functions each of which must be called in the evalution of q. Initial experiments suggest that calling in advance the integrand q with a large 'meshgrid' of points in anticipated by the quadrature scheme, is far more efficient than evaluating the integrand q the large number of times (for a relatively small number of points) required by the packaged quadrature schemes.

Second, irrespective of the forms of each function in , the integrand q is very "well-behaved", i.e. no discontinuities, no singularities and no issues of scale. In fact, through much of it exhibits minimal functional variation but for in the vicinity of a small number N of distinct locations in the interior of .

My thinking therefore is that execution times should be improved if:

1. the quadrature rule is optimised for an integrand whose nature is both well-behaved and well-understood
2. the implementation is such that function values can be computed in large quantity in advance and/or stored for future iterations to reduce the quantity (rather than size) of function calls

If anyone in the community has expertise on such a problem I welcome ideas.

Answer 1

[Torsten]

Precalculate q over the hypercube on a predefined regular mesh and use "trapz" dimensionwise.

Look at the example

"Multiple numerical integrations" under

https://de.mathworks.com/help/matlab/ref/trapz.html

This can be generalized to n dimensions.

To test whether the spacing is sufficiently fine, I suggest you compare results with m and 2*m grid points in each coordinate direction together with Richardson extrapolation.

[RickyBoy]

Thank you @Torsten, this is a very interesting proposition. It is already showing much potential in my tests so far for . For the benefit of the wider discussion it may be worth commenting that, whilst the exact functional form of q cannot be known a priori due to it taking the arbitrary vector of functions as an argument, a multivariate Taylor expansion of q implies in itself (to the first-order) a critical value of mesh spacing. I will respond after further tests.