Tengo que escribir un programa que obtenga dos funciones (curvatura y torsión) y 3 vectores del "marco" de Frenet-Serret en el punto de partida - y tengo que reconstruir la curva espacial a partir de estos datos de entrada dados.
Así que he escrito algo de código Octave, que funciona bien, por ejemplo en la curva
Phi(t) = (cos(t/sqrt(2)), sin(t/sqrt(2)), t/sqrt(2)) % natural parametrization(es una especie de espiral) y me alegré porque al trazar los resultados, la curva original era casi la misma que la reconstruida. Pero cuando probé alguna otra parametrización, había una gran diferencia entre la reconstruida y la original:
Phi(t) = (cos(t), sin(t), t) % other parametrization of the same curveRealmente no puedo encontrar lo que hice mal en mi código... Tal vez mi idea era completamente errónea, y sólo puedo reconstruir las curvas que es naturalmente parametrizado?
% note: this returns a vector
% K : curvature
% T : torsion
function equation = f(s, gvect, K, T, velocity)
gs = toStruct(gvect);
e_v = velocity(s)*K(s).* gs.n;
n_v = -velocity(s)*K(s).* gs.e + velocity(s)*T(s).* gs.b;
b_v = -velocity(s)*T(s).* gs.n;
equation = [e_v; n_v; b_v];
return;
endfunction
% returns a struct
function str = toStruct(vect)
str.e = vect(1:3);
str.n = vect(4:6);
str.b = vect(7:9);
endfunction
% returns a vector
function vect = toVect(str)
vect = [str.e; str.n; str.b];
endfunction
% numerical integration
function p = nintegrate(x, fx)
acc = 0;
for i = [1:length(x)-1]
p(i) = acc;
delta_x = x(i+1) - x(i);
acc += delta_x * fx(i);
endfor
endfunction
% K : curvature
% T : torsion
function gvals = RungeKutta4(S, g, K, T, velocity, delta)
for i = [1:length(S)-1]
g_i = toVect(g(i));
k1 = f(S(i), g_i, K, T, velocity);
k2 = f(S(i)+delta/2, g_i+delta/2*k1, K, T, velocity);
k3 = f(S(i)+delta/2, g_i+delta/2*k2, K, T, velocity);
k4 = f(S(i)+delta, g_i+delta*k3, K, T, velocity);
g(i+1) = toStruct( g_i + delta/6*(k1+2*k2+2*k3+k4) );
endfor
gvals = g;
endfunction
function drawPlots(xs, ys, zs, orig_xs, orig_ys, orig_zs)
plot3(xs, ys, zs, ";reconstructed;", orig_xs, orig_ys, orig_zs, ";original;");
endfunction
function Phi = reconstruct(K, T, e0, n0, b0, alpha, beta, start, delta = 0.1, velocity = @()1)
S = [alpha:delta:beta];
g(1).e = e0; % unit tangent vector
g(1).n = n0; % principal normal vector
g(1).b = b0; % binormal vector
g = RungeKutta4(S, g, K, T, velocity, delta);
e_x = cellfun(@(v)v(1), {g.e});
e_y = cellfun(@(v)v(2), {g.e});
e_z = cellfun(@(v)v(3), {g.e});
Phi.x = nintegrate(S, e_x) + start(1);
Phi.y = nintegrate(S, e_y) + start(2);
Phi.z = nintegrate(S, e_z) + start(3);
return;
endfunction
% works OK for this
function spiralNaturalParametrized()
L = sqrt(2)*4*pi;
e0 = [ 0, 1/sqrt(2), 1/sqrt(2)]'; % unit tangent vector
n0 = [-1, 0, 0]'; % principal normal vector
b0 = [ 0, -1/sqrt(2), 1/sqrt(2)]'; % binormal vector
Phi = reconstruct(@()1/2, @()1/2, e0, n0, b0, 0, L, [1 0 0], 0.05);
t = [0:0.05:L];
drawPlots(Phi.x, Phi.y, Phi.z,
cos(t./sqrt(2)), sin(t./sqrt(2)), t./sqrt(2));
endfunction
% another parametrization: huge error
function spiralAnother()
L = 4*pi;
e0 = [ 0, 1/sqrt(2), 1/sqrt(2)]';
n0 = [-1, 0, 0]';
b0 = [ 0, -1/sqrt(2), 1/sqrt(2)]';
Phi = reconstruct(@()1/2, @()1/2, e0, n0, b0, 0, L, [1 0 0], 0.05, @(v)sqrt(2));
t = [0:0.05:L];
drawPlots(Phi.x, Phi.y, Phi.z, cos(t), sin(t), t);
endfunctionHe incluido mi código. Se agradece cualquier ayuda.
Actualización:
He creado otra prueba que revela el hecho de que el algoritmo tiene algún error en su interior incluso cuando se trata del caso paramétrico natural:
function YetAnotherTest()
L = 10*pi;
e0 = [ 0, 2/sqrt(13), 3/sqrt(13)]';
n0 = [ -1, 0, 0]';
b0 = [ 0, -3/sqrt(13), 2/sqrt(13)]';
Phi = reconstruct(@()(2/13), % curvature
@()sqrt(3/13), % torsion
e0, n0, b0, % initial frame
0, L, % interval
[2 0 0], 0.05); % start point and step size
t = [0:0.05:L];
drawPlots(Phi.x, Phi.y, Phi.z,
2*cos(t./sqrt(13)), 2*sin(t./sqrt(13)), t*(3/sqrt(13)));
endfunctionNota: esta es una prueba para $$Phi(t) = ( 2*\cos(t/sqrt(13)), 2*\sin(t/sqrt(13)), 3t/sqrt(13) )$$
que es una curva paramétrica natural. Si echas un vistazo en esta parcela se puede ver que la curva reconstruida es totalmente diferente a la original.
Actualización2:
Tengo la sospecha de que tal vez Resuelvo mal las ecuaciones diferenciales y eso provoca el error. Tengo el sistema de ecuaciones diferenciales (mi solucionador se puede encontrar en el RungeKutta4 y las funciones f arriba):
Fórmulas Frenet-Serret :
e'(t) = K(t)*velocity(t)*n(t)
n'(t) = -K(t)*velocity(t)*e(t) + T(t)*velocity(t)*b(t)
b'(t) = -T(t)*velocity(t)*n(t)
given the initial values:
e(0) = e0
n(0) = n0
b(0) = b0
where
K : curvature
T : torsionTal vez use mal la Runge-Kutta...
Actualización3:
Si cambio mi propio algoritmo Runge-Kutta por la función incorporada de Octave lsode y reescribo el código de integración numérica para utilizar trapz (regla del trapecio), nada cambia. Así que creo que el error no está en estas funciones... ¿Pero dónde?
