I need to integrate complex, oscillatory function that is only known at equidistant grid points. Simpson's not accurate enough, what are my choices? I have a complex function and I need to integrate its real part (which oscillates highly).
I only know the function values at equidistant grid points in $[a, b]$. 
I have many such points, enough that I would normally expect a highly accurate result from trapezoid or simpson's rule ... but both schemes are only accurate up to 4th or 5th digit. This is accurate, but not accurate enough for my purposes.
What can I do to improve the accuracy? Is there some other method I can try, or adapt the above ones? 
Note again: I only know the function values at equidistant grid points. So I can't use Gaussian quadrature, nor composite Simpson's rule.
I am willing to trade computational time for more accuracy. Is there anything I can do?
 A: Newton-Cotes integration should be your friend.
A: Numerical integration of highly oscillating functions is a hard problem. There are accurate numerical integration schemes available, but they only work well in special cases. For example if the function is specified to be of the form $f(x)\exp\left[i u g(x)\right]$ and $f(x)$, $g(x)$ and $u$ are known or can be extracted accurately from your data.
I have obtained good results with an alternative non-quadrature approach to numerical integration, but this is still work in progress so it's as of yet unpublished. It can be used to tackle highly oscillating integrands and integration over high dimensional spaces, which are both regarded to be extremely hard problems. 
My method is based on Ramanujan's master theorem which is usually used to get to exact expressions for integrals. But, in principle, one can also attempt to apply this theorem to tackle numerical integration problems. According to the theorem, for a function $f(x)$ with series expansion of the form:
$$f(x) = \sum_{k=0}^{\infty}(-1)^k\frac{c_k}{k!}x^k$$
we have:
$$\int_0^{\infty} f(x)dx = c_{-1}$$
This is a special case of the expression:
$$\int_0^{\infty} x^s f(x)dx = s! c_{-(s+1)}$$
for $s=0$. Now, this theorem comes with conditions such as a limit on the growth of the modulus of the  analytically continued $c_s$ in the imaginary direction to be less than $\exp(\pi y)$. Without these conditions, you could wonder how to deal with the freedom to add a function proportional to $\sin(\pi n)$ to $c_n$, which leaves the $c_n$ invariant. Bu this is discussed in detail in the literature on Ramanujan's master theorem.
The obstacle one then faces is with extrapolating series expansion coefficients $c_n$ where $n$ an integer larger or equal to zero to $n = -1$. This has to be done extremely  accurately, otherwise there is no gain of this method relative to the conventional numerical algorithms.
In your case you would first need to extract accurate series expansion coefficients for your function around one point and then transform the integral to one from zero to infinity and then obtain the series expansion around zero of the new integrand. 
One can then, in principle, apply a suitable extrapolation method to the $c_n$, but this won't yield accurate results. If you go to higher and higher orders and include more and move values for $c_n$ for larger and larger $n$, the result may improve a bit at first but will soon become worse and worse. Such problems do have solutions, albeit using non-rigorous methods developed in theoretical physics.
The first step to tackle the problem is to write the result of applying the extrapolation of the coefficients as a perturbation of a trivial action that can be performed exactly. For the latter one can choose the identity operator. We thus consider the deformed extrapolation operator:
$$E(g) = \left(1 + g \Delta\right)^{-1}$$
where $\Delta$ is the difference operator defined as $\Delta c_n = c_{n+1} - c_n$. We then have that $E(g=1)$ applied to $c_0$ yields $c_{-1}$. But this is an infinite formal expression in powers of $\Delta$ for fractional or negative $k$. We have deformed this operator, the parameter $g$ interpolates between the trivial identity operator and the conventional extrapolation operator. This gives us a perturbation series in powers of $g$, which we want to evaluate at $g = 1$.
Directly substituting $g = 1$ in the series is obviously not going to bring us anywhere, as that's just $E(g=1)$ which gives us bad results. Instead, we can apply resummation methods developed in theoretical physics to extract extremely accurate results from (usually) divergent perturbation series. There are quite few different methods to choose from.
A method that sometimes works very well is the so-called "principle of minimum sensitivity". This works by introducing an auxiliary parameter in the perturbation series. It's then best to start with the original problem at hand and then introduce a parameter $p$ such that the result doesn't depend on $p$. Let's consider the operator:
$$E(g, k, p) = \left[1 + g \Delta + p g(1-g)\Delta^2\right]^{-k}\tag{1}$$
The parameter $k$ will be set to 1 at the end. The idea is then that because the exact result does not depend on $p$, that in cases where the extrapolation is very accurate that the dependence on $p$ will be small. The extrapolation will, of course, be accurate for small $k$. But the accuracy will also depend on $p$. If we were to plot the result for fixed $k$ as a function of $p$, we would see a large relatively flat region. The result is then the most accurate inside that region as a more accurate result implies that the result should depend less on $p$.
Suppose we make $k$ larger so that the result becomes less accurate for generic values for $p$. The approximation will then not have become uniformly worse for all values of $p$. For $p$ in certain regions the approximation will have dropped less in accuracy compared to $p$ in other regions. In such more accurate regions, the plot will still be very flat. So, what we'll see is that the large flat region for small $k$ will have become smaller and we now only get an accurate result if we choose $p$ inside the smaller flat region.
We  can therefore locate the flat region by setting the derivative of the result w.r.t. $p$ of the result to zero. It's advisable to start with small $k$ and gradually increase $k$ to 1 instead of jumping ahead and start at $k =1$. This way one can eliminate spurious results. 
Now, this will yield a polynomial equation for $p$. How do we then choose the best value? We look at the modulus of the second derivative and choose that solution that has the smallest modulus for the second derivative. 
You can improve the result by including the next best result, or if the result is complex, the complex conjugate of the value you've used. Let's write the result as a function of $p$ as $R(p)$. The zeroes of $R'(p)$ are $p_n$ for $n$ positive integers up to the degree of the polynomial, and we have:
$$\left|R''(p_1)\right| \leq \left|R''(p_2)\right|\leq\cdots$$
Then we can consider the function:
$$\widetilde{R}(a,p) = \frac{R(p) + a R\left(\frac{p_2}{p_1}p\right)}{1+a}$$
Then $\widetilde{R}(a,p)$ being an average of different values taken by $R(p)$ is also an estimate of the integral. The derivative of this function is zero at $p = p_1$ for general $a$. We can thus choose $a$ such that the second derivative will also be zero. This often yields a much better result than the original result. It can also be easily generalized to  getting to a function that has many derivatives set to zero by using many solutions.
You then pick a solution obtained this way and see how things change if you increase $k$ to 1. If the solution smoothly evolves all the way to $k =1$, then you'll  likely have found a very accurate solution. If you find that the solution seems to have a singularity as a function of $k$, then it's best to modify the term involving $p$ in expression (1).
It is also possible to work with more parameters, but in practice it become enormously more complex to do this in a systematic way. 
