1. 引言
MATLAB具备强大的符号运算功能。符号运算就是所谓的计算机代数,通俗的说就是利用计算机进行数学公式的推导。这篇文章主要总结几个MATLAB进行符号运算时的小技巧,这也是作者在进行技术研究过程中实际碰到的一些难题,希望后来者能少走写弯路。
2. 符号运算
2.1 常规符号运算
先来说最朴素的情况,我们可以通过syms或者sym定义符号,我比较喜欢用syms,感觉这个指令功能更强大一些。代码如下:
syms v m g h;
L = 1/2*m*v^2 - mgh
以上代码定义了四个符号变量v,m,g,h。利用这四个符号变量定义了一个符号函数。此时希望对符号函数进行如下处理:
里面的偏导数是比较容易计算的,MATLAB提供了一个叫做diff的函数专门用来求导数。这样可以通过如下方式计算:
lv = diff(L, v)
此时得到的结果是,这个结果是符合我们预期的,但是接下来怎么让
对时间求导数呢?我们从公式上可以看出来
并不是显含时间
的,但是从原始公式上你可能已经看出
是一个速度的概念,它一定是和时间
相关的。也就是
。那么MATLAB如何处理此处对时间的导数呢?要不试一下:
syms t;
lvt = diff(lv, t);
很不幸,这种情况下你得到的将是0。因为MATLAB并不知道你这里的是时间的函数。这就是我们要介绍的第一个技巧。
2.2 抽象符号函数
前面我们提到MATLAB并不知道速度是时间的函数,那么怎么告知MATLAB呢?答案是在定义符号的时候将
定义成时间的抽象函数:
syms v(t) m g h;
L = 1/2*m*v^2 - mgh
syms v(t)这样的定义方式matlab会自动创建一个符号变量t并认为v是符号变量t的抽象函数。这样实现方式如下:
syms v(t) m g h;
L = 1/2*m*v^2 - mgh;
lv = diff(L, v);
lvt = diff(lv, t);
这时得到的lvt将是lvt=m*diff(v(t),t),这正是我们希望得到的。当公式变得越来越复杂时,你可能会觉得diff(v(t),t)这种表达方式过于繁琐,希望能够通过一个更简洁的变量来代替它,比如新定义一个符号a,用a取代diff(v(t),t),这就是第二个技巧。
2.3 变量替换
可以通过如下方式进行变量替换操作:
syms v(t) m g h;
L = 1/2*m*v^2 - mgh;
lv = diff(L, v);
lvt = diff(lv, t);
syms a;
lvt = subs(lvt, diff(v, t), a);
得到的表达式变为lvt=m*a。
2.4 符号表达式化简
当一个符号表达式存在冗余项时,可以对符号表达式进行合并等简化操作,使用的指令为simplify,代码如下:
syms v(t) m g h;
L = 1/2*m*v^2 - mgh;
lv = diff(L, v);
lvt = diff(lv, t);
L = simplify(L);
2.5 实变量符号运算
定义如下的符号表达式:
syms v(t) m g h;
y = [m g h];
z = y';
会发现z=[conj(m); conj(g); conj(h)]。出现这种问题的原因是matlab默认的运算 ' 是共轭转置,而单纯的转置使用的是 .' 因此一种解决方式是:
syms v(t) m g h;
y = [m g h];
z = y.';
另外一种方式是我们声明符号变量的时候就告知matlab这是一个实数范围内的符号变量:
syms v(t) m g h real;
y = [m g h];
z = y';
这种方式的问题是real本身只能修饰符号变量,不能修饰抽象符号函数v(t),因此matlab会报警,如果矩阵中含有v(t)那么进行转置操作v(t)依然会带上conj,解决方法是用assumeAlso指令:
syms m g h real;
syms v(t);
assumeAlso(v(t), 'real');
y = [m g h v];
z = y';
2.6 取出符号函数体
现在有一个新的问题是由抽象符号函数引入的,假如我们输入如下指令:
syms m g h real;
syms v(t);
assumeAlso(v(t), 'real');
y = [m g h v];
现在我们希望取出矩阵y中的前三个元素要怎么做呢?你可能自然而然会想到用y(1:3),但是这个操作的结果并不如你所愿,它得到的将会是一个四维的元组(cell),每个元组有三个元素,展开后是y={{m m m} {g g g} {h h h} {v(1) v(2) v(3)}},其实这个指令的真正含义是当t分别取1,2,3时矩阵各个维度的值。如果你真的需要得到矩阵y的前三个元素,需要进行如下操作:
syms m g h real;
syms v(t);
assumeAlso(v(t), 'real');
y = [m g h v];
y = formula(y);
y = y(1:3);
3. 参考文献
[1]. Issue with reality assumption of an implicit function -
[2]. How do you declare a symbolic function of time as a real variable -