How do I change variables so that I can differentiate with respect to a derivative?

Related question: https://www.mathworks.com/matlabcentral/answers/105526-differentiation-with-respect-to-a-derivative

I am attempting something similar to the question above. I attached my code below. If I use dq1(t) = diff(q1(t), t), then I get an error on D3(t). If I define dq1(t) as a symbol, I get the error shown below. I want to differentiate a function with respect to a derivative, and then differentiate that function with respect to a variable that the derivative depends on.

% Max 3 Dof
% No Non-conservative forces
clear all; clc; close all;

% Symbols
syms q1(t) q2(t) dq1(t) dq2(t) y1 y2 m1 m2 g

% Kinetic Energy
% dq1(t) = diff(q1(t), t)
% dq2(t) = diff(q2(t), t)

T1(t) = (m1/2)*(dq1(t)^2)

T1(t) = 

T2(t) = (m2/2)*(dq2(t)^2)

T2(t) = 

T = T1 + T2

T(t) = 

% Potential Energy
V1 = m1*g*y1 % gravitational

V1 = 

V2 = m2*g*y2 % gravitational

V2 = 

V = V1 + V2

V = 

% Lagrangian Function
L = T - V

L(t) = 

% Partial Derivatives
D1(t) = functionalDerivative(L, q1(t))

D1(t) = 

0

D2(t) = functionalDerivative(L, q2(t))

D2(t) = 

0

D3(t) = functionalDerivative(L, dq1(t))

D3(t) = 

D4(t) = functionalDerivative(L, dq2(t))

D4(t) = 

% Time derivatives
D5(t) = functionalDerivative(D3(t), t)

Error using sym/functionalDerivative
Variables must be specified by function calls.

D6(t) = functionalDerivative(D4(t), t)

% Lagrangian Equation of Motion
accel1 = functionalDerivative(q1, t, 2)
accel2 = functionalDerivative(q2, t, 2)

eqn1 = simplify(D5 - D1);
eqn2 = simplify(D6 - D2);

% Solutions
q1Accel = solve(eqn1, accel1);
q2Accel = solve(eqn2, accel2);

[Bjorn Gustavsson]

My advice would be to look at the Euler-Lagrange tools on the file exchange. There are a couple of good submissions there that does just this. You can read the code and use the same tricks - if I recall correctly it is a bit fidgety.

[Curran Robertson]

Thank you for that. There are some good examples there.

Answer 1

[Paul]

Hi Curran,

The doc page functionalDerivative has related example and shows the problem can be solved like this:

% Symbols
syms q1(t) q2(t) dq1(t) dq2(t) y1 y2 m1 m2 g

% Kinetic Energy
dq1(t) = diff(q1(t), t);
dq2(t) = diff(q2(t), t);
T1(t) = (m1/2)*(dq1(t)^2);
T2(t) = (m2/2)*(dq2(t)^2);
T = T1 + T2

T(t) = 

I changed y to q in the potential energy, which seems right. Why did the original code use y?

% Potential Energy
% V1 = m1*g*y1 % gravitational
% V2 = m2*g*y2 % gravitational
% change y to q
V1 = m1*g*q1(t); % gravitational
V2 = m2*g*q2(t); % gravitational
V = V1 + V2

V = 

% Lagrangian Function
L = T - V

L(t) = 

Here, we just take the functional derivatives wrt q1(t) and q2(t)

% Lagrangian Equation of Motion
eqn = functionalDerivative(L,[q1(t) q2(t)]) == 0

eqn(t) = 

This result seems correct and can be simplified further

assumeAlso([m1 m2],'positive');
eqn = simplify(eqn)

eqn(t) = 

If you want to do the whole thing by hand, then this approach seems correct for this problem

syms q1(t) q2(t) dq1(t) dq2(t) y1 y2 m1 m2 g

% Kinetic Energy
%dq1(t) = diff(q1(t), t)
%dq2(t) = diff(q2(t), t)
T1(t) = (m1/2)*(dq1(t)^2);
T2(t) = (m2/2)*(dq2(t)^2);
T = T1 + T2

T(t) = 

% Potential Energy
% change y to q
V1 = m1*g*q1(t); % gravitational
V2 = m2*g*q2(t); % gravitational
V = V1 + V2

V = 

% Lagrangian Function
L = T - V

L(t) = 

D1(t) = functionalDerivative(L, q1(t));
D2(t) = functionalDerivative(L, q2(t));
D3(t) = functionalDerivative(L, dq1(t));
D4(t) = functionalDerivative(L, dq2(t));

Now, in order to get the correct time derivatives, we need to sub back in the proper definition of dq1(t) and dq2(t)

% Time derivatives
D5(t) = diff(subs(D3(t),dq1,diff(q1,t)), t);
D6(t) = diff(subs(D4(t),dq2,diff(q2,t)), t);
% Lagrangian Equation of Motion
eqn = [simplify(D5 - D1); simplify(D6 - D2)] == 0

eqn(t) = 

assumeAlso([m1 m2],'positive')
eqn = simplify(eqn)

eqn(t) = 

[Curran Robertson]

This seems to be the correct solution to the question I asked. The reason I used y1 and y2 is due to the physics of the problem. The potential energy is related to the height of the object. q1 and q2, the degrees of freedom, are not necessarily y1 and y2. Using the subs function does seem to work, but I think I will have to modify the code to hard code q1, q2, y1, and y2 because they are dependent on each other.

[Paul]

I think that either approach shown above would work with y1 and y2 instead of q1 and q2 for the potential energy. But, if y1 and y2 depend on q1 and q2, you'll have to make those substitutions early in the process. Probably best to make those subs right up front. The symbolic engine won't care.

Good luck with your project.