5 views (last 30 days)

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.

James Tursa
on 28 Jun 2019

Edited: 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.

James Tursa
on 1 Jul 2019

Sign in to comment.

Matt J
on 28 Jun 2019

Edited: 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.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## Direct link to this comment

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

⋮## Direct link to this comment

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

## Direct link to this comment

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

⋮## Direct link to this comment

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

## Direct link to this comment

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

⋮## Direct link to this comment

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

## Direct link to this comment

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

⋮## 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.