find $ \int_0^4\int_0^4\int_0^4 \sqrt{x^2+y^2+z^2}\,dx\,dy\,dz$

Numerical integration gives that the closest integer is $246$, since $$I=\int_{[0,4]^3}\sqrt{x^2+y^2+z^2}\,d\mu=245.91154\ldots$$ The best way to achieve a good accuracy "by hand" is to study the pdf of $\sqrt{X^2+Y^2+Z^2}$ when $X,Y,Z$ are three independent random variable with a uniform distribution over $[0,4]$, by eventually using some quantitative form of the Central Limit Theorem.

Since $\mathbb{E}[X^2]=\frac{16}{3}$, the Jensen's inequality gives $I\leq 4^4=256$.

Anyway, such integral also admits the closed form: $$ I = \frac{32}{3} \left(6 \sqrt{3}-\pi +\log\left(3650401+2107560 \sqrt{3}\right)\right).$$

$\newcommand{\+}{^{\dagger}} \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\down}{\downarrow} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\isdiv}{\,\left.\right\vert\,} \newcommand{\ket}[1]{\left\vert #1\right\rangle} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert} \newcommand{\wt}[1]{\widetilde{#1}}$ With $\ds{\pars{~\mbox{spherical coordinates}~}\ r = \root{x^{2} + y^{2} + z^{2}}}$, we have $\ds{\nabla\cdot\pars{r\,\vec{r}} = {\vec{r} \over r}\cdot\vec{r} + 3r = 4r}$:

---

\begin{align}&\int_{0}^{4}\int_{0}^{4}\int_{0}^{4}% \root{x^{2} + y^{2} + z^{2}}\,\dd x\,\dd y\,\dd z =\int_{\bracks{0,4}^{3}}r\,\dd^{3}\vec{r} =\int_{\bracks{0,4}^{3}}{1 \over 4}\,\nabla\cdot\pars{r\,\vec{r}}\,\dd^{3}\vec{r} \\[3mm]&={1 \over 4}\int_{\rm S}r\,\vec{r}\cdot\dd\vec{\rm S} ={1 \over 4}\int_{\braces{0,4}^{2}}\root{16 + y^{2} + z^{2}}\,4\,\dd y\,\dd z \\[3mm]&+{1 \over 4}\int_{\braces{0,4}^{2}}\root{x^{2} + 16 + z^{2}}\,4 \,\dd x\,\dd z +{1 \over 4}\int_{\braces{0,4}^{2}}\root{x^{2} + y^{2} + 16}\,4\,\dd x\,\dd y \end{align}

---

\begin{align}&\int_{0}^{4}\int_{0}^{4}\int_{0}^{4}% \root{x^{2} + y^{2} + z^{2}}\,\dd x\,\dd y\,\dd z =3\int_{\braces{0,4}^{2}}\root{x^{2} + y^{2} + 16}\,\dd x\,\dd y \\[3mm]&=3\int_{0}^{4}\int_{0}^{4}\root{x^{2} + y^{2} + 16}\,\dd y\,\dd x =192\int_{0}^{1}\int_{0}^{1}\root{x^{2} + y^{2} + 1}\,\dd y\,\dd x \\[3mm]&=48\int_{0}^{1}\left\{2\left[\root{x^{2} + 2} + \left(x^{2} + 1\right) \ln\left(\root{x^{2} + 2} + 1\right)\right] - \left(x^{2} + 1\right) \ln\left(x^{2} + 1\right)\right\}\,\dd x \\[3mm]&=\bbox[15px,border:1px solid navy]{\ds{ 64\left[\root{3} + \ln\left(7 + 4\root{3}\right)\right] - {32\pi \over3}}} \approx {\tt 245.9115409} \end{align}

// monteCarlo0.cc
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
using namespace std;
typedef long double ldouble;
typedef unsigned long long ullong;
const ldouble RANDMAX1=ldouble(RAND_MAX) + 1.0L;
const ullong ITER=100000000ULL; // One hundred million
ldouble X,Y;

inline ldouble f(ldouble x,ldouble y)
{
 return sqrt(x*x + y*y + 1.0L);
}

inline ldouble f()
{
 X=rand()/RANDMAX1;
 Y=rand()/RANDMAX1;
 return (f(X,Y) + f(X,1.0L - Y) + f(1.0L - X,Y) + f(1.0L - X,1.0L - Y))/4.0L;
}

int main()
{
 ldouble suma=0;

 for ( ullong n=0 ; n<ITER ; ++n ) suma+=f();
 cout<<setprecision(50)<<(192.0L*(suma/ITER))<<endl;

 return 0;
}

It runs for 7 sec.

Result: $\ds{\color{#f00}{ 245.911}07650843407733676215798368502873927354812622}$.

And now for something entirely random:

Monte Carlo integration gives

Monte Carlo (n=1000) estimate=245.876995211
Monte Carlo (n=10000) estimate=246.605709455
Monte Carlo (n=100000) estimate=245.440557832
Monte Carlo (n=1000000) estimate=245.841098533
Monte Carlo (n=10000000) estimate=245.946260986

Generated from the following:

import math
import random

def norm(x):
    return math.sqrt(sum([xi*xi for xi in x]))

# returns a uniformly distributed sample in [0, 4]^3
def unif():
    return [random.uniform(0, 4) for i in (1, 2, 3)]

area = 4*4*4
for total_count in (1000, 10000, 100000, 1000000, 10000000):
    sum_func = 0
    for i in xrange(total_count):
        x_samp = unif()
        sum_func += norm(x_samp)

    est_volume = sum_func/float(total_count)*area

    print("Monte Carlo (n=%s) estimate=%s" % (total_count, est_volume))