Solving (Frenet-Serret) differential equation system in Matlab

I'm about (trying) to solve the Frenet-Serret equation given by the known formulas, finding $e(s)$, $n(s)$, $b(s)$, where

$e'(s) = \kappa(s)v(s)n(s)$

$n'(s) = -\kappa(s)v(s)e(s) + \tau(s)v(s)b(s)$

$b'(s) = -\tau(s)v(s)n(s)$

with the initial values $e(0)$, $n(0)$, $b(0)$, and given $\kappa$, $\tau$ and $v$.

I have to solve this in Matlab, but I've got no idea, what function to use and how to parametrize it. Tried dsolve, but got syntax errors all the way.

Help appreciated!

-
I haven't worked with Frenet-Serret equations before, but from what I can tell, $\kappa$, $\tau$ and $\nu$ are constants, right? Then why do you write them as if they were function of s? – If you have problems with the implementation, you should post your code here and describe what exactly is the problem. – A. Donda Dec 6 '13 at 1:14
Nope, $\kappa$, $\tau$ and $v$ are might not constant in every scenario. Curvature and torsion, as well as velocity can change in different points. For example: Viviani's curve I also totally lack a code, because I've got no clue which function should I use and mostly, how. I'm very novice in Matlab. – Dyin Dec 6 '13 at 8:29
Well, if the parameters are a function of $s$, then it's not a linear differential equation system anymore – which changes the approach, lde have closed-form solutions while nonlinear ones do not necessarily, so you need numerical integration. – But you must have tried something, since you write about syntax errors. I'm telling you, you'll have a much better chance at getting help if you show people in detail what you've done and where exactly you got stuck. – A. Donda Dec 6 '13 at 13:44
Also, consider moving the question to stackoverflow, most Matlab programming experts are active there. – A. Donda Dec 6 '13 at 13:45
Are you assuming that the curvature is non-vanishing? There are a couple of typos in your question. For non-unit speed curves, the first equation's right-hand side is missing a factor of $v(s)$. Also, the third equation should have left-hand side $b'(s)$, not $n'(s)$ – AppliedSide Dec 6 '13 at 20:23

Unless your system's parameters satisfy specific conditions (e.g., zero curvature, zero torsion, or this more complicated case), it won't have a closed form solution and dsolve will not be helpful. You'll need to numerically integrate the differential equation along the arc length $s$. The equations are likely non-stiff so the standard ode45 solver should suffice. Here's some basic code with constant parameter functions, but it also shows how you might convert those into functions of arc length. Also, your state variables are vectors so this is a 9-D system for the basic case:

a = 2;
b = 2;
T = 2*pi*sqrt(a^2+b^2);
kappa = abs(a)/(a^2+b^2);
nu = 1;
tau = b/(a^2+b^2);
kappa = @(s)kappa;
nu = @(s)nu;
tau = @(s)tau;
A = @(s)[              0 kappa(s)*nu(s)            0;
-kappa(s)*nu(s)              0 tau(s)*nu(s);
0  -tau(s)*nu(s)            0];
f = @(s,y)reshape(A(s)*reshape(y,[3 3]),9,1); % Handle matrix ODE
sspan = linspace(0,T,20);
y0 = eye(3); % [e1 e2 e3;n1 n2 n3;b1 b2 b3]
[s,y] = ode45(f,sspan,y0(:));
y = reshape(y.',[3 3 numel(y)/9]);

% Three basis vectors as functions of arc length
E = y(:,[1 4 7]).'; % [e11 ... e1m;e21 ... e2m;e31 ... e3m]
N = y(:,[2 5 8]).'; % [n11 ... n1m;n21 ... n2m;n31 ... n3m]
B = y(:,[3 6 9]).'; % [b11 ... b1m;b21 ... b2m;e31 ... b3m]

% Circular helix
t = s.'/sqrt(a^2+b^2);
x = a*cos(t);
y = a*sin(t);
z = b*t;

figure
plot3(x,y,z,'r')
hold on
axis equal
grid on

Note that the initial conditions, y0 ($e(0)$, $n(0)$, $b(0)$), must be non-zero vectors as the origin is a fixed point of this system.

-
The curvature for a spiral is $r/(a^2+r^2)$ and torsion is $a/(a^2+r^2)$, where $r$ is the radius and $a$ would define the steep in the 3rd dimension. Both curvature and torsion are constants. – Dyin Dec 14 '13 at 18:33
The Frenet–Serret frame should not be confused with the Cartesian or inertial frame ($R^3$). All of those pretty pictures in the Wikipedia article and elsewhere, show curves in the standard Cartesian coordinates for the most part because that's what our brain understands. Some show how the Frenet–Serret frame moves along these curves as a little coordinate system, but they don't generally plot out the Frenet coordinates directly like my plots. – horchler Dec 14 '13 at 18:52
Thank you for your kind answer! How could I make your solution to create such a pretty plot of the original curve reconstructed from curvature and torsion? Could you show me an example on a spiral? – Dyin Dec 14 '13 at 20:44
@Dyin: I just updated my answer. I looked at your equations and then investigated the Frenet–Serret equations further. The state variables are not scalars as your notation indicates, but vectors, which changes things a lot. The code outputs the three basis vectors. I also plotted a helix for you. Assuming all is correct, you'll need to figure out for yourself how to plot little coordinate frames along the helix if that's what you want. Your question said nothing about helices, or plotting or converting back from the Frenet–Serret so I think all of those are best left for new questions. – horchler Dec 14 '13 at 22:58
First of all, thank you for your kind update! I see the TNB frame in y, each vector is orthogonal, so it must be fine. I bet if I want to plot the original space curve, I need to integrate $e(s)$ to receive the parametrization by arc length. Should I create a new question for that? – Dyin Dec 16 '13 at 10:03

Look at Example 1 and the otyers . Change e to y_1, n to y_2 and b to y_3, and plug in your known values for the torsion, curvature, and for v.

-