How many significant digits should I use for simulating aircraft equations of motion?

I have coded the equations of motion in Matlab, and I am using the ode45 function to simulate the nonlinear dynamics of F/A-18 Hornet without any controller (just the airframe). I have given initial conditions as trim conditions.

d2r = pi/180;           % degree to rad

V       =  438.653328;  % Airspeed , ft/s
beta    =  0*d2r;       % Sideslip Angle, rad
alpha   =  10*d2r;      % Angle-of-attack, rad

p       =  0*d2r;       % Roll rate, rad/s
q       =  0*d2r;       % Pitch rate, rad/s
r       =  0*d2r;       % Yaw rate, rad/s

phi     =  0*d2r;       % Roll Angle, rad
theta   =  10*d2r;      % Pitch Angle, rad
psi     =  0*d2r;       % Yaw Angle, rad

pN = 0;                 % X position in Earth Frame, ft
pE = 0;                 % Y position in Earth Frame, ft
h = -25000;             % Z position in Earth Frame, ft

During the first iteration, C_m has a very positive small value (1e-10), which leads to some pitching moment when multiplied by S, c, and qbar (even though it's small). But, as time keeps on increasing, this value of C_m goes on increasing and as a result, a large moment is created. Due to the high precision of Matlab, the pitching moment coefficient is not for the given trim condition.

My question is: should I round off the values of C_m to 3 or 4 digits? By rounding off, the C_m will be zero and I get M as zero.

The below figure shows how velocity is varying with time when the aforementioned initial conditions are given.

This figure is obtained after rounding off C_m to 4 significant digits. This figure is obtained without rounding off C_m to 4 significant digits. Alpha and theta values without rounding off C_m. • Wait, why does C_m (pitching moment coefficient of the whole aircraft) increase? What is the value (or at least sign) of its main component, the $C_m^{\alpha}$ (or $^{n_z}$) derivative? It should determine pitch stability, given that the rest is basically damping. If it's pitch stable (the derivative is negative), the computational errors should not be significant.
– Zeus
Mar 13 '20 at 0:18
• The sign of the pitching moment is positive, which shows that it is unstable. So, without external control, it will not be in equilibrium (irrespective of trim condition input), right? Mar 13 '20 at 5:12
• It looks to me as if the simulation rate is too low. The rounding of the numbers could create better results just because it happens to better align with the simulation rate. (I think it's called aliasing effect)
– Jan
Mar 13 '20 at 6:32
• Positive pitching moment simply means nose up. What determines stability is its derivatie (by angle of attack or normal acceleration), which is a parameter of the model itself. (It is, of course, a function of other parameters such as CG location). If the derivative is positive, the model/aircraft is unstable, and it doesn't matter much which kind of disturbance you induce. If you cheat with C_m, you will see instability the moment you start changing other things (e.g. normal fligtht control).
– Zeus
Mar 13 '20 at 6:33
• In fact, judging by the last graph you added (I understand it's airspeed - but pitch or alpha would be more indicative), the model is stable. It's just the initial trim condition is a bit off, and it accumulates over time, which results in a typical phugoid motin. By reducing precision, you are forcing it to be 'on' trim. Basically, it looks like there is nothing wrong with the model as it is. If you add control, you should be able to fly it.
– Zeus
Mar 13 '20 at 6:41

From your description, it's not an issue with the aerodynamic coefficients. If as you have described previously that you're modeling F-18, then the airframe should be unstable longitudinally. A 1e-10 deviation (assuming it applies to body rates, body speeds, etc. as well) is well-within typical trim error bound; in fact, a deviation is unavoidable when dealing with computational floating point. But if it's longitudinally unstable, any deviation will magnify over time exponentially.

• As far as I know, the 'original' F/A-18 (A-D) is statically stable. E-F is about neutral or marginally unstable. I wouldn't expect its model depart the trim condition significantly (with no external disturbances) in less than 7-10 seconds...
– Zeus
Mar 12 '20 at 23:54
• @JZYL, As you said the value of 1e-10 is well within the trim error, should I round off? Because Matlab is not considering it zero and it grows with time. Mar 13 '20 at 5:17
• @Zeus, As you have stated, velocity is constant for the first 45 sec. Mar 13 '20 at 5:19

Computers in these days do not care about "significant decimal digits" for computation as they use binary IEEE-754 representation internally. The typical number sizes are float (most often 32 bits in these days) and double (most often 64 bit). C offers 80-bit extended precision (long double) and sometimes 128 bit double is available. You can read about sizes here. In addition, some programming languages and packages can do arbitrary precision but it is usually slow.

There are often "corner cases" in model that may require unusually high precision. For instance, the expression (a-b)/(c-d) may be very tough to compute when a gets close to b while c gets close to d. Because of such cases, I usually prefer to use the highest precision available (least long double) in simulations. It may even be corner cases that require complex type arithmetic, even if both input and result are of the real type (sqrt(-1)**2).

Smaller sizes are often not even faster. They sometimes make sense if you work with data matrices so large that memory becomes a problem, of if you use a gamer GPU card for computing that does not feature double precision.

For human readable reports, it is common to leave enough decimal digits to represent the value with known accuracy, or one more, but not longer. For instance, if the value is known with the 0.1 accuracy, it could be printed as 0.6 or 0.66 for instance, but not as 0.6666666666666666. But if this text output is not for humans and will be read back by computer, there is no reason to limit the number of digits. Also, if you know more digits, rounding can only do the same or worse.

• This is all hardly relevant for the question, but I can't help but nitpick that 80-bit precision is a feature of the Intel FPU (and some others), rather than that of C. In fact, C doesn't guarantee much in this regard; specifically, in Visual C long double is 64 bit, whereas in gcc it's 80 on the same hardware.
– Zeus
Mar 16 '20 at 7:22