Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
up vote 9 down vote accepted
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...

share|improve this answer
    
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[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}
   ]
  ]

Manipulate animation


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 *)
share|improve this answer
2  
I think the sentence "NDSolve can solve vector equations" ought to be in big, bold font. – Dr. belisarius Dec 20 '13 at 12:35
1  
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 '15 at 7:15

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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