Matlab 符号计算和任意精度计算

贡献者： addis

　　 Matlab 的符号计算需要符号计算工具箱，取决于你的证书类型，可能需要额外购买。另外需要提醒的是，虽然 Matlab 和 Python 都有符号计算功能，但在符号计算领域 Mathemtica 才更为主流，其默认界面也更适合符号计算。

1. 符号变量和符号表达式

• Matlab 中用于储存符号计算表达式的变量类型为 sym。可以用 syms 变量1 变量2 ... 声明变量类型为 sym。例如 syms x y z;。
• syms f(x) 定义一个符号函数 f。该命令会同时创建符号变量 x 和符号函数 f(x)（symfun 类型，没有具体定义）。接下来可以用 f(x) = x^2 对其定义，这相当于运行 syms x; f(x) = x^2;。
• 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')。注意在较新版的 Matlab 中，sym('pi') 将得到名为 pi 的符号变量，而不是圆周率。
• 对表达式求导如 diff(sym('x')^2)。若要求偏导，diff(sym('x')^2 * sym('y')^3) 默认对 x 求偏导（x 是通过 symvar(expr, 1) 决定的），结果是 $2x y^3$。也可以声明对 y 求偏导，如 diff(sym('x')^2 * sym('y')^3, sym('y'))，或者更简洁地，diff(sym('x')^2 * sym('y')^3, 'y')。这是因为，如果若函数的一个形参需要符号表达式，但实参使用了字符串或者数值，那么 Matlab 就会自动将其用 sym() 函数进行转换。建议总是声明对哪个变量求偏导。若不声明，则会通过 symvar(expr, 1) 决定。expr 也可以是 symfun 类型的。
• diff(f, var, n) 可以求高阶导数。diff(f, y, z) 求混合偏导。
• 要把一个数字作为符号，可以使用例如 sym('2') 或 sym('7/3')，但这仅限于分式，不允许诸如 sym('sqrt(2)') 这样的用法，应该用 str2sym('sqrt(2)') 或者 sqrt(sym(2))。
• 另一种方法是使用 sym(数值)。例如 sqrt(sym(2)) 的结果是表达式 $\sqrt 2$，而不同于数值计算中的 sqrt(2) = 1.414...。由于符号计算是精确无误差的，无理数 $\sqrt{2}$ 并不会被自动转换为小数形式。
• sym(数值) 会自动猜测 数值 所代表的根式，例如 sym(0.866025403784439) 的结果是表达式 $\sqrt{3}/2$，又例如 sym(pi) 的结果是精确的圆周率符号 $\pi$。但 sym(sqrt(pi)) 就只能猜出一个近似的分数。
• 如果不希望猜测，那么可以用 “双精度和变精度浮点数测试（Matlab）” 中的 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} $。

2. 变精度计算

　　 注意变精度计算功能往往和符号计算一同使用，但这两个功能从实现来看是完全不同的。变精度计算和双精度一样本质上还是数值浮点计算，计算过程存在数值误差。只不过我们可以规定每个变量的有效数字为任意多，而不是统一使用双精度类型的 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$ 该数无论参与任何计算都是绝对精确的。
• 变精度浮点数（vpa）的计算并不会追踪误差，例如 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) 结果相同。
• 如果不希望猜测，那么可以用 “双精度和变精度浮点数测试（Matlab）” 中的 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$。

3. 特殊函数

　　 由于一些特殊函数若用双精度计算，往往在一些区间产生较大误差。所以 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) 来求任意位有效数字。