I'm trying to use the Frenet–Serret formulas to find the curve that matches the torsion and curvature I specify numerically with an InterpolatingFunction.

The problem is that the system is overdetermined, and have been getting all kinds of errors from Mathematica. If I don't tell NDSolve about the definition of N as the normal vector and B as the cross product of T and N, then I have a 9x9 system, which will solve on occasion, but with results don't seem right.

The curve should be a Mobius strip with length $2\pi$, so I also have the boundary conditions of $T(0)=T(2\pi), N(0)=N(2\pi), B(0)=B(2\pi)$, but when I include those conditions, I get

"There are significant errors {0., 0., 0., 0., 0., 0., 2.49764*10^-6, 0., 0., 
-2.06667*10^-6, 0., 0.} in the boundary value residuals. Returning the best solution 

Some errors are to be expected, and $10^{-6}$ is not significant to me, but the resulting curve is a straight line from $(0,0,0)$ to $(\pi/2,0,0)$.

When I try specifying $N(s)$ as $\frac{T'(s)}{\|T'(s)\|}$ and $B(s)$ as $T\times N$, and only using the first 2 equations on the Wikipedia page, I get:

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. 
NDSolve will try solving the system as differential-algebraic equations.`

NDSolve::ndcf: Repeated convergence test failure at s == 0.; unable to continue.

As far as I can tell, NDSolve doesn't support vector equations, so I have been breaking everything down into the component equations.

I'd be happy to share the points that form the interpolating functions for torsion and curvature, but I don't know where to post them: it's almost 1500 points.


2 Answers 2

x''[t] == k1[t] nx[t], nx'[t] == -k1[t] x'[t] - k2[t] bx[t], bx'[t] == k2[t] nx[t],
y''[t] == k1[t] ny[t], ny'[t] == -k1[t] y'[t] - k2[t] by[t], by'[t] == k2[t] ny[t],
z''[t] == k1[t] nz[t], nz'[t] == -k1[t] z'[t] - k2[t] bz[t], bz'[t] == k2[t] nz[t],
x[ti] == iniPos[[1]], y[ti] == iniPos[[2]], z[ti] == iniPos[[3]],
x'[ti] == iniDir[[1]], y'[ti] == iniDir[[2]], z'[ti] == iniDir[[3]],
nx[ti] == iniNor[[1]], ny[ti] == iniNor[[2]], nz[ti] == iniNor[[3]],
bx[ti] == iniBin[[1]], by[ti] == iniBin[[2]], bz[ti] == iniBin[[3]]},
{x, y, z, nx, ny, nz, bx, by, bz}, {t, ti, te}]

k1[t] is curvature, k2[t] is torsion, {x[t],y[t],z[t]} is path, {nx[t],ny[t],nz[t]} is normal, {bx[t],by[t],bz[t]} is binormal, arc length is (te-ti)*norm[iniDir].

This is signed curvature...

  • $\begingroup$ You solution gave me a curve I had actually seen before, but discarded because it wasn't closed. I went back and saw that I was actually using unsigned curvature. I fixed that and got a great result. Thanks!! $\endgroup$
    – 0xFE
    Commented Dec 20, 2013 at 7:43

NDSolve can solve vector equations. To the Frenet-Serret equations, in which t, b, n are the tangent, normal, and binormal resp., I added r'[s] == t[s], which will cause the parametrization of the curve r[s] to be integrated along with the frame. NDSolve will recognize the system as a system of vector equations if the initial values are vectors.

Clear[r, s, t, n, b, κ, τ, r0, t0, n0, b0];
eqns = {
   t'[s] == κ[s] n[s],
   n'[s] == -κ[s] t[s] + τ[s] b[s],
   b'[s] == -τ[s] n[s],
   r'[s] == t[s],
   t[0] == t0,
   n[0] == n0,
   b[0] == b0,
   r[0] == r0};

κ[s_] := 10.08 Sin[Pi s]^2;
τ[s_] := 3 Cos[Pi s];
{t0, n0, b0} = Orthogonalize[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}];
r0 = {0, 0, 0};
sol = First @ NDSolve[eqns, {r, t, n, b}, {s, 0, 12}];

With[{s1 = (r /. sol)["Domain"][[1, 1]], 
   s2 = (r /. sol)["Domain"][[1, 2]]},
    ParametricPlot3D[Evaluate[r[s] /. sol], {s, s1, s2},
     PlotStyle -> {Thick, Brown}, PlotRangePadding -> 1],
     Dynamic@Translate[{Thick, Red, Arrow[{{0, 0, 0}, t[s0]}], Green, 
         Arrow[{{0, 0, 0}, n[s0]}], Blue, 
         Arrow[{{0, 0, 0}, b[s0]}]} /. sol, r[s0] /. sol]
   {s0, s1, s2}

Manipulate animation

One can also enter the equations in matrix form and NDSolve will sort them out. (Update: Packaged as a function.)

fsIntegrate[κ_, τ_, {r0_, frame0_}, {r_, t_, n_, b_}, {s_, s0_, s1_}] :=
    fsmat = {κ, 0, -τ} \[TensorWedge] {0, 1, 0} // Normal,  (* Frenet-Serret RHS *)
    frame = {t[s], n[s], b[s]}},                            (* TNB frame *)
     eqns = {D[frame, s] == fsmat.frame, r'[s] == t[s],     (* ODE *)
       frame == frame0 /. s -> 0, r[0] == r0}},             (* ICs *)
    NDSolve[eqns, {r, t, n, b}, {s, s0, s1}]

 10.08 Sin[Pi s]^2, 3 Cos[Pi s],  (* curvature, torsion *)
 {{0, 0, 0},                      (* initial point *)
  IdentityMatrix[3]},             (* initial frame *)
 {r, t, n, b},                    (* variables for curve, tangent, normal, binormal *)
 {s, 0, 12}]                      (* interval of integration *)
  • 2
    $\begingroup$ I think the sentence "NDSolve can solve vector equations" ought to be in big, bold font. $\endgroup$ Commented Dec 20, 2013 at 12:35
  • 2
    $\begingroup$ Back when version 5 was new, they actually promoted NDSolve[]'s brand-new ability to handle vector-valued DEs by using the Frenet-Serret equations as an example. $\endgroup$ Commented May 23, 2015 at 7:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.