贡献者: 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) 来求任意位有效数字.