MATLAB: Numerical integration of the differential equation of motion of the two body problem

integration

I have to write a code that integrates the differential equation of motion of the 2-body problem numerically, starting from initial values of position and velocity in the three-dimensional space, using this equation:
Initial values of a Geostationary satellite.
I have created a model, integration and main code file.
Model: model.m
function xdot = model(x,u);
global mu r
xdot(1,1) = x(2);
xdot(2,1) = (-mu/r^3)*x(1)
Integration: rk4int.m
function x = rk4int(model, h, x, u)
k1 = h*feval(model, x, u);
k2 = h*feval(model, x+k1/2, u);
k3 = h*feval(model, x+k2/2, u);
k4 = h*feval(model, x+k3, u);
x = x + (k1 + 2*k2 + 2*k3 + k4)/6;
Main code:
Define Constants and initial conditions
global mu r v
mu = 398600; % km^3/s^2
R0 = [42164 0 0]
V0 = [3.0746 0 0]
r = norm(R0)
v = norm(V0)
t = 0
E = v*v/2 – mu/r
a = -mu/(2*E)
n = sqrt(mu/(a*a*a))
T = 2*pi/n %period
Define the parameters
stepsize = T/1000 % Integration step size
comminterval = 0.01 % Communications interval
EndTime = T % Duration of the simulation (final time)
i = 0 % Initialise counter for data storage
Initial states
u = 0
x = [42164,0,0]'
xdot = [3.0746,0,0]'
Time
for time = 0:stepsize:EndTime
if rem(time,comminterval)==0
i = i+1 % increment counter
tout(i) = time % store time
xout(i,:) = x % store states
xdout(i,:) = xdot % store state derivatives
end
xdot = model(x,u)
x = rk4int('model',stepsize,x,u)
end
However, when I try to run the main code, it stops at the ''x = rk4int('model',stepsize,x,u)'' line gives me the following error:
What is the issue? How can I fix it?

Best Answer

Your biggest problem is that you don't carry enough states in your derivative function. Your ODE is a 3D 2nd order equation, so that means you will need 3*2 = 6 states to carry. But your derivative function only calculates two scalar derivatives when it should be calculating six scalar derivatives. Also, you need to calculate r from the current state vector ... not pass in a constant for this. So, first let's decide what goes into your six element state vector x:
x(1) = x position
x(2) = y position
x(3) = z position
x(4) = x velocity (i.e., xdot)
x(5) = y velocity (i.e., ydot)
x(6) = z velocity (i.e., zdot)
Given those definitions, your derivative function should look something like this instead:
function xdot = model(x,u);
global mu
xdot = zeroes(size(x));
xdot(1:3) = x(4:6);
r = norm(x(1:3));
xdot(4:6) = (-mu/r^3)*x(1:3);
end
Then for initial state you would simply have
R0 = [42164 0 0];
V0 = [0 3.0746  0];  % <-- changed!  The velocity needs to be in the ydot element for circular orbit
x = [R0,V0]';
Related Question