Adaptive numerical integration of a univariate vector integrand

Question

Background & Problem formulation

I'm trying to write a simple program in C++ that performs adaptive numerical integration of vector valued integrands (in one variable), i.e.

$$\int_a^b \bar{f}(x) \ dx $$

where $\bar{f} = \{f_1, f_2, ..., f_n\}$.

I know that there are good libraries that can do this (for example GSL or Cubature) but I cannot use any of these because they are under GNU GPL. On the other hand I think these libraries are a bit to involved for my simple problems. But if there is library that I can use I will probably use it!

What I need is something simple like an adaptive Gauss-Kronrod or Gauss-Lobatto solver, but Newton-Cotes probably suffices.

I have been using a Gauss-Lobatto algorithm before that is inspired by this code and ideally I would just modify that one to be able to handle vector integrands.

Before doing so I tried to "vectorize" a simpler algorithm, namely an adaptive mid-point rule algorithm (i.e. is a rectangular Newton-Cotes). While doing so I realized that it is very difficult to "vectorize" the algorithm and maintain tail-recursion. And I noticed significant loss of performance (~30%) when removing the tail-recursion in this simple algorithm. So I imagine the impact can be even bigger in an adaptive Gauss-Lobatto algorithm.

So to summarize, here are some of my questions:

Q1

When writing my own adaptive algorithm should I avoid using recursion if I cannot maintain the tail-recursion?

Q2

Is there a library I can use that solves my problem and is under a less restrictive license than GNU GPL? (LGPL should be fine).

Thanks in advance!

Answer 1

I think you are trying to do something inherently impossible (or at least inefficient): how can you be "adaptive" if you want to integrate different functions at the same time?

The essence of adaptive quadrature is that your algorithm tries to detect the hard parts of the integrand and subdivides the integration interval accordingly. How can the algorithm do this efficiently if you ask it to do this for different functions at the same time?

If you are looking for a less restricting library, why not use the mother of all quadrature libraries: Quadpack by Piessens et al.? It's written in FORTRAN77 but you can link to it in C++.

Answer 2

Since your problem can be expressed as a set of independent ordinary differential equations with $x$ being the time variable, you could use the Boost Odeint library. You should find in easy to modify the examples given in the tutorial for the purposes of your problem.