Вычислить производную вектора

Question

У меня есть следующие функции (кривая Вивиани):

Phi  = @(t)[ cos(t)^2, cos(t)*sin(t), sin(t) ] 

Просто проверить, что это действительно:

s = linspace(0,T,1000); 
plot3(cos(s).^2, cos(s).*sin(s), sin(s)); 

Как дериват функцию Phi (возможно, несколько раз), который представляет кривую Вивиани в точке t, где t идет от 0 до 2*pi? Я определил Phi, подходящий для такого производного? Я пробовал diff, но он не оставил Phi, поскольку мне это нужно.

Если вторая производная будет Phi_d2, мне нужно получить ее значение (например, в t = 0).

Как я могу это достичь?

Answer 1

Вот три способа выполнить это. Первый использует subs, второй использует symfun, а третий использует complex step differentiation:

% Using subs 
syms t 
Phi = [cos(t) cos(t).*sin(t) sin(t)]; 
Phi_d2 = diff(Phi,t) 
double(subs(Phi_d2,t,0)) 

% Using symfun 
syms t 
Phi(t) = [cos(t) cos(t).*sin(t) sin(t)]; 
Phi_d2 = diff(Phi,t) 
double(Phi_d2(0)) 

% Using complex step differentiation 
Phi = @(t)[cos(t) cos(t).*sin(t) sin(t)]; 
h = 2^-28; 
cdiff = @(f,x)imag(f(x(:)+1i*h))/h; 
Phi_d2 = cdiff(Phi,0) 

Вы можете найти функцию для выполнения первого и второго порядка, сложный этап дифференциации on my GitHub: cdiff. Обратите внимание, что комплексная дифференциация шагов не будет работать хорошо для производных более высокого порядка. Лучше всего, когда у одного есть недифференцируемая функция или нужны быстрые численные первые производные.

Answer 2

Ради полноты картины, численное решение без использования каких-либо дополнительных инструментарии:

N = 999; 
t = linspace(0,2*pi,N+1); 
Phi = [cos(t); cos(t).*sin(t); sin(t)]; 
dPhi = gradient(Phi,2*pi/N) 

Для неравномерно разнесенных векторов аргумента, второй параметр gradient определяется промежутком вектора вместо скаляр. (подходящий вектор времени или угла) - в этом случае, очевидно, необходимо разделить размеры. (Хотя я не знаю, почему.)

Так в качестве альтернативы, хотя и не нужно:

dX = gradient(Phi(1,:),t); 
dY = gradient(Phi(2,:),t); 
dZ = gradient(Phi(3,:),t); 
dPhi = [dX; dY; dZ];