# Finding a 3d curve from torsion and curvature with NDSolve

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.

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[], y[ti] == iniPos[], z[ti] == iniPos[],
x'[ti] == iniDir[], y'[ti] == iniDir[], z'[ti] == iniDir[],
nx[ti] == iniNor[], ny[ti] == iniNor[], nz[ti] == iniNor[],
bx[ti] == iniBin[], by[ti] == iniBin[], bz[ti] == iniBin[]},
{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...

• 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!!
– 0xFE
Dec 20 '13 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 == t0,
n == n0,
b == b0,
r == 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 == 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},             (* initial frame *)
{r, t, n, b},                    (* variables for curve, tangent, normal, binormal *)
{s, 0, 12}]                      (* interval of integration *)

• I think the sentence "NDSolve can solve vector equations" ought to be in big, bold font. Dec 20 '13 at 12:35
• 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. May 23 '15 at 7:15