Unexpected numerical error in built-in cross product
4 views (last 30 days)
Show older comments
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.
4 Comments
Matt J
on 28 Jun 2019
@Otis,
You haven't shown us anything to prove that cross() is at fault. Maybe your hT was calculated incorrectly.
Answers (2)
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.
2 Comments
James Tursa
on 1 Jul 2019
Edited: James Tursa
on 1 Jul 2019
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.)
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
0 Comments
See Also
Categories
Find more on Get Started with MATLAB in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!