Dustinsfl
- 2,217
- 5
How can I extract time data from a system 2nd order ODEs in Mathematica?
This discussion focuses on extracting time data from a second-order ordinary differential equation (ODE) using Mathematica. The provided code simulates the 2-body problem, specifically modeling the trajectory of a spacecraft influenced by Earth's gravity. Key functions utilized include NDSolve for solving the ODEs and ParametricPlot3D for visualizing the trajectory. Additional suggestions include using Table to generate numerical results alongside graphical outputs for better data representation.
PREREQUISITESResearchers, physicists, and engineers involved in computational physics, particularly those working with orbital simulations and numerical analysis in Mathematica.
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]