Generate a Runge-Kutta Integrator solver for the orbital motion equation, when given the initial conditions. Execute for 500 seconds.

The Problem statement is as followed.

Using the constant for km/s. and the inital conditions

km

and

km/s

The orbital motion equation is:

integrate using runge kutta RK4 for an h =10 seconds and an end time of 500 seconds

So I am working on this problem and I am stuck on the state variable transformation to make it a first order system and subsequently applying the Runge Kutta integrator. I am not sure how to work through it seeing that we are dealing with vectors instead of a scalar . I did the Ode45 approach as well as an attempt as well but i am not certain if i am correct. Any pointers would be greatly appreciated.

function xdot = exfunc(t,x)
mu = 398600.64;
xdot(1) = x(2);
xdot(2) = (-mu/((x(1))^(3)))*x(1);
xdot(3) = x(4);
xdot(4) = (-mu/((x(3))^(3)))*x(3);
xdot(5) = x(6);
xdot(6) = (-mu/((x(5))^(3)))*x(5);
xdot = xdot';
end

clear
clc

[t,x] = ode45('exfunc',[0 1],[5371.844151186472, 3.230645742388105,...
    -4141.062031065303, 3.522344029484922, 460.1400917227622,...
    -5.981911152962826])
plot(t,x(:,[1,3,5]))


R_end = x(end,[1,3,5])
V_end = x(end,[2,4,6])

[James Tursa]

FYI, you are missing a negative sign in your rdotdot equation above.

[Von Clarke]

yeah thanks, i noticed some inconsistancies with my question a little too late

Answer 1

[James Tursa]

First, I find it easier to keep the position elements together, and the velocity elements together, in the state vector. So I would make this definition:

x(1) = x position

x(2) = y position

x(3) = z position

x(4) = x velocity

x(5) = y velocity

x(6) = z velocity

Then in your derivative function you can work with the vectors more easily. Be sure to use the full position vector for calculating r. And I would advise appending comments with units on all lines with constants. E.g.,

function xdot = exfunc(t,x)
mu = 398600.64; % km^3/s^2
r = x(1:3);
rdot = x(4:6);
rmag = norm(r);
rdotdot = -(mu/rmag^3)*r;
xdot = [rdot;rdotdot];
end

And of course that means you would need to rearrange your initial conditions that you feed to ode45:

r0 = [5371.844151186472; -4141.062031065303; 460.1400917227622]; % km
v0 = [3.230645742388105; 3.522344029484922; -5.981911152962826]; % km/s
t0 = 0; % s
tend = 500; % s
[t,x] = ode45(@exfunc,[t0 tend],[r0;v0]);
plot(t,x(:,1:3))
R_end = x(end,1:3)
V_end = x(end,4:6)