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

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?

Answer 1

[James Tursa]

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]';

[Elza Zandere]

Thank you James, it works now!

[BHOOMIKA]

Hello ! I have the same question and I want to solve it using 0de45, could you help me with the code when we use that function.