Unexpected numerical error in built-in cross product

I encountered an issue when using the built in cross product function in the Command Window.

rTAN = [6.284251782248913e+06; 2.176644923429249e+06; -2.615238551283982e+06]
vTAN = [-5.225826867852853e+02; 6.171595435460577e+03; 4.262061413416244e+03]
hT = 5.371945921706661e+10
hTN = cross(rTAN, vTAN)

hT is the specific angular momentum of the orbit, and hTN is its vector representation in the inertial reference frame. The magnitude of hTN should be identical to hT within machine precision, but, when I checked it, I got

>> norm(hTN) - hT
ans =

  -2.1230e+03
>> ans/hT
ans =

  -3.9521e-08

Can anyone explain to me why I'm getting this fairly large relative error and how I could go about fixing it? I'm using version R2018a.

[the cyclist]

How did you calculate that value of hT?

[Otis Hudnut]

hT was calculated using the formula hT = sqrt(aT*mu*(1-eT^2)), where aT is the semimajor axis of the transfer orbit and eT is the eccentricity of the transfer orbit.

aT = 7.243582592195826e+06
eT = 0.022905923178296
mu = 3.9860044e+14

aT and eT were calculated from an analytical solution to the node-to-node non-Hohmann transfer equations.

[Matt J]

@Otis,

You haven't shown us anything to prove that cross() is at fault. Maybe your hT was calculated incorrectly.

[Otis Hudnut]

@Matt J

The values of rTAN and vTAN, which both agreed with magnitudes for those values calculated using the same aT, eT, and mu as given in my previous comment within machine precision, were calculated using the same input values I used to calculate hT. rTA, vTA, and hT were all calculated with classical Keplerian orbital relations, which are correct for the two-body problem which I am trying to solve.

thetaNT = 3.048417059938535e+02 % degrees
eT = 0.022905923178296
aT = 7.243582592195826e+06
mu = 3.9860044e+14

rTA = aT*(1-eT^2)/(1+eT*cosd(thetaNT))
vTRA = (mu/hT*eT*sind(thetaNT)
vTTanA = mu/hT*(1+eT*cosd(thetaNT))
vTA = sqrt(vTRA^2 + vTTanA^2)
hT = sqrt(aT*mu*(1-eT^2))

The conversion of rTA, vTRA, vTTanA into rTAN and vTAN was achieved with verified unit vectors and direction cosine matrices, all of which were confirmed to have a norm of 1.

Therefore, I cannot see another place where the error could have come from besides cross(). If you can suggest one, please do, however.

Answer 1

[James Tursa]

When I calculate things from scratch, everything works to the expected precision. E.g., run this code:

% From post
rTAN = [6.284251782248913e+06; 2.176644923429249e+06; -2.615238551283982e+06]; % (m)
vTAN = [-5.225826867852853e+02; 6.171595435460577e+03; 4.262061413416244e+03]; % (m/s)
hT = 5.371945921706661e+10; % (m^2/s)
aT = 7.243582592195826e+06; % (m)
eT = 0.022905923178296;
mu = 3.9860044e+14; % m^3/s^2

% Raw calculations
r = norm(rTAN);
v = norm(vTAN);
hvec = cross(rTAN, vTAN);
h = norm(hvec);
E = (v^2)/2 - mu/r;
a = -mu/(2*E);
evec = ( (v^2 - mu/r)*rTAN - dot(rTAN,vTAN)*vTAN ) / mu;
e = norm(evec);
h_from_ae = sqrt(a*mu*(1-e^2));

% Look at results of h
disp(' ');
disp('h raw calculations results:');
fprintf('h from norm(cross(r,v))   = %20.8f\n',h);
fprintf('h from sqrt(a*mu*(1-e^2)) = %20.8f\n',h_from_ae);
fprintf('Difference                = %20.8f\n',h-h_from_ae);
fprintf('Relative Difference       = %20.8g\n',(h-h_from_ae)/h);

% Look at a vs post
disp(' ');
disp('a raw calculations vs post:');
fprintf('a from calculations       = %20.8f\n',a);
fprintf('a from post               = %20.8f\n',aT);
fprintf('Difference                = %20.8f\n',a-aT);
fprintf('Relative Difference       = %20.8g\n',(a-aT)/a);

% Look at e vs post
disp(' ');
disp('e raw calculations vs post:');
fprintf('e from calculations       = %20.8f\n',e);
fprintf('e from post               = %20.8f\n',eT);
fprintf('Difference                = %20.8f\n',e-eT);
fprintf('Relative Difference       = %20.8g\n',(e-eT)/e);

% Look at h vs post
disp(' ');
disp('h raw calculations vs post:');
fprintf('h from calculations       = %20.8f\n',h);
fprintf('h from post               = %20.8f\n',hT);
fprintf('Difference                = %20.8f\n',h-hT);
fprintf('Relative Difference       = %20.8g\n',(h-hT)/h);

To get this result:

h raw calculations results:
h from norm(cross(r,v))   = 53719457094.01964600
h from sqrt(a*mu*(1-e^2)) = 53719457094.01966100
Difference                =          -0.00001526
Relative Difference       =       -2.8404585e-16
 
a raw calculations vs post:
a from calculations       =     7243582.59219583
a from post               =     7243582.59219583
Difference                =           0.00000000
Relative Difference       =        3.8571628e-16
 
e raw calculations vs post:
e from calculations       =           0.02290765
e from post               =           0.02290592
Difference                =           0.00000172
Relative Difference       =        7.5275803e-05
 
h raw calculations vs post:
h from calculations       = 53719457094.01964600
h from post               = 53719459217.06661200
Difference                =       -2123.04696655
Relative Difference       =       -3.9521006e-08

It looks like we differ in the eccentricity value. We only agree to less than 5 decimal digits, and this is propagating into your h calculation, so that is the source of the differences you are seeing. I would advise looking into how you are calculating eccentricity. When I use norm(eccentricity_vector) to calculate e I get consistent results as shown above, with relative differences in the e-16 range which is what you would expect from double precision.

P.S. I would strongly advise that you put physical units in all of your code! It really, really helps when you (or someone else) has to look at your code in the future.

[Otis Hudnut]

I see where you're going with this, but since vTAN was calculated from the posted hT (from aT, eT) and eT values, I'm unsure about where the disparity in eccentricity values is appearing from. Could you comment on that? Is it possible that the issue is appearing because of a loss of significant figures in mu?

[James Tursa]

Bottom line is that your eccentricity and your angular momentum magnitude don't match your pos & vel vectors to the precision you are expecting. The semi-major axis does match. So it is really up to you to examine the methods you are using to derive the position and velocity from the orbital elements, because they don't seem to work to the precision you expect. Since you didn't post any of that code I can't comment on it. For all I know the problem could involve other orbital elements as well that you didn't post (e.g., true anomaly, inclination, etc.)

Answer 2

[Matt J]

To build trust in cross(), I compute norm(hTN) below from first principles, with only elementary operations:

x = [6.284251782248913e+06; 2.176644923429249e+06; -2.615238551283982e+06]; %rTAN
y = [-5.225826867852853e+02; 6.171595435460577e+03; 4.262061413416244e+03]; %vTAN


Nx=sqrt(sum(x.^2)); %norm of x
Ny=sqrt(sum(y.^2)); %norm of y

theta=acosd(x.'*y/Nx/Ny);


>> normhTN=Nx*Ny*sind(theta) %Only elementary functions

normhTN =

     5.371945709401965e+10


>> norm(cross(x,y))  %Using cross()

ans =

     5.371945709401965e+10