Vector ODE solution is not periodic/ as expected

3 views (last 30 days)
Lucas Parkins on 9 Mar 2022
Commented: Lucas Parkins on 10 Mar 2022
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);

James Tursa on 9 Mar 2022
Edited: James Tursa on 10 Mar 2022
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 on 10 Mar 2022
Yes! This worked and was easily altered to fit myh problem, thank you!

Torsten on 9 Mar 2022
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
Torsten on 9 Mar 2022
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);
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))
DU = norm(r )
for all times t ?

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!