Problem using subs function for multiple variables on 2023b

Hi, I have the following code:

syms ml1 ml2 a1 a2 q1 q2 q1p q2p q1pp q2pp l1 l2 g real
syms theta1(t) theta2(t) 
type_sym= [q1 q1p q1pp q2 q2p q2pp];
type_symfun= [...
    theta1(t) diff(theta1(t),t) diff(theta1(t),t,t) ...
    theta2(t) diff(theta2,t) diff(theta2,t,2)];

X0p =[(a1^2*cos(theta1(t))^2 + a1^2*sin(theta1(t))^2)*diff(theta1(t), t, t) ; 0] ;
X0p=subs(X0p,type_symfun,type_sym);

disp(X0p)

And I obtain:

 [0 ;0]

which is wrong. I its happening on Matlab 2023b version. However, on previous versions, I didn’t get this problem. Any idea?

[Dyuman Joshi]

"which is wrong."

It is correct.

When you substitute q1 for theta1(t), the expression changes to

a1^2*diff(q1, t, 2)

And as you q1 is not a function of t, it is treated as a constant, thus it's 2nd derivative w.r.t t is zero.

[Jose]

I have just tried on 2022a and I obtain the appropiate result:

q1pp*(a1^2*cos(q1)^2 + a1^2*sin(q1)^2)

Answer 1

[Walter Roberson]

Your theta2 is a symfun so diff(theta2,t) is a symfun, and that makes the entire type_symfun expression into a symfun that returns an array with 6 components. The meaning of substituting a vector of variables for a single symfun is not well defined. You should either use {} around the components or else you should treat theta2 the same way you treated theta1, putting the (t) in explicitly. Or both.

syms ml1 ml2 a1 a2 q1 q2 q1p q2p q1pp q2pp l1 l2 g real
syms theta1(t) theta2(t) 
type_sym= {q1 q1p q1pp q2 q2p q2pp};
type_symfun= {...
    theta1(t) diff(theta1(t),t) diff(theta1(t),t,t) ...
    theta2(t) diff(theta2(t),t) diff(theta2(t),t,2)};

X0p =[(a1^2*cos(theta1(t))^2 + a1^2*sin(theta1(t))^2)*diff(theta1(t), t, t) ; 0] ;
X0p=subs(X0p,type_symfun,type_sym);

disp(X0p)

As @Dyuman Joshi indicates, this 0 result is correct. After you substitute a symbolic name for theta1(t) then you are taking the derivative of the name with respect to t, and that is going to be 0 unless the name happens to be t.

syms ml1 ml2 a1 a2  l1 l2 g real
syms q1(t) q2(t) q1p(t) q2p(t) q1pp(t) q2pp(t)
syms theta1(t) theta2(t) 
type_sym= {q1 q1p q1pp q2 q2p q2pp};
type_symfun= {...
    theta1(t) diff(theta1(t),t) diff(theta1(t),t,t) ...
    theta2(t) diff(theta2(t),t) diff(theta2(t),t,2)};

X0p =[(a1^2*cos(theta1(t))^2 + a1^2*sin(theta1(t))^2)*diff(theta1(t), t, t) ; 0] ;
X0p=subs(X0p,type_symfun,type_sym);

simplify(X0p)

ans = 

The appearance of the second derivative here has to do with the order of substitutions: the diff() operations are active operations so after q1 gets substituted for theta1, derivatives are calculated. Ideally you would want the substitution of q1pp to happen "after" the derivatives are taken -- but subs() has never guaranteed the order of operations. Internally it makes the most sense for the substitutions to be done in a single pass. For example,

syms x y
subs(x+y, {x, y}, {2*y, 3*x})

ans = 

If you did the substitution for x first and evaluated, you would get 2*y+y which would simplify to 3*y and then you would subs 3*x for each y getting 9*x as the result. Likewise if you did the substitution for y first you would get a result in a single variable. In order to be able to "swap" variables like this, you have to do a single pass without evaluating and then let it evaluate once afterwards.

[Dyuman Joshi]

Walter, I guess the question is that why is the result in R2022a different from the result in R2023b (for the same code).

[Jose]

Exactly, because I have never had these kinds of problems until I updated to 2023b...

[Walter Roberson]

The question gets down to the order of substitutions. You have to substitute for the diff() before you substitute for theta1().

The order of substitutions is not documented.

It would not be unreasonable to substitute for the "largest" sub-tree first, but it looks like R2023b might now process in order.

syms ml1 ml2 a1 a2 q1 q2 q1p q2p q1pp q2pp l1 l2 g real
syms theta1(t) theta2(t) 
type_sym= {q1 q1p q1pp q2 q2p q2pp};
type_symfun= {...
    theta1(t) diff(theta1(t),t) diff(theta1(t),t,t) ...
    theta2(t) diff(theta2(t),t) diff(theta2(t),t,2)};

X0p =[(a1^2*cos(theta1(t))^2 + a1^2*sin(theta1(t))^2)*diff(theta1(t), t, t) ; 0] ;
X0p=subs(X0p, fliplr(type_symfun), fliplr(type_sym));

disp(X0p)

[Jose]

OK, now it works, thank you!