Vector ODE solution is not periodic/ as expected
1 view (last 30 days)
Show older comments
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);
Thanks in advance!
0 Comments
Accepted Answer
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
James Tursa
on 10 Mar 2022
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
More Answers (1)
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
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 ?
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!