Vector ODE solution is not periodic/ as expected

Hi,

I am coding a solver for the orbital differential equations, which of course i expect to be periodic as the orbit is nearly circular. Instead The resulting orbit makes little physical sense and all parameters either diverge or tend to unrealistic constants.

The 6 elements of the ODE represent (x,y,z) position and velocity components' derivatives.

This is the differential equation function:

function dxdt = gravity(x, DU)
    dxdt = zeros(6,1);
    dxdt(1) = x(4);  
    dxdt(2) = x(5);
    dxdt(3) = x(6);

    dxdt(4) = -DU^2* x(1)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3; 
    dxdt(5) = -DU^2* x(2)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3; 
    dxdt(4) = -DU^2* x(3)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3;
return
end

To set up the solver and solve the equations:

r = [2408.8; -6442.7;-0000.0];
v = [4.2908; 1.6052; 6.0795]; %rDot
x = [r;v]; 
DU = norm(r);

y0 = [r;v];
span = [0 20*pi]; %Represents 10 periods in non-dimensional coords
f = @(t,y) gravity(x, DU);
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[ts, ys] = ode45(f, span, y0, opts);

Thanks in advance!

Answer 1

[James Tursa]

This index 4

dxdt(4) = -DU^2* x(3)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3;

needs to be index 6:

dxdt(6) = -DU^2* x(3)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3;

So you would have something like:

r = norm(x(1:3));
dxdt(1:3) = x(4:6);
dxdt(4:6) = -G*(m1+m2) * x(1:3) / r^3;

As long as your DU^2 gives the same value as G*(m1+m2) in dimensionless coords then you would be OK. I haven't checked into that part of it, nor have I checked to see if your initial conditions match what would be expected for a circular orbit in a dimensionless system. Frankly, just coding things up to use units might be simpler than the dimensionless system you seem to be setting up.

[Lucas Parkins]

Thanks for the advice, the x => y swap and indexing changed my solutions but still not as expected, so i suspect my equations are at fault.

The dimensionless solution improves the accuracy of the solution significantly, at least in theory, however i have run the code as a dimensional prblem without any luck.

[James Tursa]

An example using units:

% Simple example of circular orbit at 10,000 m radius
mu = 3.98574405096E14; % Earth (m^3/s^2)
r = 10000; % (m)
v = sqrt(mu/r); % circular orbit velocity
y0 = [r;0;0;0;v;0]; % Initial circular orbit
n = sqrt(mu/r^3); % mean motion
period = 2*pi/n; % orbit period
tspan = [0 period];
f = @(t,y) gravity(t,y,mu);
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[ts, ys] = ode45(f, tspan, y0, opts);
plot(ys(:,1),ys(:,2));
grid on
axis square

function dydt = gravity(t,y,mu)
R = y(1:3);
dydt = [y(4:6);-mu*R/norm(R)^3];
end

[Lucas Parkins]

Yes! This worked and was easily altered to fit myh problem, thank you!

Answer 2

[Torsten]

r = [2408.8; -6442.7;-0000.0];
v = [4.2908; 1.6052; 6.0795]; %rDot
x = [r;v]; 
DU = norm(r);
y0 = [r;v];
span = [0 20*pi]; %Represents 10 periods in non-dimensional coords
f = @(t,y) gravity(y,DU);
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[ts, ys] = ode45(f, span, y0, opts);

function dxdt = gravity(x,DU)
dxdt = zeros(6,1); 
dxdt(1) = x(4);  
dxdt(2) = x(5);
dxdt(3) = x(6);
dxdt(4) = -DU^2* x(1)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3; 
dxdt(5) = -DU^2* x(2)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3; 
dxdt(4) = -DU^2* x(3)/sqrt(x(1)^2+x(2)^2+x(3)^2)^3;
end

[Lucas Parkins]

This changes nothing... i see that you have swapped an x for a y when declaring the function but this does not alter the results at all.

[Torsten]

If you write

f = @(t,y) gravity(x,DU);

you will solve your system with the initial DU and x-vector inserted in the ODEs, thus

dxdt(1) = v(1);  
dxdt(2) = v(2);
dxdt(3) = v(3);
dxdt(4) = -DU^2* r(1)/sqrt(r(1)^2+r(2)^2+r(3)^2)^3; 
dxdt(5) = -DU^2* r(2)/sqrt(r(1)^2+r(2)^2+r(3)^2)^3; 
dxdt(4) = -DU^2* r(3)/sqrt(r(1)^2+r(2)^2+r(3)^2)^3;

The solution should change substantially if you use

f = @(t,y) gravity(y,DU);

instead.

But if the results are not as expected, you should check your equations.

Maybe DU also has to be updated during computation as

DU = norm(x(1:3))

instead of using

DU = norm(r )

for all times t ?