# 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[[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...

-
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], where r is the curve. 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.

Clear[r, s, t, n, b, κ, τ, r0, frame0];
(*fsmat={{0, κ[s], 0}, {-κ[s], 0, τ[s]}, {0, -τ[s], 0}};*)   (* explicit matrix *)
fsmat = {κ[s], 0, -τ[s]} \[TensorWedge] {0, 1, 0} // Normal; (* computed matrix *)
frame = {t[s], n[s], b[s]};
eqns = {
D[frame, s] == fsmat.frame,
r'[s] == t[s],
frame == frame0 /. s -> 0,
r[0] == r0};

κ[s_] := 10.08 Sin[Pi s]^2;
τ[s_] := 3 Cos[Pi s];
frame0 = Orthogonalize[RandomReal[1, {3, 3}]]; (* random orientation *)
r0 = {0, 0, 0};

sol = First@NDSolve[eqns, {r, t, n, b}, {s, 0, 12}];

-
I think the sentence "NDSolve can solve vector equations" ought to be in big, bold font. – belisarius has settled 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. – J. M. May 23 at 7:15