Approximate $|x|$ with a smooth function

A bit late to the party, but you can smoothly approximate $f(x) = |x|$ by observing $\partial f/\partial x = \mbox{sgn}(x)$ for $x \neq 0$. Therefore approximating the $\mbox{sgn}(x)$ function by

$$ f(x) = 2\mbox{ sigmoid}(kx)-1$$

(the $k$ being a parameter that allows you to control the smoothness), we get

$$ \partial f/\partial x = 2\left (\frac{e^{kx}}{1+e^{kx}} \right ) - 1 $$ $$ \Rightarrow f(x) = \frac{2}{k}\log(1+e^{kx})-x-\frac{2}{k}\log(2)$$

where the constant term was chosen to ensure $f(0) = 0 $.

Included are plots of the function for $x \in [-5,5]$, where $k=1$ is red, $k=10$ is blue and $k=100$ is black.

Note that if you wanted to use the more smooth $k=1$ in the interval $[-5,5]$ it may be worth applying a further linear transformation i.e. ~$5f(x)/3.6$ to ensure the values at the edge of the interval are correct.

I have used your $\sqrt{x^2+\epsilon}$ function once before, when an application I was working on called for a curve with a tight radius of curvature near $x=0$. Whether it is best for your purpose might depend on your purpose.

(Side note: The derivative of $|x|$ does not exist at $x=0$. In physics, we sometimes cheat and write ${d\over dx}|x| = \rm{sgn}(x)$, the "sign" function that gives $-1$ when $x<0$, $0$ when $x=0$, and $+1$ when $x > 0$. But only when we know that the value at $x=0$ will not get us into trouble!)

Another solution would be the function $x \rightarrow x \frac{e^{kx}}{e^{kx} + e^{-kx}} -x\frac{e^{-kx}}{e^{kx} + e^{-kx}} $ with higher the $k$, closer you will be to the absolute value function. Furthermore it is equal to 0 in 0 wich is sometimes a desirable feature (compared to the solution with the root square stated above). However you lose the convexity.