Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
add comment

2 Answers 2

up vote 5 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
add comment

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}
   ]
  ]

Manipulate animation


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

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.