Asked by Otis Hudnut
on 28 Jun 2019

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.

Answer by James Tursa
on 28 Jun 2019

Edited by James Tursa
on 28 Jun 2019

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
on 30 Jun 2019

James Tursa
on 1 Jul 2019

Sign in to comment.

Answer by Matt J
on 28 Jun 2019

Edited by Matt J
on 28 Jun 2019

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## the cyclist (view profile)

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/469429-unexpected-numerical-error-in-built-in-cross-product#comment_719535

## Otis Hudnut (view profile)

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/469429-unexpected-numerical-error-in-built-in-cross-product#comment_719537

## Matt J (view profile)

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/469429-unexpected-numerical-error-in-built-in-cross-product#comment_719542

## Otis Hudnut (view profile)

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/469429-unexpected-numerical-error-in-built-in-cross-product#comment_719547

Sign in to comment.