将符号表达式转换为函数时如何避免精度损失?
How to avoid loss of accuracy when converting symbolic expressions to functions?
tl;dr: 我观察到在使用 matlabFunction
将符号表达式转换为函数句柄时完全失去了准确性。不知道有没有什么办法可以提高转换率,避免精度损失。
我有一个符号变量 x
,1 <= x && x <= 2
成立:
syms x real
assumeAlso(1 <= x & x <= 2);
我在该变量中有许多计算机生成的冗长符号表达式。其中一个看起来像这样:
expression = ( ...
1015424780204960143215273323910078528628754663952913658657288835138029704791686717561885487746105223164496264397062144000 ...
*( ...
60345244216851610523130575942127473698515638085026410114070496399315754724985170479034760688327679099512793302302720*2^(1/2) ...
- 85341062796188128379389251264141456937242003828816791711306294886439805159655912635773490018204138847163921224639041 ...
)*(x - 2)^2 ...
) ...
/31504288346872372712061562812941419427167561153216213605146864117586323331155800294021400857537041202039565879856190821729945216935526707438840242294801134287960934949194864204307116208568439380511702602881;
编辑: 你必须构造它(粘贴上面的内容将使 expression
成为常量 0)
expression = sym('(1015424780204960143215273323910078528628754663952913658657288835138029704791686717561885487746105223164496264397062144000*(60345244216851610523130575942127473698515638085026410114070496399315754724985170479034760688327679099512793302302720*2^(1/2) - 85341062796188128379389251264141456937242003828816791711306294886439805159655912635773490018204138847163921224639041)*(x - 2)^2)/31504288346872372712061562812941419427167561153216213605146864117586323331155800294021400857537041202039565879856190821729945216935526707438840242294801134287960934949194864204307116208568439380511702602881');
调用 double(subs(expression, x, 1))
将此表达式计算为大约 -5.9492e+03
没有任何问题,这实际上是正确的值。但是,评估花费的时间太长,并且是我的应用程序中的一个巨大(或者说很小?)瓶颈。这就是为什么我打算将表达式转换为对双精度运算并且速度更快的匿名函数,如下所示:
evaluator = matlabFunction(expression, 'Vars', x);
结果是
@(x) (sqrt(2.0).*6.034524421685161e115-8.534106279618813e115) ...
.*(x-2.0).^2.*3.223131940086397e-86
我过去用这种方法取得了很大的成功。不幸的是,在这种情况下,evaluator(x)
对于 x 的任何值都计算为 0,因为第一行恰好为零。显然,这与双精度数的有效数字有限有关。
有没有办法解决这个问题?我可以告诉 MATLAB 考虑 x
的范围,以便它可以找到更好的常量表示吗?
问题可以通过 variable-precision arithmetic 解决:
>> expr = vpa(expression)
expr =
-5949.2156936801140978790460880789*(x - 2.0)^2
将其转换为函数
>> func =
function_handle with value:
@(x)(x-2.0).^2.*-5.949215693680114e3
然后评价一下
>> func(1)
ans =
-5.9492e+03
tl;dr: 我观察到在使用 matlabFunction
将符号表达式转换为函数句柄时完全失去了准确性。不知道有没有什么办法可以提高转换率,避免精度损失。
我有一个符号变量 x
,1 <= x && x <= 2
成立:
syms x real
assumeAlso(1 <= x & x <= 2);
我在该变量中有许多计算机生成的冗长符号表达式。其中一个看起来像这样:
expression = ( ...
1015424780204960143215273323910078528628754663952913658657288835138029704791686717561885487746105223164496264397062144000 ...
*( ...
60345244216851610523130575942127473698515638085026410114070496399315754724985170479034760688327679099512793302302720*2^(1/2) ...
- 85341062796188128379389251264141456937242003828816791711306294886439805159655912635773490018204138847163921224639041 ...
)*(x - 2)^2 ...
) ...
/31504288346872372712061562812941419427167561153216213605146864117586323331155800294021400857537041202039565879856190821729945216935526707438840242294801134287960934949194864204307116208568439380511702602881;
编辑: 你必须构造它(粘贴上面的内容将使 expression
成为常量 0)
expression = sym('(1015424780204960143215273323910078528628754663952913658657288835138029704791686717561885487746105223164496264397062144000*(60345244216851610523130575942127473698515638085026410114070496399315754724985170479034760688327679099512793302302720*2^(1/2) - 85341062796188128379389251264141456937242003828816791711306294886439805159655912635773490018204138847163921224639041)*(x - 2)^2)/31504288346872372712061562812941419427167561153216213605146864117586323331155800294021400857537041202039565879856190821729945216935526707438840242294801134287960934949194864204307116208568439380511702602881');
调用 double(subs(expression, x, 1))
将此表达式计算为大约 -5.9492e+03
没有任何问题,这实际上是正确的值。但是,评估花费的时间太长,并且是我的应用程序中的一个巨大(或者说很小?)瓶颈。这就是为什么我打算将表达式转换为对双精度运算并且速度更快的匿名函数,如下所示:
evaluator = matlabFunction(expression, 'Vars', x);
结果是
@(x) (sqrt(2.0).*6.034524421685161e115-8.534106279618813e115) ...
.*(x-2.0).^2.*3.223131940086397e-86
我过去用这种方法取得了很大的成功。不幸的是,在这种情况下,evaluator(x)
对于 x 的任何值都计算为 0,因为第一行恰好为零。显然,这与双精度数的有效数字有限有关。
有没有办法解决这个问题?我可以告诉 MATLAB 考虑 x
的范围,以便它可以找到更好的常量表示吗?
问题可以通过 variable-precision arithmetic 解决:
>> expr = vpa(expression)
expr =
-5949.2156936801140978790460880789*(x - 2.0)^2
将其转换为函数
>> func =
function_handle with value:
@(x)(x-2.0).^2.*-5.949215693680114e3
然后评价一下
>> func(1)
ans =
-5.9492e+03