Matlab ode45 Help: Solving Orbital Propagation Issues

[Hashmeer]

Homework Statement

I need to create an orbit propagator for the orbit of a specified body and graph it. I've written a code that works and correctly gives me a circular orbit for a satellite in geosynchronous orbit and it appears to correctly display the orbit of the Earth for one year. I'm running into trouble when I attempt to propagate further into the future. My graph then begins to spiral inwards instead of maintaining its path.

The Attempt at a Solution

I am using ode45 for this. Here is the code (first is the ode function, second is the run script):

function dP = orbitProp( t, P )
% Equations of motion for a 2D orbit
% P = [x position, y position, x velocity (vx), y velocity (vy)] 
 
mu = 1.32712440018e+11; % gravitational parameter in km^3/s^2
 
r = sqrt( P(1)^2 + P(2)^2);
 
dP = zeros(4,1);
dP(1) = P(3);
dP(2) = P(4);
dP(3) = -(mu/r^3)*P(1);
dP(4) = -(mu/r^3)*P(2);
end% set initial conditions, (position in km and velocity in km/s)
x0 = [ -1.458240050323911e+08, 3.114907245044354e+07, -6.715255132769428, -29.26108860586488];
 
% integrate forward in time 1 day (1 year = 3.15581e+07)
[T, Y] = ode45(@orbitProp,[0,3.15581e+07],x0);
% plot results
plot(Y(:,1),Y(:,2));
xlabel('X(km)');
ylabel('Y(km)');

The data comes from JPL's HORIZONS and are the values as of 03/09/2011 at 12:00am.

My instructor said that using odeoptions should allow me to clear this up, but I'm not sure exactly how to use them or what I should even be looking for to clean the graph up.

Any help/hints are appreciated. Thanks!

 

[sozkan]

Well I guess it is quite late but I am working on a similar project and you just needed to adjust tolerances using odeset. Using the following code gives a circle even for a year of integration.

options = odeset('RelTol',1e-10,'AbsTol',1e-11);
[T, Y] = ode45(@orbitProp,[0,365 * 3.15581e+07],x0,options);