Double integrals with variable upper limit

[niteOwl]

How can you compute
$F(k) = k\int^{\infty}_{0}dy\int^{y}_{0}dx f(kx,y)$

in C. I know about Python's scipy.integrate.dblquad function but it's just too slow. I have written some Cython code with a 2D gaussian quadrature function in C but it only takes doubles as limits. I think C doesn't have anything like Python's lambda, so how can this be done in C?

 

[pasmith]

How can you compute
$F(k) = k\int^{\infty}_{0}dy\int^{y}_{0}dx f(x,y)$

in C. I know about Python's scipy.integrate.dblquad function but it's just too slow. I have written some Cython code with a 2D gaussian quadrature function in C but it only takes doubles as limits. I think C doesn't have anything like Python's lambda, so how can this be done in C?

This is $k$ times a constant, which happens to be $\int_0^\infty \int_0^y f(x,y)\,dx\,dy$. You can obtain a rectangular, finite domain of integration by making the substitution $(x,y) = (\mathrm{arctanh}(s) \cos\theta, \mathrm{arctanh}(s) \sin \theta)$ to obtain $$\int_0^1 \int_{\frac14\pi}^{\frac12\pi} \frac{\mathrm{arctanh}(s)}{1 - s^2} f(\mathrm{arctanh}(s) \cos\theta, \mathrm{arctanh}(s) \sin \theta)\,d\theta\,ds.$$ Probably best to use a scheme of integration which doesn't require values of the integrand on $s = 1$.

 

[niteOwl]

This is $k$ times a constant

i forgot to say the integrand also depends on k, question is edited.

 

[chiro]

Hey niteOwl.

You could theoretically code a routine that loops through the one dimensional integrals and adds them up to achieve a double integral.

If you can integrate one "strip" (as it were) you just add up the strips in the double integral.

The accuracy will be determined with respect to the kind of function being integrated and if you know various derivative properties you could optimize it so that you only do so many computations per evaluation and that gets processed in some finite time interval. If you do this so many times within a couple of loops then you will be able to gauge the computational complexity of the integration routine.

If you want to optimize further then try learning floating point or GPU programming as an exercise.