MATLAB符号运算小技巧

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。利用这四个符号变量定义了一个符号函数。此时希望对符号函数进行如下处理:

\frac{d}{dt}\left( \frac{\partial L}{v} \right )

里面的偏导数是比较容易计算的,MATLAB提供了一个叫做diff的函数专门用来求导数。这样\frac{\partial L}{v}可以通过如下方式计算:

lv = diff(L, v)

此时得到的结果是lv=m*v,这个结果是符合我们预期的,但是接下来怎么让lv对时间求导数呢?我们从公式上可以看出来lv并不是显含时间t的,但是从原始公式上你可能已经看出v是一个速度的概念,它一定是和时间t相关的。也就是v=v\left(t \right )。那么MATLAB如何处理此处对时间的导数呢?要不试一下:

syms t;
lvt = diff(lv, t);

 很不幸,这种情况下你得到的将是0。因为MATLAB并不知道你这里的v是时间的函数。这就是我们要介绍的第一个技巧。

2.2 抽象符号函数

        前面我们提到MATLAB并不知道速度v是时间的函数,那么怎么告知MATLAB呢?答案是在定义符号的时候将v定义成时间的抽象函数:

syms v(t) m g h;
L = 1/2*m*v^2 - mgh

syms v(t)这样的定义方式matlab会自动创建一个符号变量t并认为v是符号变量t的抽象函数。这样\frac{d}{dt}\left( \frac{\partial L}{v} \right )实现方式如下:

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 -

[3]. Return body of symbolic function - MATLAB formula

[4]. Use a symfun as index of a matrix -

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:大白 设计师:CSDN官方博客 返回首页
评论 1

打赏作者

hitgavin

你的鼓励是我创作的动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值