Numerical sign problem software?

[mollwollfumble]

I'm an applied mathematician with lots of experience in classical physics but none in quantum physics.
What's the simplest piece of code (eg C++ on Cygwin) that illustrates the https://en.wikipedia.org/wiki/Numerical_sign_problem
in quantum physics?

 

[mollwollfumble]

Have posted question on 3 forums but no replies. Perhaps I'm asking the wrong question. To simplify answer just choose a,b,c,... From this list.
A. There is no such code because you can't get a probability from complex numbers.
B. Ask the author of a journal paper because it's a very advanced topic.
C. Use a full blown lattice qcd package.
D. Can't run on PC because it's too slow.
E. Not in c++. (I'll accept any language).
F. Write from scratch using ideas from "lattice qcd for novices" or similar.
G. The question is ambiguous.

 

[A. Neumaier]

Feynman path integrals (over $e^{iS/\hbar}$) in quantum mechanics (and even more in QFT) are so highly oscillatory that they can hardly be evaluated numerically. The remedy is to work in the Euclidean formulation where (due to the time integration in the action $S$) the constant $i/\hbar$ responsible for the wild oscillations becomes a harmless negative real constant, and analytically continue the results back to real time. In particular, this is used in all of lattice gauge theory, which is feasible numerically only in its Euclidean version.

 

[atyy]

See p8 of http://www.ncts.ncku.edu.tw/phys/qis/091219/files/07Naoki_Kawashima.pdf

 

[A. Neumaier]

see [above post] - its not simple to code, though. But as a toy example you can easily encode the trapeziodal rule for approximating $\int_{-a}^a e^{i x^2} dx$ in any language of your choice, and observe its misbehavior as $a$ increases!

 

[mollwollfumble]

> See p8 of http://www.ncts.ncku.edu.tw/phys/qis/091219/files/07Naoki_Kawashima.pdf
LOL at that. I can see from that how $(3-2)^{10}$ becomes difficult to accurately evaluate using the binomial theorem.
Triangular lattices would be interesting. In my computational fluid dynamics we use finite element or finite volume on lattices where a lattice cell can be any shape; that makes it easy to locally refine it in areas of interest.

see [above post] - its not simple to code, though. But as a toy example you can easily encode the trapeziodal rule for approximating $\int_{-a}^a e^{i x^2} dx$ in any language of your choice, and observe its misbehavior as $a$ increases!

Thanks. Will try that. Also will try for $a$ tends to infinity, because I've been looking at equations like this from a pure maths perspective.

 

[A. Neumaier]

I can see from that how $(3-2)^{10}$ becomes difficult to accurately evaluate using the binomial theorem.

You can also try to evaluate $e^{-100}$ using the power series expansion of $e^x$. (But these examples have nothing to do with quantum physics.)

 

[A. Neumaier]

Also will try for a tends to infinity

The natural setting for such limits of oscillating integrals in finite dimensions is the Henstock integral.

 

[mollwollfumble]

I really like the toy example $\int_{-a}^a e^{i x^2} dx$. The trapezoidal rule blows up spectacularly, as does Simpsons rule. Also, it's just difficult enough to be non-trivial. But with a transformation $\int_{-\infty}^\infty e^{i x^2} dx = \int_0^\infty \cos(x)/\sqrt{x}\,dx$, integration schemes like the Newton-Coates schemes or the Legendre-Gauss quadrature so beloved of Finite Element analysis rapidly produce a result to high accuracy, when coupled with an iterative extrapolation $I_n=\int_{n\pi}^{(n+1)\pi} \cos(x)/\sqrt{x}\,dx$, $I_{n+1}=I_n(I_n-2I_{n-1})/(2I_n+5I_{n-1})$ for large enough $n$.

I need to try an example in multidimensions next.

 

[A. Neumaier]

it's just difficult enough to be non-trivial

This still corresponds to a noninteracting situation, hence is physically trivial (and can be solved exactly for $a=\infty$). For a nontrivial interacting version consider $\int e^{ip(x)}dx$ with an even degree polynomial $p(x)$ such as $p(x)=x^2+ \epsilon x^4$ and multivariate generalizations.

 

[mollwollfumble]

The natural setting for such limits of oscillating integrals in finite dimensions is the Henstock integral.

I'm still reading about the Henstock integral. I haven't yet convinced myself that it is a) consistent in the sense that the answer is independent of the gauge, b) necessary for QM.

On the pure maths side, I looked at historical systems of infinite numbers. There are a dozen or so, including transseries, Hahn series, Cantor cardinals, Cantor ordinals and Conway's surreals. These are all different, but related. The two systems of relevance to QM are those of Du Bois Reymond's infinitary calculus (1877) and Robinson's hyperreals (1960s). Du Bois Reymond's infinitary calculus ties in with the QM cancellation of infinities and renormalization/regularization. Robinson's hyperreals tie in with the QM concepts of "spooky action at a distance" and the "many-world interpretation".

Starting from Du Bois Reymond, (whose theory generates limits for smoothly diverging sequences) I further developed the theory to include limits of oscillating divergent sequences (and by extension divergent series and Riemann integrals). Hence my interest in the numerical sign problem, I ought to be able to generate unique answers even when the oscillations are divergent. My past experience also includes quadrature and low-discrepancy sequences (for speeding up Monte-Carlo).

This still corresponds to a noninteracting situation, hence is physically trivial (and can be solved exactly for $a=\infty$). For a nontrivial interacting version consider $\int e^{ip(x)}dx$ with an even degree polynomial $p(x)$ such as $p(x)=x^2+ \epsilon x^4$ and multivariate generalizations.

I'll follow that advice. I'll try the multivariate generalization.
$\int_{-\infty}^{\infty} e^{i((x_1-1)^2+(x_2-x_1)^2+(x_3-x_2)^2+(x_4-x_3)^2+(x_5-x_4)^2+(x_6-x_5)^2+(x_7-x_6)^2+(1-x_7)^2 +\sum_{i=1}^7p(x_i))}dx_1dx_2dx_3dx_4dx_5dx_6dx_7$
This may take some time, but if I can manage that then I'll be well on the way to handling any integral numerical sign problem. If I succeed then the next hurdle would be in tackling cases which involve a ratio of complex integrals.

 

[A. Neumaier]

If I succeed then the next hurdle

If you succeed, you'd tell us something about the techniques you employ; they are interesting even for the toy examples.
But if you want to get really serious about applications in quantum mechanics, you should consider the true Feynman path integral, not the toy versions we have been discussing so far! See, e.g., the following paper:

Albeverio, S., & Mazzucchi, S. (2006). The time-dependent quartic oscillator—a Feynman path integral approach. Journal of Functional Analysis, 238(2), 471-488.

See also the 2008 book Mathematical Theory of Feynman Path Integrals by Albeverio et al..

transseries, Hahn series, Cantor cardinals, Cantor ordinals and Conway's surreals.

I don't see how these relate to oscillating integrals. The generalization of the number concept most useful in the present context seems to be nonstandard analysis; see ref [10] in the above paper.

 

[A. Neumaier]

when coupled with an iterative extrapolation

This works in the particular case since the oscillations are equidistant, a property peculiar to the noninteracting situation. A general purpose method cannot rely on such accidents.

 

[mollwollfumble]

Trying that again with equations, and with nonlinearity added to the end.

Here's how the numerical solution of an oscillating divergent integral in QM might work. Equal intervals are used but are not an essential part of the solution method. The first is a toy problem, far too simplistic.
Solve $$\int_{-\infty}^\infty e^{i(x_1^2+...+x_7^2)}dx_1...dx_7$$

The real part equals $$2\int_0^\infty \cos(x_1^2+...+x_7^2)\,dx_1...dx_7$$

Convert to hyperspherical coordinates, that equals $$2S_7\int_0^\infty r^6\cos(r^2)\,dr$$
where $$S_7= 16\pi^3/15$$ is the surface ofr a 7-ball with unit radius.

Substitute $$y=r^2$$ to get the integral component
$$\int_0^\infty y^{2.5}\cos(y)dy$$

To solve this divergent oscillating integral numerically, start by integrating by parts on a finite domain.
$$\int_0^x y^{2.5}\cos(y)dy=x^{2.5}\sin(x)+2.5x^{1.5}\cos(x)-\frac{15}{4}x^{0.5}\sin(x)+\frac{15}{8}\int_0^xy^{-0.5}\sin(y)dy$$

Here's where the study of infinite numbers comes in. It's OK to redefine the infinite limit of a pure oscillating function to be equal to zero. Examples of pure oscillating functions include $$f(x)\cos(g(x))$$ for a range of non-oscillating functions f and g including non-negative powers, finite polynomials, exponential and logarithmic functions. So
$$\lim_{x\to\infty}(x^{2.5}\sin(x)+2.5x^{1.5}\cos(x)-\frac{15}{4}x^{0.5}\sin(x))=0$$
and
$$\int_0^\infty y^{2.5}\cos(y)\,dy=\frac{15}{8}\int_0^\infty y^{-0.5}\sin(y)\,dy$$

The integral on the right hand side is oscillating convergent, so can be directly solved numerically. For faster numerical convergence, do the direct integral only from 0 to π and carry out two more applications of integration by parts to get the term $$\int_{\pi}^{\infty}y^{-2.5}\sin(y)\,dy$$

This can, of course, all be considered moot because the original integral has an analytic solution. But my hope is that a similar numerical method can be used for oscillating integrals when no analytic solution is available.

A word of warning. In general,
$$\lim_{x\to\infty}\frac{\int^xA}{\int^xB}\ne\frac{\int^\infty A}{\int^\infty B}$$
When calculatating expected values, the left hand side must be used not the right hand side.

The nonlinear oscillator is not much more difficult. Start with
$$\int_{-\infty}^\infty e^{i\sum_{n=1}^7p(x_n)}dx_1...dx_7$$
This is separable
$$\int_{-\infty}^\infty e^{i\sum_{n=1}^7p(x_n)}dx_1...dx_7=\prod_{n=1}^7{\int_{-\infty}^\infty e^{ip(x_n)}dx_n}$$

Separate out the real and imaginary component of the one-dimensional integral then, after solving each, multiply the seven complex numbers. For the real component (the imaginary component is treated in an identical way) this is
$$2\int_0^\infty \cos(p(x_n))dx_n$$
Change the variable to $$y=p(x_n)$$
$$dx_n=dy/p'^{-1}$$

When $$p(x_n)=x^2+ax^4$$ that gives
$$p'^{-1}=\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}$$
so we're integrating
$$\int_0^\infty \cos(y)/\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}dy$$

We might as well split that into two parts, direct numerical integration for
$$\int_0^\pi \cos(y)/\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}dy$$
and integrate by parts twice (to suppress numerical oscillations) for
$$\int_\pi^\infty \cos(y)/\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}dy$$

$$\frac{d}{dy}\frac{\sin(y)}{p'^{-1}}=\frac{\cos(y)}{p'^{-1}}+\sin(y)A$$
where
$$A=\frac{4-6\sqrt{1+4ay}}{((2/a)(1+4ay)(\sqrt{1+4ay}-1)^{3/2}}$$

Integrating by parts once more
$$\frac{d}{dy}(-\cos(y)A)=\sin(y)A-\cos(y)\frac{dA}{dy}$$
where
$$\frac{dA}{dy}=\frac{\sqrt{1+4ay}(84ay+33)-120ay-30}{\sqrt{2}\sqrt{1+4ay}(((1+4ay)^{3/2}-1-4ay)/a)^{5/2}}$$

Using our definition of the limit of pure oscillatories at infinity.
$$\lim_{y\to\infty}\sin(y)/p'^{-1}=\lim_{y\to\infty}\cos(y)A=0$$

Then just integrate numerically in 1-D out to whichever multiple of π is sufficient to give the accuracy desired, for two decimal places, 6π may be large enough to use instead of ∞.

In the general nonlinear oscillator where $$p(x_n)\ne x^2+ax^4$$ the procedure is the same, but A and dA/dy are obtained by numerically differentiating 1/p'^-1 rather than analytically differentiating it.

 

[A. Neumaier]

The nonlinear oscillator is not much more difficult.

Thanks for providing the details.

Yes, the separable case is not very difficult since it factors into univariate integrals. Also the general quadratic case is not difficult since after a linear transformation it becomes separable.

But in the multivariate interacting case, the exponent has the form $i(x^TKx/2m+V(x))$ where $m>0$ is a mass parameter, the potential $V(x)$ is separable and bounded below but $K$, the stiffness matrix of the oscillating system, is symmetric positive definite and non-diagonal. The tridiagonal case is already nontrivial since there is no simple transformation to separable form.

Thus you need to extend your current techniques.

 

[A. Neumaier]

A word of warning.

In practice, expectation values are most often calculated via derivatives of the free energy $F=i\log Z$, where $Z$ denotes the integral, after approximating the free energy by an appropriate analytic expression.

 

[mollwollfumble]

Got it. I said before that I wanted to be able to numerically integrate the path integral
$$I=\int_{-\infty}^\infty e^{i((x_1-1)^2+(x_2-x_1)^2+(x_3-x_2)^2+(x_4-x_3)^2+(x_5-x_4)^2+(x_6-x_5)^2+(x_7-x_6)^2+(1-x_7)^2+\sum_{n=1}^7p(x_n))}\,dx_1dx_2...dx_7$$
Write this as
$$I=\int_{-\infty}^\infty e^{iq}\,dx_1dx_2...dx_7$$
The antisymmertic parts of $q$ don't contribute to $I$ so can be omitted, leaving
$$I=\int_{-\infty}^\infty e^{i(2+\sum_{n=1}^7(2x_n^2+p(x_n)))}\,dx_1dx_2...dx_7$$
This is separable, giving
$$I=e^{2i}(\int_{-\infty}^\infty e^{i(2x^2+p(x))}\,dx)^7$$
Only an integral in 1-D is needed. Change variable using $y=2x^2+p(x)$ to give
$$\int_{-\infty}^\infty e^{i(2x^2+p(x))}\,dx=2\int_0^\infty e^{iy}(dx/dy)dy$$
Here for simple $p(x)$ the formula for $dx/dy$ can be found analytically and for messy $p(x)$ the formula for $dx/dy$ can be found numerically. This then resolves to the case of the nonlinear oscillator described above and the same methods are used. In a nutshell, keep integrating by parts until the oscillations are suppressed enough to get an accurate answer using standard quadrature methods. The differentials for the integral by parts $$d^2x/dy^2,\,d^3x/dy^3$$ etc. can be evaluated either analytically or numerically, depending on the original $p(x)$.

This method also goes a long way towards solving the more general integral
$$ \int_{-\infty}^\infty f(\mathbf{x})e^{ig(\mathbf{x})}\,d\mathbf{x}$$
when $f(\mathbf{x})$ is separable and $g(\mathbf{x})$ does not contain symmetric cross correlations such as $x^2y^2$.
Here, it's $f(x)(dx/dy)$ that is differentiated again for each iteration of integration by parts.

Thanks for providing the details.

Yes, the separable case is not very difficult since it factors into univariate integrals. Also the general quadratic case is not difficult since after a linear transformation it becomes separable.

But in the multivariate interacting case, the exponent has the form i(xTKx/2m+V(x))i(xTKx/2m+V(x))i(x^TKx/2m+V(x)) where m>0m>0m>0 is a mass parameter, the potential V(x)V(x)V(x) is separable and bounded below but KKK, the stiffness matrix of the oscillating system, is symmetric positive definite and non-diagonal. The tridiagonal case is already nontrivial since there is no simple transformation to separable form.

Thus you need to extend your current techniques.

In practice, expectation values are most often calculated via derivatives of the free energy F=ilog⁡Z" F=ilogZF=ilog⁡Z, where Z denotes the integral, after approximating the free energy by an appropriate analytic expression.

Excellent. I'll see what I can do on that basis, as said before this may take some time. (sorry about loss of equations in quote, I don't know why).

 

[A. Neumaier]

The antisymmetric parts of $q$ don't contribute to $I$

There is no antisymmetric part in a non-diagonal quadratic form. $\int dx^n e^{ix^TKx}$ is proportional to $(\det K)^{-1/2}$, but your formula with $p(x)=0$ gives something different!

By the way, at PF, use two hashes # in place of each $. Morever, when quoting formulas, you need to edit out the artifacts created by the system! (You can edit it even days afterwards).

 

[mollwollfumble]

Oops. Thank you very very much for correcting my error. (It's not by accident that "fumble" is part of my username).

All along, from the very beginning, I've had three ultimate goals in mind. Three general-purpose numerical methods for evaluating the non-separable asymmetric n-dimensional integral
$$\int_{-\infty}^\infty f(\mathbf{x})e^{ig(\mathbf{x})}\,d\mathbf{x}$$
(And by extension, expected values using the ratio of two such integrals).
All three start the same way. Find the minimum of $g(\mathbf{x})$. Call the xector at that minimum $\mathbf{x}_0$. This minimum can be found either analytically, or by some variation of the conjugate gradient method, Polak-Ribiere if analytical derivatives of $g(\mathbf{x})$ are available or Powell's method if they are not.
From here on evaluate
$$I=\int_{-\infty}^\infty f(\mathbf{x})e^{ig(\mathbf{x}-\mathbf{x}_0)}\,d\mathbf{x}$$

Method 1, suitable when $g(\mathbf{x}-\mathbf{x}_0)$ is dominated by the quadratic term $(\mathbf{x}-\mathbf{x}_0)^TK(\mathbf{x}-\mathbf{x}_0)$. This has high speed and low accuracy.

The conjugate gradient method generates, as a byproduct, the principal axes of the ellipsoid in n dimensions around the minimum point. Using that, or the principal axes found analytically, numerically integrate $I$ along each of the n principal axes in turn, using integration by parts to reduce oscillations. Multiply the resulting n integrals to get the final answer $I$.

Method 2, suitable when $g(\mathbf{x}-\mathbf{x}_0)$ has a significant quadratic term component. This has intermediate speed and intermediatre accuracy.
Numerically integrate $I$ along $n(n+1)/2$ one dimensional vectors given in cartestian coordinates by $x_1-x_0, x_2-x_0, x_2-x_1, x_3-x_0, x_3-x_1, x_3-x_2, x_4-x_0,...,x_n-x_{n-1}$, using integration by parts in each case to reduce oscillations.

From this derive an approximate matrix $K$ using $x_ix_j=(1/2)((x_i-x_0)^2+(x_j-x_0)^2-(x_i-x_j)^2)$ for off-diagonal terms $i\ne j$. Evaluate $I$ from $\sqrt{\det(K)}$.

Method 3, suitable for nasty $f$ or $g$. This is slow and has an accuracy that ranges from poor to good depending on the number of vectors chosen.

Define on the surface of a unit sphere in n-D, the endpoints of $m>>n$ unit vectors using an algorithm that spaces these endpoints $\mathbf{x}_1$ to $\mathbf{x}_m$ as far apart from each other as possible. One way to do this would be to define a grid similar to that used to define the nodes of a geodesic dome. A second way would be use use additive quasirandom numbers (see wikipedia for algorithm on a hypercube and then map to the hypersphere surface using the Sag-Szekeres transformation).

Then use a higher-dimensional extension of the method of Veronoi polyhedra (or a faster approximation) to assign a weight $w_1,...,w_m$ to each vector, endpoints further away from their neighbours are allocated a greater weight.

Numerically integrate $I$ along the $m$ one-dimensional vectors, using integration by parts in each case to reduce oscillations.

Evaluate $I$ as the nth root of the product of the $m$ integrations, each weighted using the assigned weights $w_1,...,w_m$.

A correction for kurtosis is necessary, there's no correction when the 1-D distribution $\int_{-\infty}^\infty f(x)e^{-g(x)}\,dx$ is Gaussian (zero excess kurtosis). When the excess kurtosis is -1.2 (uniform distribution) the correction is $S_n/(n2^n)$ and when the excess Kurtosis is 3 (Laplace distribution) the correction is $S_n(n-1)!/2^n$ where $S_n$ is the surface of an n-ball.

 

[A. Neumaier]

Method 1, suitable when $g(\mathbf{x}-\mathbf{x}_0)$ is dominated by the quadratic term $(\mathbf{x}-\mathbf{x}_0)^TK(\mathbf{x}-\mathbf{x}_0)$. This has high speed and low accuracy.

The conjugate gradient method generates, as a byproduct, the principal axes of the ellipsoid in n dimensions around the minimum point. Using that, or the principal axes found analytically, numerically integrate $I$ along each of the n principal axes in turn, using integration by parts to reduce oscillations. Multiply the resulting n integrals to get the final answer $I$.

This is a variant of the saddle point method used to find approximations to the integral. it is much used in quantum mechanics. The assumption is not satisfied when there are quartic terms, which always dominate the quadratic term for large $x$. It is often used nevertheless as an asymptotic series.

Method 2, suitable when $g(\mathbf{x}-\mathbf{x}_0)$ has a significant quadratic term component. This has intermediate speed and intermediatre accuracy.
Numerically integrate $I$ along $n(n+1)/2$ one dimensional vectors given in cartestian coordinates by $x_1-x_0, x_2-x_0, x_2-x_1, x_3-x_0, x_3-x_1, x_3-x_2, x_4-x_0,...,x_n-x_{n-1}$, using integration by parts in each case to reduce oscillations.

From this derive an approximate matrix $K$ using $x_ix_j=(1/2)((x_i-x_0)^2+(x_j-x_0)^2-(x_ix_j)^2$ for off-diagonal terms $i\ne j$. Evaluate $I$ from $\sqrt{\det(K)}$.

I don't see why this should give an approximation of the integral. Please explain for $n=2$.

Method 3, suitable for nasty $f$ or $g$. This is slow and has an accuracy that ranges from poor to good depending on the number of vectors chosen.

Define on the surface of a unit sphere in n-D, the endpoints of $m>>n$ unit vectors using an algorithm that spaces these endpoints $\mathbf{x}_1$ to $\mathbf{x}_m$ as far apart from each other as possible. One way to do this would be to define a grid similar to that used to define the nodes of a geodesic dome. A second way would be use use additive quasirandom numbers (see wikipedia for algorithm on a hypercube and then map to the hypersphere surface using the Sag-Szekeres transformation).

This is best done via spherical $t$-designs and related constructions; see my paper
Measures of strength 2e, and e-optimal designs, Sankhya 54 (1992), 299-309.
and its citations.

Then use a higher-dimensional extension of the method of Veronoi polyhedra (or a faster approximation) to assign a weight $w_1,...,w_m$ to each vector, endpoints further away from their neighbours are allocated a greater weight.

Numerically integrate $I$ along the $m$ one-dimensional vectors, using integration by parts in each case to reduce oscillations.

Evaluate $I$ as the nth root of the product of the $m$ integrations, each weighted using the assigned weights $w_1,...,w_m$.

A correction for kurtosis is necessary, there's no correction when the 1-D distribution $\int_{-\infty}^\infty f(x)e^{-g(x)}\,dx$ is Gaussian (zero excess kurtosis). When the excess kurtosis is -1.2 (uniform distribution) the correction is $S_n/(n2^n)$ and when the excess Kurtosis is 3 (Laplace distribution) the correction is $S_n(n-1)!/2^n$ where $S_n$ is the surface of an n-ball.

Again I don't see how you can replace the integration over $R^n$ by $m$ univariate integrations.

 

[mollwollfumble]

I shouldn't have mentioned kurtosis, using that is nowhere near good enough, and I already have a better correction. Note earlier my integral:
$$2S_7\int_0^\infty r^6\cos(r^2)dr$$
In n dimensions for the third method, the integral needed is the n-1th moment of the absolute value (not the 4th moment used for kurtosis).
$$\frac{nS_n}{m}\int_{-\infty}^\infty |r|^{n-1}f(r)e^{ig(r)}dr$$
For higher order moments, more iterations of integration by parts are needed to suppress the oscillations.

 

[mollwollfumble]

This is a variant of the saddle point method used to find approximations to the integral. it is much used in quantum mechanics. The assumption is not satisfied when there are quartic terms, which always dominate the quadratic term for large $x$. It is often used nevertheless as an asymptotic series.I don't see why this should give an approximation of the integral. Please explain for $n=2$.This is best done via spherical $t$-designs and related constructions; see my paper
Measures of strength 2e, and e-optimal designs, Sankhya 54 (1992), 299-309.
and its citations.Again I don't see how you can replace the integration over $R^n$ by $m$ univariate integrations.

Saddle point method. I've heard of that. I'll look it up. That had to do with limiting oscillation by only considering trajectories in the complex domain that are connected to the minimum, or do I have that wrong?

I'm exceedingly interested in your paper on spherical t-designs.

As for getting an n-dimensional integral from m 1-dimensional integrals, I haven't worked out the details yet, or checked for errors, but there have to be two ways, one based on the weighted arithmetic mean and the other based on the weighted harmonic mean.

 

[mollwollfumble]

I've done it!

I'm now convinced that I have an algorithm that can evaluate oscillatory integrals in n dimensions
$$\int_{-\infty}^\infty\,g(\mathbf{x})\cos(f(\mathbf{x}))\,d\mathbf{x}$$
as fast as it can calculate the equivalent non-oscillatory integrals
$$\int_{-\infty}^\infty\,g(\mathbf{x})e^{-f(\mathbf{x})}\,d\mathbf{x}$$
Replacing cos by sin in order to compute the imaginary part doesn't slow the algorithm down.

A slight problem so far is that it doesn't evaluate the non-oscillatory integral particularly quickly, because unlike for example the vegas algorithm it doesn't have adaptive gridding, yet.

I've tested this using my previously mentioned 7-D path integral:
$$I=\int_{-\infty}^\infty e^{i((x_1-1)^2+(x_2-x_1)^2+(x_3-x_2)^2+(x_4-x_3)^2+(x_5-x_4)^2+(x_6-x_5)^2+(x_7-x_6)^2+(1-x_7)^2+\sum_{n=1}^7 p(x,n))}\,dx_1dx_2...dx_7$$
I began by subtracting off the value of $\mathbf{x}$ for which $f$ is a minimum, but this is probably an unnecessary operation.

Define on the surface of a unit sphere in n-D, the endpoints of $m>>n$ unit vectors using an algorithm that spaces these endpoints $\mathbf{x}_1$ to $\mathbf{x}_m$ as far apart from each other as possible. To do this I'm using additive quasirandom numbers. Define starting vector
$$\mathbf{v}_1=(\sqrt{2},\sqrt{3}\sqrt{5},\sqrt{7},\sqrt{11},\sqrt{13},\sqrt{17},...)$$
Map it back onto the n-D hyperecube using $\mathbf{v}_1:=\mod(\mathbf{v}_1,1)$
and iterate
$$\mathbf{v}_j=\mod(\mathbf{v}_{j-1}+\mathbf{v}_1,1)$$
Map onto the n-D interval $[-1,1]$ using $\mathbf{x}=2\mathbf{v}-1$ and reject those values of $\mathbf{x}$ outside the unit n-sphere to get $|\mathbf{x}|<=1$ . For numerical reasons I also reject values of $\mathbf{x}$ whose norm is too close to zero to get $|\mathbf{x}|>=0.001$. Then I map the result onto the surface of the unit n-sphere using $\mathbf{x}:=\mathbf{x}/|\mathbf{x}|$. This gives a set of $m$ uniformly distributed vectors on the surface of an n-sphere.

Quasirandom numbers have two slight advantages over polyhedron methods. One is that with a polyhedron there's a risk of a false result due to alignment of the symmetries of the problem with one of the symmetries of the polyhedron. The second is that quasirandom numbers work almost equally well with any arbitrary number of vectors whereas with polyhedron methods it's much more difficult to achieve that.

I don't use a weighted sum, just a normnal Riemann sum over these vectors. If adaptive gridding is used then a weighted sum would be needed. In mathematics this is a sum of 1-D integrals. Let's write this using symbol $y$ for a radial distance along the direction defined by $\mathbf{x}_j$, so $g(\mathbf{x})=g(y\mathbf{x}_j)$ and $f(\mathbf{x})=f(y\mathbf{x}_j)$. $S_n$ is the surface area of a unit n-ball. For example, $S_7=16\pi^3/15$. For arbitrarily large $n$,
$$S_{n-1}=2\pi^{n/2}/\Gamma(n/2)$$

The desired answer is given by
$$\frac{S_n}{m}\sum_{j=1}^m{\int_0^\infty{y^{n-1}g(y)\cos(f(y))\,dy}}$$
I split this into three parts, a near zero part $y<\pi$ that is evaluated directly, a constant part at $y=\pi$ that comes from integration by parts of the integral for $y\ge\pi$, and an integral for $y>\pi$ that comes from the integration by parts. For my 7-D path integral, integrating by parts four times suppresses oscillations beautifully.

Then
$$I=I_1+I_2+I_3$$
$$I=\int_0^\infty{y^{n-1}g(y)\cos(f(y))\,dy}$$
$$I_1=\int_0^{f^{-1}(\pi)}y^{n-1}g(y)\cos(f(y))\,dy$$
Values of $f^{-1}()$ are found using Newton's method.
$$I_2=-y^{(0)}\sin(\pi)-y^{(1)}\cos(\pi)+y^{(2)}\sin(\pi)+y^{(3)}\cos(\pi)=y^{(1)}-y^{(3)}$$
This doesn't work for every possible choice of $g()$, but does work for most choices of $g()$ that are non-oscillatory. The third component of the integral is
$$I_3=\int_\pi^\infty y^{(4)}\cos(f)df$$

For best numerics, differentiate $f()$ analytically and define $y^\prime=1/(df/dy)$. Then iteratively define
$$y^{(0)}=y^{n-1}g(y)y^\prime$$
$$y^{(k+1)}=y^\prime\frac{d}{dy}y^{(k)}$$
$y^{(0)}$ is evaluated analytically and $y^{(k)}$ for $k\ge 1$ can be found analytically or numerically. I find $y^{(k)}$ for $k\ge 1$ numerically by first expanding the iterative formula using divided differences, which gives $y^{(k)}$ directly from $y^{(0)}$. For example
$$y^{(1)}=\frac{y^\prime_0}{\Delta y}(y^{(0)}_{0.5}-y^{(0)}_{-0.5})$$
$$y^{(3)}=\frac{y^\prime_0}{(\Delta y)^3}(y^{(0)}_{1.5}y^\prime_1 y^\prime_{0.5}-y^{(0)}_{0.5}(y^\prime_1 y^\prime_{0.5}+y^\prime_{0.5} y^\prime_{0}+y^\prime_0 y^\prime_{-0.5})+y^{(0)}_{-0.5}(y^\prime_{0.5} y^\prime_{0}+y^\prime_0 y^\prime_{-0.5}+y^\prime_{-0.5}y^\prime_{-1})-y^{(0)}_{-1.5}y^\prime_{-1}y^\prime_{-0.5})$$
The formula for $y^{(4)}$ is twice as long as the formula for $y^{(3)}$, reply here if you want the full version.
The choice of $\Delta y$ is a balance, too large and errors appear due to function nonlinearity, too small and roundoff errors are amplified. I use $\Delta y=\pi/200$ in the 7-D integral.
The subscript notation is shorthand for the following. From given real number $f$, which comes from the integration point in $I_3$, define $y_0=f^{-1}(f)$, then for any subscript $a$, $y_a=y_0+a\Delta y$, $y^\prime_a$ is $1/(df/dy)$ where $df/dy$ is evaluated analytically at $y=y_a$. $y^{(0)}_a=y_a^{n-1}g(y_a)y^\prime_a$. I haven't yet said how to analytically evaluate $df/dy$, from the definition of $f(\mathbf{x})$ define the vector $df/d\mathbf{x}$ and substitute $\mathbf{x}=y\mathbf{x}_j$. If not used correctly, calculating $f^{-1}()$ could also be time consuming, but the algorithm can be and is set so that this is only called once per vector in evaluating $I_1$ and called once per integration point per vector in evaluating $I_3$.

To evaluate $I_3$, replace $\infty$ with some multiple of $\pi$. $3\pi$ gives errors of order $10^{-5}$ on my 7-D path integral. I suggest $4\pi$, so use three intervals.
$$I_3=\sum_{j=1}^3\int_{j\pi}^{(j+1)\pi} y^{(4)}\cos(f)df$$
Across each interval for $I_1$ and $I_3$ I use a single application of the 10 point Lagrange-Gauss quadrature (see Press "Numerical Recipes" for details) which has a maximum error of about $10^{-9}$. So each evaluation of $I_3$ requires $10\times 3=30$ integration points. Overall, I get
$$\int_{-\infty}^\infty\cos\left((x_1-1)^2+(1-x_7)^2+\sum_{k=1}^6(x_{k+1}-x_k)^2+\sum_{k=1}^7(x_k^2+x_k^4)\right)dx_1dx_2dx_3dx_4dx_5dx_6dx_7=-0.627$$

For the imaginary part, the three components of the integral are
$$I=I_1+I_2+I_3$$
$$I_1=\int_0^{f^{-1}(\pi/2)}y^{n-1}\sin(f(y))\,dy$$
$$I_2=y^{(1)}-y^{(3)}$$
$$I_3=\sum_{j=1}^3\int_{(j-1/2)\pi}^{(j+1/2)\pi} y^{(4)}\sin(f)df$$
The integration limits are chosen to help suppress oscillations because, for all $j$,
$$ \int_{j\pi}^{(j+1)\pi} \cos(f)df=\int_{(j-1/2)\pi}^{(j+1/2)\pi} \sin(f)df=0$$
The amount of oscillation suppression obtained by choosing these limits is roughly equivalent to one step of integration by parts.

The most numerically time consuming part of the algorithm is $I_3$, and this only accounts for 4% of the total value in my test case, so evaluate it using fewer vectors than used for $I_1$ and $I_2$. One tenth of the number of vectors works well.

---------------

An outline for a tentative possible method for grid adaption is as follows. (Here I use the words "point" and "vector" interchangeably because a point on the surface of an n-sphere defines a normalised vector from the origin to that point). From the cumulative distribution of evaluated integrals $I_1+I_2$ (not $I_3$) along each vector, attract points to where the slope of the cumulative distribution is largest or repel points from where the slope is smallest. The process of attraction or repuslion is as follows. Given an attractor as a normalised vector $\mathbf{a}$, and using $\mathbf{x}=2\mathbf{v}-1$ as above for random points inside the n-sphere (no need to normalise yet, but it also works for other point generation methods where the vectors are normalized), map $\mathbf{x}$ using $\mathbf{x}:=\mathbf{x}+\alpha(\mathbf{x}\cdot\mathbf{a})\mathbf{a}$ followed by normalization. When $0<\alpha$ this is an attaction towards both $\mathbf{a}$ and its antipodal point $-\mathbf{a}$ of maximum strength $(1+\alpha)^{n-1}$ in n-D. When $-1<\alpha<0$ this is a repulsion away from $\mathbf{a}$ and its antipodal point of maximum strength $(1+\alpha)^{1-n}$.

This is a very long way from being a working algorithm. Issues still to be considered include how to avoid distorting the cumulative distribution once mapping is enacted, how many attraction/repulsion points should be used and with what values of $\alpha$, and how to keep track of the strength of attraction/repulsion for every vector (for $I_3$ as well) along which the integral is calculated. Once that's done, the weighting for the Riemann sum is trivial to calculate, it's identical to the net strength of repulsion.

 

[mollwollfumble]

You were asking earlier where the study of systems of infinite numbers comes in. It’s difficult to describe the purpose of subroutine intgcosf() for the evaluation of the definite integral g(x)cos(f(x)) without first digressing into a bit of pure mathematics.

While doing a bit of pure mathematics I found that every function that I could construct could be decomposed into a pure oscillatory part and a non-oscillatory part. The pure oscillatory part could be periodic or random or chaotic, the non-oscillatory part is none of these. This means that a function that is oscillating divergent can sometimes have the pure oscillatory component subtracted off to yield a convergent result. That’s what this subroutine is doing, using integration by parts to suppress the pure oscillatory component of an oscillating divergent integral to get a real convergent result. We see this process elsewhere in applied mathematics, in the evaluation of “asymptotic integrals”, and in the evaluation of “divergent series”.

Some results, a mix of theoretical and numerical evaluations of the integral of g(x) cos(x) over the interval 0 to infinity, give the following table of integrals. You won't find this table anywhere else on the web. Accuracy of the numerical evaluations is not guaranteed, particularly in the last decimal place.

 

[mfb]

Most of your integrals are not even well-defined, so it is not surprising they are not listed anywhere.
The others are easy to get with the usual tools. Sometimes even analytically.

This forum is not the right place for software development, I closed the thread.