Vector ODE solution is not periodic/ as expected

3 views (last 30 days)
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!

Accepted Answer

James Tursa
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.
  3 Comments
Lucas Parkins
Lucas Parkins on 10 Mar 2022
Yes! This worked and was easily altered to fit myh problem, thank you!

Sign in to comment.

More Answers (1)

Torsten
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
  2 Comments
Torsten
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);
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 ?

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!