Relativistiv kinetic energy and floating point

Question

My function is $E(v)=mc^2(\frac{1}{\sqrt{1-v^2/c^2}} - 1)$, (c=3e8, m=1) and I have to calculate it for values of v between 1e-6 and 2.99e8. The point of this problem is floating point precision.

For very small values of v the result is E=0. I figured that's because of the root: $\sqrt{1-v^2/c^2}=1$.

I tried using Taylor: $mv^2/2 + 3mv^4/8c^2 + 5mv^6/16c^4 + 35mv^8/128c^6$, and that works well for small v but gives me wrong results for very big values of v.

How do I transform E(v) to a better formula that works for all values of v?

Answer 1

I'm assuming you don't have a square root routine available. Something you could try...

For simplicity I'll assume $m=c=1$, and call the answer $u$, so we have

$$E(v)=u=\frac{1}{\sqrt{1-v^2}}-1$$

We can rearrange this to get:

$$(u+1)^2(1-v^2)=1$$

And we can write it as a root finding problem in $u$:

$$(1-v^2) u^2 + 2(1-v^2)u + (v^2 - 2v) = 0$$

Or, perhaps it's clearer to let $s=1-v^2$ and rearrange,

$$s(u^2 + 2u + 1)-1 = 0$$

You could use any root finding algorithm, to find $f(u)=0$, Newton's method iterates:

$$u_{n+1} = u_n - \frac{f(u_n)}{f'(u_n)} = u_n - \frac{s (u_n^2 + 2u_n + 1)-1}{2s(u_n + 1)}$$

With a few iterations that should give you the answer to a decent accuracy, at least.

Answer 2

You can not have a formula which is numerically reliable for all values of $v$.

The fundamental problem is that relativistic $\gamma$ factor, i.e., $$ \gamma(x) = \frac{1}{\sqrt{1 - x^2}}, \qquad x = \frac{v}{c}$$ is ill-conditioned in the vicinity of $x = 1$. Specifically, the condition number of $\gamma$ satisfies $$ \kappa_{\gamma}(x) = \left|\frac{x \gamma'(x)}{\gamma(x)}\right| = \frac{x^2}{1-x^2} \rightarrow \infty, \quad x \rightarrow 1, \quad x \in [0,1)$$ If $\hat{x}$ is the computed value of $x$, then $$ \left|\frac{\gamma(x) - \gamma(\hat{x})}{\gamma(x)}\right| \approx \kappa_\gamma(x) \left|\frac{x - \hat{x}}{x}\right| $$ Any rounding errors committed when computing $x$ will be magnified by $\kappa_\gamma(x)$ when computing $\gamma(x)$. The effect is dramatic if $v \approx c$.

However, you can get pretty close to $c$ before this becomes an issue. For $v = 0.999c$ the condition number is still less than 500!

The expression $\gamma(x) - 1$ suffers from subtractive cancellation when $x \approx 0$. The cure is to do a Taylor expansion. The order hinges on the accuracy which you seek.

When is the expression $\gamma(x) - 1$ safe to evaluate? The subtraction $d = a - b$ can be computed safely, when $|a| \ge 2 |b|$ or $2|a| \leq |b|$. In our case the relevant point is $x_0 =\sqrt{\frac{3}{4}} \approx 0.8660$. You can use a Taylor expansion on the interval $[0,x_0]$ and the original expression on the interval $(x_0,1)$.

It is frequently necessary to use multiple expression to cover the domain of a given function. Having a single expression is a rare luxury.