Search code examples
matlabsymbolic-math

How to avoid loss of accuracy when converting symbolic expressions to functions?


tl;dr: I observe a complete loss of accuracy when converting a symbolic expression to a function handle with matlabFunction. I wonder if there are ways to improve the conversion in order to avoid the loss of accuracy.


I have a symbolic variable x for which 1 <= x && x <= 2 holds true:

syms x real
assumeAlso(1 <= x & x <= 2);

I have many computer-generated, lengthy symbolic expression in that variable. One of them looks like this:

expression = ( ...
  1015424780204960143215273323910078528628754663952913658657288835138029704791686717561885487746105223164496264397062144000 ...
  *( ...
     60345244216851610523130575942127473698515638085026410114070496399315754724985170479034760688327679099512793302302720*2^(1/2) ...
   - 85341062796188128379389251264141456937242003828816791711306294886439805159655912635773490018204138847163921224639041 ...
   )*(x - 2)^2 ...
) ...
/31504288346872372712061562812941419427167561153216213605146864117586323331155800294021400857537041202039565879856190821729945216935526707438840242294801134287960934949194864204307116208568439380511702602881;

Edit: You have to construct it with (pasting the above will make expression a constant 0)

expression = sym('(1015424780204960143215273323910078528628754663952913658657288835138029704791686717561885487746105223164496264397062144000*(60345244216851610523130575942127473698515638085026410114070496399315754724985170479034760688327679099512793302302720*2^(1/2) - 85341062796188128379389251264141456937242003828816791711306294886439805159655912635773490018204138847163921224639041)*(x - 2)^2)/31504288346872372712061562812941419427167561153216213605146864117586323331155800294021400857537041202039565879856190821729945216935526707438840242294801134287960934949194864204307116208568439380511702602881');

Calling double(subs(expression, x, 1)) evaluates this expression to approximately -5.9492e+03 without any problems, which actually is the correct value. However, the evaluation takes far too long and is a huge (or rather tiny?) bottleneck in my application. That’s why I intend to convert the expression to an anonymous function that operates on doubles and is much faster, like so:

evaluator = matlabFunction(expression, 'Vars', x);

The result being

@(x) (sqrt(2.0).*6.034524421685161e115-8.534106279618813e115) ...
.*(x-2.0).^2.*3.223131940086397e-86

I had much success with this approach in the past. Unfortunately in this case, evaluator(x) evaluates to 0 for any value of x, because the first row is exactly zero. Obviously, this has to do with the limited number of significant digits for doubles.

Are there ways to work around this? Can I tell MATLAB to consider the range of x so that it can find a better representation of the constants?


Solution

  • The problem can be fixed with variable-precision arithmetic:

    >> expr = vpa(expression)
    expr =
    -5949.2156936801140978790460880789*(x - 2.0)^2
    

    Convert it to a function

    >> func =
      function_handle with value:
        @(x)(x-2.0).^2.*-5.949215693680114e3
    

    Then evaluate it

    >> func(1)
    ans =
      -5.9492e+03