I'm having issues obtaining the numerical values of various integrals similar to the one below: $$\int_0 ^L \int_0^L\log(n(x,y))\sin{\left( i\dfrac{\pi}{L}x\right)}\sin{\left( \mu\dfrac{\pi}{L}x\right)}\cos{\left( j\dfrac{\pi}{L}y\right)}\cos{\left( \nu\dfrac{\pi}{L}y\right)}dxdy$$ Here, $n(x,y)$ is a user-defined function that has a step-like behaviour, meaning that it takes the value $n_1$ in certain region and $n_2$ in the rest of the square region defined by $x,y \in [0,L]$.
$i,j,\mu,\nu$ are integers that range from 1 to another user-defined value. I have tried to solve the integrals letting them go up to 10.
I've tried solving those integrals using scipy.integrate.dblquad and mpmath.quad. However, I cannot trust the results. Let me explain why. Both integration routines also calculate an error (uncertainty) estimation of the integration result. This tends to be the case when $i,j,\mu$ or $\nu$ have high-values, between 7 and 10. In some of those cases, I'm getting relative uncertainties bigger than 1. To get over this, I've tried to subdivide the interval $[0,L]$ in smaller steps, integrate the functions in those smaller regions and add the results. In each of the integrals over a smaller domain I can define an accuracy parameter that bounds above the relative uncertainty of the result. Say, for example, that each of those integrals is precise up to 3 significant digits.
However, if I fix the values of $i,j,\mu,\nu$ and compute one of those integrals, but changing the number of divisions I'm taking of $[0,L]$, the integral results vary wildly, even though it was always the case that each step was supposedly accurate up to 3 significan digits. The same happens with the integrals that did not give any warning when being evaluated over the whole domain. When I then evaluate them in pieces and add the results, the values change, sometimes by orders of magnitude. I don't know how to tackle this problem.
I've tried to look for specific algorithms to integrate oscillating functions, but so far I have found nothing that I can use. Any hint or help would be greatly appreciated.
Thank you in advance
Just a Thought:
You could perhaps change the order of the double integral
$$\int_0 ^L \sin{\left( i\dfrac{\pi}{L}x\right)}\sin{\left( \mu\dfrac{\pi}{L}x\right)}\cos{\left( j\dfrac{\pi}{L}x\right)}\cos{\left( \nu\dfrac{\pi}{L}x\right)}\int_0^L\log(n(x,y))\, dy\,dx\tag{1}$$
and then use an identity to give up to 8 simpler integrals, which can each be numerically integrated (by parts) over $x$.
$$\sin \left(i\frac{\pi}{L}x\right) \sin \left(\mu \frac{\pi}{L}x\right) \cos \left(j\frac{\pi}{L}x\right)\cos \left(\nu\frac{\pi}{L}x\right)\equiv \frac{1}{8} \left( -\cos \left(\frac{\pi i x}{L}+\frac{\pi j x}{L}+\frac{\pi \mu x}{L}+\frac{\pi \nu x}{L}\right) -\cos \left(\frac{\pi i x}{L}-\frac{\pi j x}{L}+\frac{\pi \mu x}{L}+\frac{\pi \nu x}{L}\right) +\cos \left(\frac{\pi i x}{L}+\frac{\pi j x}{L}-\frac{\pi \mu x}{L}+\frac{\pi \nu x}{L}\right) +\cos \left(\frac{\pi i x}{L}-\frac{\pi j x}{L}-\frac{\pi \mu x}{L}+\frac{\pi \nu x}{L}\right) -\cos \left(\frac{\pi i x}{L}+\frac{\pi j x}{L}+\frac{\pi \mu x}{L}-\frac{\pi \nu x}{L}\right) -\cos \left(\frac{\pi i x}{L}-\frac{\pi j x}{L}+\frac{\pi \mu x}{L}-\frac{\pi \nu x}{L}\right) +\cos \left(\frac{\pi i x}{L}+\frac{\pi j x}{L}-\frac{\pi \mu x}{L}-\frac{\pi \nu x}{L}\right) +\cos \left(\frac{\pi i x}{L}-\frac{\pi j x}{L}-\frac{\pi \mu x}{L}-\frac{\pi \nu x}{L}\right)\right)$$
(Identity generated using Mathematica TrigReduce[] function to save me time)
Update after correction to question:
Question is now
$$\int_0 ^L \cos{\left( j\dfrac{\pi}{L}y\right)}\cos{\left( \nu\dfrac{\pi}{L}y\right)}\int_0^L \sin{\left( i\dfrac{\pi}{L}x\right)}\sin{\left( \mu\dfrac{\pi}{L}x\right)}\log(n(x,y))\, dx\,dy\tag{2}$$
Analogous identities are:
$$\cos{\left( j\dfrac{\pi}{L}y\right)}\cos{\left( \nu\dfrac{\pi}{L}y\right)}\equiv \frac{1}{2} \left(\cos \left(\frac{\pi j y}{L}+\frac{\pi \nu y}{L}\right)+\cos \left(\frac{\pi j y}{L}-\frac{\pi \nu y}{L}\right)\right)$$
$$\sin{\left( i\dfrac{\pi}{L}x\right)}\sin{\left( \mu\dfrac{\pi}{L}x\right)}\equiv \frac{1}{2} \left(\cos \left(\frac{\pi i x}{L}-\frac{\pi \mu x}{L}\right)-\cos \left(\frac{\pi i x}{L}+\frac{\pi \mu x}{L}\right)\right)$$
I think this will be of help in terms of integrating numerically; at least the method has the attribute of separating out the component frequencies to help with the diagnosis of what is going wrong with the original numerical integral.