# Computing numerical derivatives

I am trying to create a sweeping surface, for which I need the frenet frame of a curve. I am trying to compute this for arbitrary curves but for testing I am just using the parametric unit half circle. i.e. $$C(t) = (\cos(t), \sin(t), 0)$$

I am computing my derivatives using central differences, i.e.:

pub fn frenet_trihedron(x: f32, f: &dyn Fn(f32) -> Vec3) -> (Vec3, Vec3, Vec3)
{
let h = 1e-6 as f32;

let normal = (-2.0 * f(x) + f(x + h) + f(x - h)) / (h * h);
let tangent = (f(x + h) - f(x - h)) / (2.0 * h);

let tangent = tangent.normalize();
let normal = normal.normalize();
let binormal = tangent.cross(&normal);

println!("t {}, n {}, b {}", tangent, normal, binormal);

(tangent, normal, binormal)
}


This however seems to yield very impossible results, for example:

t
┌            ┐
│ -0.8274758 │
│ 0.56150144 │
│          0 │
└            ┘

, n
┌   ┐
│ 0 │
│ 1 │
│ 0 │
└   ┘

, b
┌            ┐
│          0 │
│          0 │
│ -0.8274758 │
└            ┘


That normal is not orthogonal to the tangent not even close What am I doing wrong?

This is the calling code:

let f = |v| Vec3::new(f32::cos(v * PI), f32::sin(v * PI), 0.0);
let (t, n, b) = frenet_trihedron(v, &|x| f(x));


It seems the error is being cause by floating point imprecision. Using 1e-3 asmy epsilon yields remarkably better results. so maybe f32 is just not precise enough.

• Show us the implementation of the function you are trying to differentiate, as well as the driver code that shows how you call the frenet_trihedron function. Commented Feb 16, 2023 at 0:43
• I don't think that's the source of your bug, but if you are doing numerical differentiation good practice is using $h \approx eps^{1/2} x$, to minimize the numerical error. Your $h$ is much smaller than that, so you will incur in a larger error from the cancellation in the numerator. Commented Feb 16, 2023 at 8:45
• Could you add the value of v that produced your example output? This will tell at least tell us (and you) if the bug is the wrong calculation of t or n (or both). This is an important debugging principle: narrow down the position of the error. Commented Feb 16, 2023 at 8:46
• @Federico Poloni : For the symmetric differences the error is quadratic, so the optimal $h$ is $\epsilon^{1/3}\max(1,|x|)$ with relative error of the derivative of about $\epsilon^{2/3}$. So $h\approx 10^{-3}$ as you gave, and up to $h\approx 2^{-8}=0.0039$, times the size of $x$. Commented Feb 16, 2023 at 10:08
• @Federico In the numerical evaluation you get $$\Big|\frac{f(x+h)-f(x-h)}{2h}-f'(x)\Big|\le \frac{2(|f|(x)+|f'(x)x|)\mu}{2h}+\frac{h^2}{6}|f'''(x)|+O(h^4),$$ and the minimum of the two first terms on the right is at $h=(C\mu)^{1/3}$, where $C$ captures the function-dependent factors. For large $x$ this constant becomes dominantly linear in $|x|$. So the actual estimate for the best $h$ should look a bit different again. Its not an exact science. A bit of error can be mitigated by computing $\frac{f(x+h)-f(x-h)}{(x+h)-(x-h)}$, so that the denominator is the actual difference of the arguments. Commented Feb 16, 2023 at 11:15

The unit tangent vector is $$\hat t(x)=\frac{C'(x)}{\|C'(x)\|}$$. The normal direction is the derivative of this unit tangent, not just the second derivative $$n(x)=\hat t'(x)=\frac{(\|C'(x)\|^2I-C'(x)C'(x)^T)C''(x)}{\|C'(x)\|^3} =\frac{(I-\hat t(x)\hat t(x)^T)C''(x)}{\|C'(x)\|}$$ then again the unit normal is the normalization after this projection $$\hat n(x)=\frac{n(x)}{\|n(x)\|}$$. The scalar denominator cancels in the normalization, so can be left out from the start

The winding direction after that is then indeed the cross product.

In Python this can look like

from numpy import array, cross, cos, sin, pi
from numpy.linalg import norm
def C(x): return np.array((cos(pi*x),sin(pi*x),0*x))

def dreibein(C,x):
h = 0.005*max(1,abs(x))
Cm,C0,Cp = C(x-h),C(x),C(x+h)
#tangent
t = Cp-Cm; t/=norm(t)
# normal
n = Cp-2*C0+Cm; n -= n.dot(t)*t; n /= norm(n)
# winding
b = cross(t,n)
return t,n,b

t,n,b =dreibein(C,0.31019)
print(t,n,b)


returns

[-0.82741594  0.56158959  0.        ] [-0.56158959 -0.82741594  0.        ] [0. 0. 1.]


which visibly is an orthonormal basis.

The same basis can be found from the columns of $$Q$$ in a suitable QR decomposition

Q,R = np.linalg.qr(np.transpose([Cp-Cm,Cp-2*C0+Cm]),"complete")
Q = Q @ np.diag([*np.sign(np.diag(R)),1])
t,n,b = Q.T