Dustinsfl
- 2,217
- 5
How can I extract time data from a system 2nd order ODEs in Mathematica?
The discussion revolves around extracting time data from a system of second-order ordinary differential equations (ODEs) in Mathematica, specifically in the context of simulating the 2-body problem. Participants explore methods for visualizing and retrieving numerical data from the simulation.
There is no explicit consensus on the best method to extract time data, as participants propose different approaches and modifications without resolving which is superior.
The discussion includes various assumptions about the simulation parameters and the desired output format, which may affect the applicability of the proposed solutions.
Ackbach said:Can you post the ODE's and your Mathematica code so far?
Numerical Simulation of the 2-Body Problem
ClearAll["Global`*"];
We begin with of the the necessary data for this problem...
M = 5974*10^21; (* mass of Earth, kg *)
m = 1000; (* mass of spacecraft , kg *)
\[Mu] = 3.986*10^5; (* gravitaional parameter, based on km units of length, \
km/s for velocity *)
Rearth = 6378; (* radius of the Earth, km *)Simulation Inputs
r0 = {3950.55, 43197.9, 0};(* initial position vector, km *)
v0 = {3.3809, -7.25046, 0}; (* initial velocity vector, km *)
Days = 1/10; (* elapsed time of simulation days *)
\[CapitalDelta]t = Days*24*3600;(* convert elapsed days to seconds *)
s = NDSolve[
{
x1''[t] == -(\[Mu]/(Sqrt[x1[t]^2 + x2[t]^2 + x3[t]^2])^3)*x1[t],
x2''[t] == -(\[Mu]/(Sqrt[x1[t]^2 + x2[t]^2 + x3[t]^2])^3)*x2[t],
x3''[t] == -(\[Mu]/(Sqrt[x1[t]^2 + x2[t]^2 + x3[t]^2])^3)*x3[t],
x1[0] == r0[[1]], (* intial x-position of satellite *)
x2[0] == r0[[2]],(* intial y-position of satellite *)
x3[0] == r0[[3]],(* intial y-position of satellite *)
x1'[0] == v0[[1]],(* intial vx-rel of satellite *)
x2'[0] == v0[[2]],(* intial vy-rel of satellite *)
x3'[0] == v0[[3]](* intial vy-rel of satellite *)
},
{x1, x2, x3},
{t, 0, \[CapitalDelta]t}
];
Plot of the Trajectory Relative to Earth
g1 = ParametricPlot3D[
Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, \[CapitalDelta]t},
PlotStyle -> {Red, Thick}];
g2 = Graphics3D[{Blue, Opacity[0.6], Sphere[{0, 0, 0}, Rearth]}];Show[g2, g1, Boxed -> False]