Solving (Frenet-Serret) differential equation system in Matlab

Question

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!

Answer 1

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.

Answer 2

Use ODE15s: http://www.mathworks.com/help/matlab/ref/ode15s.html

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.