matlabFunction produces file with insufficient precision?

I am using the symbolic toolbox to generate equations of motion for a dynamics simulation, then I use matlabFunction to create a file for the resulting function. I then use the equations of motion (and their gradients, computed using automatic differentiation) in an optimal control problem. My problem is set up correctly and is well scaled, but the solver is running into numerical difficulties. Insufficient precision of the equations of motion files occurred to me as a possible cause. For example, here is an excerpt from the code generated by matlabFunction:

function EOM = eom_flight(t,in2,tau)
%EOM_FLIGHT
%    EOM = EOM_FLIGHT(T,IN2,TAU)

%    This function was generated by the Symbolic Math Toolbox version 5.6.
%    16-Aug-2011 21:24:50

...<more code>...

    t30 = t25+t26+t27+t28+t29-t32-t33-1.266084114070474e72;
    t31 = 1.0./t30;
    t34 = t15.*t16.*tau.*4.259680785351676e71;
    t35 = t17.*t18.*tau.*4.259680785351676e71;
    EOM = [dx;t31.*(t15.*t17.*3.089190909270927e52+t16.*t18.*3.089190909270927e52-t15.*t23.*5.363399619039539e67-...

...<more code and end of file>...

The fixed precision of the printed numbers and the wildly large exponents seem suspicious. Any suggestions? For instance, is there any way to increase the precision of the numbers printed in the file in attempt to avoid round-off error? Also, it would be nice to normalize the values used in the calculations, as I doubt they actually need to span tens/hundreds of orders of magnitude....

Answer 1

[Kai Gehrs]

Hi,

without knowing all the details and understanding the overall problem, I would guess it could help to set up the data strict symbolically. When you define your equations, you could try to use exact numbers (fractions like e.g. 1/10 instead of floats 0.1, pi instead of 3.1416, sqrt(2) instead of 1.414 etc.).

Just as an example: when you set up equations in the form

syms x y;
eq1 = sqrt(2)*x^2 + 1/10*log(y)

then you can compute the derivatives and work on the equations by doing arithmetic and the data will remain exact (in terms of a SYM object). This can help to avoid numerical instabilities (cancellation etc.) when doing arithmetics.

Finally, when calling matlabFunction and switching to the numerical ODE solver you will of course end up with floats.

But please note that I am just guessing. It may be worth a try.

Best regards,

-- Kai

Answer 2

[Walter Roberson]

Unfortunately, symbolic mathematics does tend to generate these large-exponent multipliers when you have floating point values in the formula. The floating point numbers tend to get converted to rational numbers, which usually involves multiplying by 2^53 and then factoring out the gcd from the top and bottom of the ratio. Square such a value and you are working with numbers on the order of 2^100... and so on.

Use integers or rationals when possible; e.g., 0.37 -> sym('37/100'). And as Kai suggests, use the symbolic form of constants. Though sqrt(sym('2')) or sym('sqrt(2)') should be used rather than sqrt(2) as sqrt(2) by itself will not be a symbolic constant.

[Matt]

Thanks for your responses, Kai, Walter. Now that I re-defined my constants as symbolics, the largest number printed by matlabFunction is on the order 10^22 as opposed to ~10^80. While that's better, is there anything else I can do, such as increase the number of digits printed?