Matlab Maximum Valid Number

Hello, I encountered the following issues while solving numerical issues with Matlab.

Explained by the code below, there is a difference between the results obtained by substituting the same expression in parentheses and the results obtained by generally calculating the same expression in parentheses.

syms dl x1 y1 dx1 dy1 x2 y2 dx2 dy2 E A s G M rho

dl = 100;
E = 131*10^9;                 
A = 2*10^-6;                   

pel =[6599999.99999543;10;1.00000000000139;0;6600099.99999543;10;1.00000000000139;0];

N111 = double(subs((3*(dl^2*dx1^2 + dl^2*dx2^2 + dl^2*dy1^2 + dl^2*dy2^2 - 14*dl^2 + 6*dl*dx1*x1 ...
                - 6*dl*dx1*x2 + 6*dl*dx2*x1 - 6*dl*dx2*x2 + 6*dl*dy1*y1 - 6*dl*dy1*y2 + 6*dl*dy2*y1 - 6*dl*dy2*y2 ...
                + 24*x1^2 - 48*x1*x2 + 24*x2^2 + 24*y1^2 - 48*y1*y2 + 24*y2^2))/(35*dl^4), ...
                [x1, y1, dx1, dy1, x2, y2, dx2, dy2], pel'));

x1 = pel(1,1);
y1 = pel(2,1);
dx1 = pel(3,1);
dy1 = pel(4,1);
x2 = pel(5,1);
y2 = pel(6,1);
dx2 = pel(7,1);
dy2 = pel(8,1);

N11 =             (3*(dl^2*dx1^2 + dl^2*dx2^2 + dl^2*dy1^2 + dl^2*dy2^2 - 14*dl^2 + 6*dl*dx1*x1 ...
                - 6*dl*dx1*x2 + 6*dl*dx2*x1 - 6*dl*dx2*x2 + 6*dl*dy1*y1 - 6*dl*dy1*y2 + 6*dl*dy2*y1 - 6*dl*dy2*y2 ...
                + 24*x1^2 - 48*x1*x2 + 24*x2^2 + 24*y1^2 - 48*y1*y2 + 24*y2^2))/(35*dl^4);

When comparing N111 and N11, you might think that both result values are very small, and this difference is negligible, but this difference seems to work by multiplying them by very large values in the future.

(In fact, the above equation is the result of integration. It was closer to the result of N111 when we performed the numerical integration.)

Looking at forums, it is expected to be related to the maximum effective number of MATLAB.

I wonder what causes these problems, which one is closer to the correct answer.

[Aquatris]

Not an expert in this but I think the issue is that evaluating the symbolic expression does not use double precision until the final step. Where as N11 stores and evaluates everything in double precision from the get go, where the chance of floating point error occuring and accumulating is high for values that are close to the highest precision matlab can handle.

[Yohan Ahn]

Thank you for your reply, I am currently using the [m] unit scale, but if I change it to the [km] unit scale, will it be an alternative to solving the problem? (e.g. switching from 6600000m to 6600km)

[Aquatris]

No, changing the magnitude will not help you. It is not about the magnitude but the precision. Your numbers require too high of a precision to be operated on digitally. You can read more about it in the answer here.

Some examples of effects of this precision can be seen by simply adding a small number to 1, and subtracting 1 again to see that result is not always 0 (notice the more precision we use, the more errors we get):

a = 1+1e-2-1;
a

a = 0.0100

b = 1+1e-7-1;
b

b = 1.0000e-07

c = 1+1e-15-1;
c

c = 1.1102e-15

d = 1+1e-16-1;
d

d = 0