Finding a 3d curve from torsion and curvature with NDSolve

Question

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 
found."

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.

Answer 1

First@NDSolve[{
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...

Answer 2

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]]},
  Manipulate[
   Show[
    ParametricPlot3D[Evaluate[r[s] /. sol], {s, s1, s2},
     PlotStyle -> {Thick, Brown}, PlotRangePadding -> 1],
    Graphics3D[
     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}
   ]
  ]

---

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

ClearAll[fsIntegrate];
fsIntegrate[κ_, τ_, {r0_, frame0_}, {r_, t_, n_, b_}, {s_, s0_, s1_}] :=
  With[{
    fsmat = {κ, 0, -τ} \[TensorWedge] {0, 1, 0} // Normal,  (* Frenet-Serret RHS *)
    frame = {t[s], n[s], b[s]}},                            (* TNB frame *)
   With[{
     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}]
    ]];

fsIntegrate[
 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 *)