Matlab符号数学工具箱使用指南：符号运算与数值计算的完美结合
Matlab符号数学工具箱使用指南：符号运算与数值计算的完美结合
符号运算与数值计算
缘，妙不可言
在使用Matlab进行代码开发时，经常遇到将符号运算与数值计算混合使用的情况。从概念上看，这种做法很有吸引力。开发者通常会先使用syms定义符号变量，将公式输入到符号领域进行推导，最后代入数值得到结果。
然而，这种混合使用的方式可能会带来一些问题。虽然Matlab支持这种工作流程，但为了确保计算效率，建议将符号计算和数值计算分开处理。
建议的实践
结论先行：将符号计算和数值计算分开处理。
推荐的工作流程是：将符号推导部分专门放在一个LiveScript文件中，推导的结果生成一系列函数文件。这些函数文件就是不同的Matlab函数，后续可以进行正常的数值计算工程应用。
理由
主要理由是：符号运算中的符号会污染所有的计算，将一切涉及符号的表达式变成符号变量和表达式，这会使得Matlab在矩阵数值计算方面的优势无法充分发挥。符号计算的速度通常比数值计算慢得多。
我们举一个简单的例子：
syms f(x,y,z)
f(x,y,z) = y*z*sin(x) + x*sin(z)*cos(y) - z^3;
[xDouble,yDouble,zDouble] = meshgrid(1:20,1:50,1:20);
接下来比较符号计算和数值计算的时间：
>> timeit(@()f(xDouble,yDouble,zDouble), 1)
ans =
2.0569
上述代码中，timeit函数用于评估fDouble = f(xDouble,yDouble,zDouble)的时间。现在来看数值计算的时间：
>> f1 = matlabFunction(f)
f1 =
包含以下值的 function_handle:
@(x,y,z)-z.^3+y.*z.*sin(x)+x.*cos(y).*sin(z)
>> timeit(@()f1(xDouble,yDouble,zDouble), 1)
ans =
9.3786e-04
可以看出，时间差距超过3个数量级。即使使用点运算符，结果也类似：
>> syms f3(x,y,z)
>> f3(x,y,z) = y.*z.*sin(x) + x.*sin(z).*cos(y) - z.^3;
>> timeit(@()f3(xDouble,yDouble,zDouble), 1)
ans =
2.1108
综上所述，不要将任何符号带入大规模的数值计算中，比如大范围的网格计算或循环数据处理。
接下来，介绍如何将符号推导的结果转化为数值计算中可以高效使用的函数。
符号计算工具箱
在前面的时间评估例子中，使用了matlabFunction函数，该函数是2008b版本之后引入的符号计算工具箱的一部分。符号计算工具箱主要处理以下对象：
- 符号常量：准确表达的数值
- 符号变量：数学上的代数变量
- 符号表达式：符号表达式
符号常量和符号变量
通过sym函数可以定义符号常量和符号变量：
a = sym(1/3);
piSym = sym(pi);
区别在于：
>> sin(piSym)
ans =
0
>> sin(pi)
ans =
1.2246e-16
可以通过vpa函数将浮点数转化为符号常量，该函数默认截取32位有效数字，也可以通过第二个参数指定有效数字位数：
>> vpa(pi, 100)
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
新建符号变量可以使用sym函数的两种调用方式：
>> sym('x', [3, 3])
ans =
[x1_1, x1_2, x1_3]
[x2_1, x2_2, x2_3]
[x3_1, x3_2, x3_3]
或者使用语法糖，直接列出符号变量：
>> syms x y z [3 3]
这样可以创建三个符号变量矩阵，矩阵的元素都是符号变量，约定向量的元素是x1,x2,x3,...，矩阵的元素是x1_1,x1_2,x1_3,...。
符号表达式
符号表达式是符号变量的组合，可以通过涉及符号变量和符号常量的运算得到：
syms x y z
f = x^2 + y^2 + z^2;
这样定义的f是一个sym类型的对象，可以进行各种运算。也可以使用symfun函数定义一个函数：
syms x y z
g = symfun(x^2 + y^2 + z^2, [x, y, z]);
或者更简洁的方式：
syms h(x,y,z)
h(x,y,z) = x^2 + y^2 + z^2;
sym和symfun的区别
sym定义的是一个符号表达式，而symfun定义的是一个符号函数。实际使用中，symfun可以像函数一样调用，将形式符号替换为实参对应的符号表达式，最终结果是一个符号表达式（sym）。在将符号表达式转换为函数句柄时，无论是sym还是symfun都是可以的。从继承关系来看，symfun是sym的子类。
符号表达式的转化
基础调用方式：转化为匿名函数
matlabFunction函数有两种调用语法：
ht = matlabFunction(f)
ht = matlabFunction(f1,...,fN)
可以是一个sym对象，也可以是多个sym对象。前一种情况得到一个输出量个数为1的函数句柄，后一种情况下得到一个输出量个数为N的函数句柄。
示例：
>> syms x y z
f = x^2 + y^2 + z^2;
>>
>> ht = matlabFunction(f)
ht =
包含以下值的 function_handle:
@(x,y,z)x.^2+y.^2+z.^2
多输出函数句柄的例子：
>> syms g(a, b, c)
>> g(a,b,c) = a + b + c
g(a, b, c) =
a + b + c
>> ht = matlabFunction(f, g)
ht =
包含以下值的 function_handle:
@(a,b,c,x,y,z)deal(x.^2+y.^2+z.^2,a+b+c)
使用示例：
>> ht(1,2,3,4,5,6)
错误使用 deal
输入的数目应与输出的数目匹配。
出错 symengine>@(a,b,c,x,y,z)deal(x.^2+y.^2+z.^2,a+b+c)
>> [x1, x2] = ht(1,2,3,4,5,6)
x1 =
77
x2 =
6
进阶调用：转化为函数文件
matlabFunction函数还可以将符号表达式转化为函数文件：
matlabFunction(f, 'File', 'myfun1');
matlabFunction(f, g, 'File', 'myfun2');
生成的myfun1.m文件示例：
function f = myfun1(x,y,z)
%MYFUN1
% F = MYFUN1(X,Y,Z)
% This function was generated by the Symbolic Math Toolbox version 23.2.
% 2024-10-19 23:19:02
f = x.^2+y.^2+z.^2;
end
生成的myfun2.m文件示例：
function [f,g] = myfun2(a,b,c,x,y,z)
%MYFUN2
% [F,G] = MYFUN2(A,B,C,X,Y,Z)
% This function was generated by the Symbolic Math Toolbox version 23.2.
% 2024-10-19 23:19:02
f = x.^2+y.^2+z.^2;
if nargout > 1
g = a+b+c;
end
end
实际应用
考虑一个分段函数，需要推导其导数并在数值计算中同时求得其值和导数：
syms p(x)
p(x) = piecewise(x<0, x^2-8, x>=0, -x)
函数图像：
fplot(p, [-10, 10])
推导并输出函数文件：
dp = diff(p, x);
matlabFunction(p, dp, 'File', 'piecewiseFunc');
生成的piecewiseFunc.m文件：
function [p,dp] = piecewiseFunc(x)
%piecewiseFunc
% [P,DP] = piecewiseFunc(X)
% This function was generated by the Symbolic Math Toolbox version 23.2.
% 2024-10-19 23:42:41
t2 = (x < 0.0);
if ~all(cellfun(@isscalar,{t2,x}))
error(message('symbolic:sym:matlabFunction:ConditionsMustBeScalar'));
end
if (t2)
p = x.^2-8.0;
elseif (0.0 <= x)
p = -x;
else
p = NaN;
end
if nargout > 1
if ~all(cellfun(@isscalar,{t2,x})) error(message('symbolic:sym:matlabFunction:ConditionsMustBeScalar'));end;if (t2) dp = x.*2.0;elseif (0.0 < x) dp = -1.0;else dp = NaN;end;
end
end
由于这个函数不能直接处理向量，可以再封装一层：
function [p, dp] = piecewiseFuncVec(x)
[p, dp] = arrayfun(@piecewiseFunc, x);
end
使用示例：
x = -10:0.1:10;
[p, dp] = piecewiseFuncVec(x);
tiledlayout(2,1)
nexttile
plot(x, p)
nexttile
plot(x, dp)
这样就可以得到函数图像和导数图像。
总结
- 一定要分离符号计算和数值计算，以免符号计算的低效影响数值计算的效率。
fplot函数可以直接绘制符号表达式的图像。- 符号函数是符号表达式的一种，可以直接调用，定义方式为
syms f(x,y,z); f(x, y, z)=...。 - 符号表达式可以通过
matlabFunction函数转化为函数句柄或函数文件，这样就可以在数值计算中使用。