贡献者: addis
Matlab 的符号计算需要符号计算工具箱,取决于你的证书类型,可能需要额外购买。另外需要提醒的是,虽然 Matlab 和 Python 都有符号计算功能,但在符号计算领域 Mathemtica 才更为主流,其默认界面也更适合符号计算。
sym
。可以用 syms 变量1 变量2 ...
声明变量类型为 sym
。例如 syms x y z;
。Matlab 的大部分自带算符和函数支持 sym
类型的变量,例如 x^2
就是 sym
类型的表达式 $x^2$,$x$ 并不是一个数值而是符号。若此时令 expr = x^2
,那么用 class(expr)
可以验证 expr
的类型也是 sym
。
syms
一次声明几个符号变量,也可以使用 sym(字符串)
其中字符串只能是变量名或数字。例如 syms x; expr = x^2;
得到的 expr
和 expr = sym('x')^2;
得到的 expr
是等效的。
str2sym(字符串)
可以把字符串转换为符号表达式。例如 str2sym('pi')
或 str2sym('sqrt(3)/2')
。
diff(sym('x')^2)
。若要求偏导,diff(sym('x')^2 * sym('y')^3)
默认对 x
求偏导,结果是 $2x y^3$。也可以声明对 y
求偏导,如 diff(sym('x')^2 * sym('y')^3, sym('y'))
,或者更简洁地,diff(sym('x')^2 * sym('y')^3, 'y')
。这是因为,如果函数的一个参数需要符号表达式,但如果输入时使用了字符串或者数值,那么 Matlab 就会自动将其用 sym()
函数进行转换。建议总是声明对哪个变量求偏导。
sym('2')
或 sym('7/3')
,但这仅限于分式,不允许诸如 sym('sqrt(2)')
这样的用法,应该用 str2sym('sqrt(2)')
。
sym(数值)
。例如 sqrt(sym(2))
的结果是表达式 $\sqrt 2$,而不同于数值计算中的 sqrt(2) = 1.414...
。由于符号计算是精确无误差的,无理数 $\sqrt{2}$ 并不会被自动转换为小数形式。另一个例子,d = 3.1; sqrt(sym(d))
得到的是表达式 $961/100$。
sym(数值)
会自动猜测 数值
所代表的根式,例如 sym(0.866025403784439)
的结果是表达式 $\sqrt{3}/2$,又例如 sym(pi)
的结果是精确的圆周率符号 $\pi$。注意在较新版的 Matlab 中,sym('pi')
将得到名为 pi
的普通变量,而不是圆周率。
num2sym(双精度)
。
sym(...)
作用在一个数组上(可以是任何上述类型),那么则逐个元素作用,并生成 sym
类型的数组。
sym
类型和一个 double 类型进行 +, -, *, /, ^
等运算时,double
类型的数会自动被 sym()
函数转换为 sym
类型。例如 sym(1)/3
和 sym(1)/sym(3)
是等效的。这样可以让输入更简洁。
subs(符号表达式,符号变量, 新表达式)
可以把表达式中的所有 符号变量
替换为 新表达式
。例如 subs(sin(sym('x')), 'x', 'y')
相当于 subs(sin(sym('x')), sym('x'), sym('y'))
,结果是 $\sin y$。又例如 subs(sin(sym('x')), 'x', pi/4)
的结果是 $\sqrt 2/2$。
新表达式
是一个数组,则依次对每个元素进行替换,输入一个 sym
类型的数组。
double()
,例如 double(sqrt(sym(2)))
结果是 1.414...
,是一个双精度数值,误差就是双精度类型的相对误差 eps
,约为 $2.2 \times 10^{-16} $。
注意变精度计算功能往往和符号计算一同使用,但这两个功能从实现来看是完全不同的。变精度计算和双精度一样本质上还是数值浮点计算,计算过程存在数值误差。只不过我们可以规定每个变量的有效数字为任意多,而不是统一使用双精度类型的 15-16 位有效数字。
Matlab 的变精度计算并不自带误差追踪功能,例如两个十分相近的数相减,返回的小数位中可能大部分都是错的。而 Mathematica 的做法是返回更少小数位,但保证最后一位是对的。
double(符号表达式)
返回双精度结果,vpa(符号表达式, 位数)
可以计算符号表达式的任意位有效数字结果,例如 vpa(sqrt(sym('2')), 50)
计算 $\sqrt{2}$ 的 50 位有效数字,并能保证四舍五入后最后一位有效数字正确。返回的值是一个 'sym'
类型的对象而不是 'double'
类型。
vpa
返回的类型也是 'sym'
,但它是有误差的,本质上是一个类似 double
的浮点类型。vpa()
就是 variable precision arithmetic,变精度计算。我们姑且把这种有误差的 sym
数字称为变精度浮点数。我们可以把变精度浮点数当作 double
类型的拓展。double
在十进制下只有约 15-16 位有效数字,而变精度浮点数的有效数字位数可以任意指定。
sym
类型的对象 x
是否有误差,只需要在命令行中把它显示出来(可以使用 disp()
函数,也可以直接运行 x
不加分号)。如果显示中出现了小数点,那么他就是变精度浮点数。sym
类型的整数,无论有多少位,如果没有误差则会把每一位显示出来而不是用科学计数法(因为科学计数法带小数点,表示有误差)。
vpa()
函数会输出变精度浮点数,另一种方法如 sym('12.3')
或者 sym('1.23e1')
(用 vpa
也一样),运行后显示结果为 12.3
,存在小数点,所以是变精度浮点数。相比之下,sym('123/10')
结果显示为 $123/10$ 该数无论参与任何计算都是绝对精确的。
sym('1.2345678902234567890323456789042345678905234567890') - sym('1.2345678902234567890323456789042345678905234500000')
返回的是 6.789000000000011835767836416814e-45
而不是 6.789e-45
。而 Mathematica 就会返回 6.789e-45
。
vpa()
和 sym()
中如果有效数字超过 15-16,要用引号,否则会先转换为 double,超出 15-16 位的部分丢失。
vpa(数值)
和 sym(数值)
一样会自动猜测 数值
所代表的根式或符号,结果相当于 vpa(sym(数值))
。例如 vpa(0.866025403784439)
的结果是表达式 $\sqrt{3}/2$,又例如 sym(pi)
的结果是精确的圆周率符号 $\pi$。例如 vpa(0.866025403784439) - vpa(str2sym('sqrt(3)/2'))
结果为零,又例如 vpa(pi, 1000)
可以得到圆周率的 1000 位有效数字。
vpa
猜测,则应该使用引号,例如 vpa('0.866025403784439') - vpa(str2sym('sqrt(3)/2'))
不为零。使用引号可以在当前 digits
设置的精度内最精确地表示引号内的数字。验证:vpa('1.866025403784439') - vpa('1.866025403784438')
vpa(1+1e-14)
得 1.0
,sym(1+1e-14)
得 1
。sym(pi+1e-14)
得 pi
。这有时候可能反而会造成麻烦。
vpa(双精度)
会先把 双精度
变为 ieee 二进制表示,然后在后面添零拓展精度。验证:num2bin2(vpa(1.234567890223456789))
和 num2bin2(1.234567890223456789)
结果相同。
num2vpa(双精度, 有效数字)
。
gamma
,hypergeom
)会保证输出结果的最后一位正确(可以对比 Mathematica 结果),但由于上一条中误差,自己写的函数则无法保证。
hypergeom
只有符号工具箱中有,所以参数直接输入 double 就行,gamma
函数有普通的版本(只支持 double 实数),如果输入 double 复数会出错,这是只需要用 gamma(vpa(复数))
即可。
1 + sym('1e-40') - 1
的结果显示为 0.0
,而 1 + sym(10)^(-40) - 1
的结果显示 1/10000....
(40 个 0)。
digits(整数)
设置,如果不设置,则默认是 32 位。也可以使用不含自变量的 digits()
查看当前设置的位数。实际上的有效位数比设置的要多 8 位,所以默认是 40 位。这可以用 “双精度和变精度浮点数测试(Matlab)” 中的 digits2
验证。
digits
函数是控制全局的,如果设置了 digits
以后调用某函数,那么函数内部的变精度计算也都会采用同样精度。
digits(32)
后,vpa(sym(pi)/2, 50) + vpa(sym(pi)/2, 100)
的结果仍然是 32 位有效数字。
digits
以后,字符串 vpa 如 vpa('1.2345...')
相当于 vpa('1.2345...', 有效数字)
(实际上还要多 8 位,也就是最后一位后面有 8 个零),若两个精度高于 digits
设置的 vpa('数字字符串')
进行运算,就先把高精度的有效数字减少为低精度的再进行运算。但如果其中一个 vpa('数字字符串')
精度低于 digits
的设置,就先变为 digits
设置的精度。例子:先设置 digits(32)
,则 vpa('1.2345678902234567890323456789042345678905234567890623456789') - vpa('1.23456789022345678903234567890423456')
的结果有 8 位有效数字。另一个例子 vpa('1.2345678902234567890323456789') - vpa('1.2345678902234567890323456788')
结果有 12 位有效数字。
digits(50)
,再运行 (1 + sym('1e-40')) - 1
就会得到 0.000...000999...99938892...
,这就比默认的 32 位有效数字计算精确多了,但还是存在误差。
sym('1.2')*sqrt(2*sym('x'))
,结果是 $1.2 \sqrt{2x}$,其中 $1.2$ 是变精度浮点数。相比之下,sym(12/10)*sqrt(2*sym('x'))
得到精确的 $6\sqrt{2x}/5$。
由于一些特殊函数若用双精度计算,往往在一些区间产生较大误差。所以 Matlab 决定只使用任意精度来计算。以复数域的 $\Gamma$ 函数 为例,若直接输入 gamma(1+1i)
则会出错,因为 Matlab 中非符号计算版本的 gamma
函数只支持 double
类型的实数输入。所以要调用符号计算工具箱提供的 gamma()
,就输入一个 sym
类型的变量即可。例如 gamma(vpa(1+1i))
返回 0.49801... - 0.15494...i
,又例如 gamma(sym(1+1i))
返回表达式 gamma(1 + 1i)
。这是因为 vpa(1+1i)
是变精度浮点数,可以显示为小数,而 sym(1+1i)
不存在误差,不能显示为小数。当然,我们也可以用例如 vpa(gamma(sym(1+1i)), 50)
来求任意位有效数字。