Matlab符号计算求导与化简

依托于Maple的符号计算引擎，Matlab可以完成高等数学中几乎一切符号计算。这里简单记录一下蔡自兴《机器人学》第三版中第四章的双连杆的动能推导过程。

在书本上我们已知的信息是：

$x_2 = d_1\sin\theta_1 + d_2\sin(\theta_1+ \theta_2)$
$y_2 = -d_1\cos\theta_1 - d_2 \cos(\theta_1 + \theta_2)$
$v_2^2 = \dot{x}^2 + \dot{y}^2$

其中 $\theta_i (i =1,2)$ 是时间的函数，我们通过演算可以得到与书中一样的不错的表达式：

$$v_2^2 = d_1^2\dot{\theta_1}^2 + d_2^2(\dot{\theta_2}^2+2\dot{\theta_1}\dot{\theta_2} + \dot{\theta_2}^2) + 2d_1d_2\cos\theta_2(\dot{\theta_1}^2 + \dot{\theta_1}\dot{\theta_2})$$

下面就以这个表达式为目标，探索如何利用Matlab进行带有导数的表达式的运算和化简

clc,clear,close all
syms th1(t) th2(t) d1 d2 Dth1 Dth2

x2 = d1*sin(th1) + d2*sin(th1 + th2);
y2 = -d1*cos(th1) - d2*cos(th1 + th2);
v22 = diff(x2)^2 + diff(y2)^2

v22 指的是 $v_2^2$ ，上述代码的输出结果为：

v22(t) =

(d2*cos(th1(t) + th2(t))*(diff(th1(t), t) + diff(th2(t), t)) + d1*cos(th1(t))*diff(th1(t), t))^2 + (d2*sin(th1(t) + th2(t))*(diff(th1(t), t) + diff(th2(t), t)) + d1*sin(th1(t))*diff(th1(t), t))^2

简化表达式

这个结果显然不是我们想要的，下面逐步进行变换化简

1. 用变量 Dth1 和 Dth2 替换 diff(th1(t), t) 和 diff(th2(t), t) ，简化表达式

   v22_pretty = subs(v22,[diff(th1(t), t) diff(th2(t), t)],[Dth1 Dth2])
   
   v22_pretty(t) =
   
   (d2*cos(th1(t) + th2(t))*(Dth1 + Dth2) + Dth1*d1*cos(th1(t)))^2 + (d2*sin(th1(t) + th2(t))*(Dth1 + Dth2) + Dth1*d1*sin(th1(t)))^2

2. 用 collect 函数按照 d1 的降幂次合并同类项

   v22_collect1 = collect(v22_pretty,d1)
   
   v22_collect1(t) =
   
   (Dth1^2*cos(th1(t))^2 + Dth1^2*sin(th1(t))^2)*d1^2 + (2*Dth1*d2*cos(th1(t))*cos(th1(t) + th2(t))*(Dth1 + Dth2) + 2*Dth1*d2*sin(th1(t))*sin(th1(t) + th2(t))*(Dth1 + Dth2))*d1 + d2^2*cos(th1(t) + th2(t))^2*(Dth1 + Dth2)^2 + d2^2*sin(th1(t) + th2(t))^2*(Dth1 + Dth2)^2

3. 表达式中出现了 $\cos^2\theta + \sin^2\theta$ 的表达式没有进行化简，下面用 simplify 函数进行化简

   v22_simplify = simplify(v22_collect1)
   
   v22_simplify(t) =
   
   Dth1^2*d1^2 + 2*cos(th2(t))*Dth1^2*d1*d2 + Dth1^2*d2^2 + 2*cos(th2(t))*Dth1*Dth2*d1*d2 + 2*Dth1*Dth2*d2^2 + Dth2^2*d2^2

4. 用 collect 依照 [d1 d2] 降幂次合并同类项

   v22_collect2 = collect(v22_simplify,[d1 d2])
   
   v22_collect2(t) =
   
   Dth1^2*d1^2 + (2*cos(th2(t))*Dth1^2 + 2*Dth2*cos(th2(t))*Dth1)*d1*d2 + (Dth1^2 + 2*Dth1*Dth2 + Dth2^2)*d2^2

5. 上面的表达式中，还有 $\cos\theta_2$ 没有化简

   v22_collect3 = collect(v22_collect2,cos(th2(t)))
   
   v22_collect3(t) =
   
   d1*d2*(2*Dth1^2 + 2*Dth2*Dth1)*cos(th2(t)) + (Dth1^2 + 2*Dth1*Dth2 + Dth2^2)*d2^2 + Dth1^2*d1^2

4 和5与下面的语句等价

collect(expr,[d1 d2 cos(th2(t))])

最后得到了文章伊始的表达式

后记

在已知结果情况下，省去探索的步骤，1-5的语句可以直接化简为：

v22_pretty = subs(v22,[diff(th1(t), t) diff(th2(t), t)],[Dth1 Dth2]);
v22_simplify = simplify(v22_pretty);
v22_collcet = collect(v22_simplify,[d1 d2 cos(th2(t))]);

依照一般的习惯，表达式的化简多依赖于 simplify 和 collect ，像类似 expand (展开) horner (嵌套) 等函数通常是我们不希望遇到的。