I'm currently trying to perform a triple integration using C++. I sucessfully used the Matlab function integral3 which yields the correct result (I know that because it fits the analytical result in test cases). Now to do it on C++, I used the midpoint and the trapezoidal rule. For the midpoint rule we have :
$$\int_{a}^{b} \int_{c}^{d} \int_{e}^{f} \phi(x,y,z) dz dx dy = (b-a)(d-c)(f-e)\sum_{i} \sum_{j} \sum_{k} \phi(x_{i},y_{j},z_{k})$$
where $x_{i}, y_{j}, z_{k}$ are the midpoints of the 1D functions (I simply applied the composite midpoint rule using $g(y,z) = \int_{a}^{b} \phi(x,y,z) dx$. I did it similarly with the trapezoid rule.
My issue is that whatever subdivision I take to apply my composite rules, the integrals seem to converge to a value that underestimates the analytical value ($2.6 \times 10^{6}$ instead of $2.8 \times 10^{8}$).
What could be the origin of such a convergence issue ? How could the Matlab integral3 algorithm converge to the analytical value and not the different rules I used ?
I know my question is a little vague but the function I am integrating cannot be written with the usual analytical functions and I don't really know its behavior. I just know that it its gets "really small, really fast" with given $y,z$ as $x$ grows. I was hoping I could get some insight on what could be causing this convergence error.
Thank you in advance for any help or remark you could provide