……不是沒有道理的。——常規的高效的迭代方法不行的全域性優化，本質上是一撥又一撥的蒙特卡洛。

這個例子可以讓常規的數值積分方法統統歇菜。比如， 如果用matlab中的integral3，肯定不行。需要非常特殊的技巧。Mathematica預設的方法一向通用性強，碰上這個問題也十分吃緊。

下面這兩個結果，蒙特卡洛方法的結果更接近準確值。

精度要求低或者精度無法保證高而計算量大的情況下，蒙特卡洛可以試試。隨機全域性優化大致如此。

求不規則形狀體積：

隱函式曲面方程

這個例子的較為精確的數值結果是：67.59353836486223222964871：

Matlab integral3歇菜

這個函式也是比較新的版本才有的。但是功能就比較老舊了。
無法得到有用結果，返回NaN和警告資訊。

Mathematica NIntegrate

用Mathematica有幾種不同實施方法，傳統的方法差太遠。

比如：

NIntegrate[Boole[Cos[x]+Cos[y]+Cos[z]+0.51` (Cos[x] Cos[y]+Cos[y] Cos[z]+Cos[z] Cos[x])+1<=0],{x,0,2Pi},{y,0,2Pi},{z,0,2Pi}]

得到：
66.70206691379813`
差距比較明顯；

NIntegrate
 
[1,{x,y,z}\[Element]ImplicitRegion[1+Cos[x]+Cos[y]+Cos[z]+51/100 (Cos[x] Cos[y]+Cos[x] Cos[z]+Cos[y] Cos[z])<=0&& 0<=x<=2Pi &&0<=y<=2 Pi &&0<=z<=2Pi,{x,y,z}]]

最新版本
67.0421386455621
這是稍好的近似值。

Monte Carlo方法，效率低，精度更差而且不太可控

ClearAll["Global`*"]; n = 125000; SeedRandom
 
[0];
f[x_] := (Cos /@ x // Total) + 
  51/100*Dot[Cos /@ RotateRight@x, Cos /@ x] + 1
Count[Flatten@(f /@ RandomReal[{0, 2 Pi}, {n, 3}]), u_ /; u <= 0]/
  n*(2.0 Pi)^3

結果
67.34464074875744`