1
$\begingroup$

Im modelling mercury, venus,earth and mars orbiting the sun for a 2 body problem.

Below are the initial position velocity conditions for the planets in question. Distance units are in AU Velocity units are in AU/day

Initial conditions

T = {88.0, 224.7, 365.2, 687.0};(*period in days*)
Ecc = {0.20563069, 0.00677323, 0.01671022, 0.09341233};(*eccentricity of orbit*)   
\[Mu] = 2.9591220823*10^-4;(*standard gravitational parameter of sun*)
M = {174.79439000000002`, 0.`, -2.48284000000001`, 
19.412480000000016`};(*Mean Anomoly*)
\[Psi] = M + 2*Ecc*Sin[M Degree] + 5/4*Ecc^2*Sin[2*M Degree] + 
13/12*Ecc^3*Sin[3*M Degree];(*True Anomoly*)
a = {0.38709893, 0.7233319899999999, 1.00000011, 1.52366231};
r = Table[a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc [[i]] Cos[\[Psi] [[i]] Degree ]), {i, 4}] (*orbital position*)
rx = Table[r[[i]] Cos[\[Psi][[i]] Degree], {i, 4}](*x component of position*)
ry = Table[r[[i]] Sin[\[Psi][[i]] Degree ], {i, 4}](*y component of position*)
v = Table[Sqrt[\[Mu] (2/r[[i]] - 1/a[[i]])], {i, 4}](*orbital velocity*)
vx = v Cos[\[Psi] Degree](*x component of velocity*)
vy = v Sin[\[Psi] Degree](* y component of velocity*)

Solving using NDSolve

Here i am using NDsolve to solve the ODE's for the equations of motion in x,y direction

planet = Table[NDSolve[{(x''[t] + (\[Mu]*x[t])/(x[t]^2 + y[t]^2)^(
   3/2) == 0, (y''[t] + (\[Mu]*
    y[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, x[0] == rx[[i]], 
 y[0] == ry[[i]], x'[0] == vx[[i]], 
 y'[0] == vy[[i]]}, {x, y}, {t,800, 0.1}], {i,4}];

The problem

If i were to set the intial conditons to the magnitudes of postions and velocity, the simulation works fine, however, its when i use the components of (x,y) for position and velocity that i encounter the error message.

NDSolve::ndsz: At t == 76.8007476284912`, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 204.5293849935109`, step size is effectively zero; singularity or stiff system suspected.

NDSolve::ndsz: At t == 333.0376141002715`, step size is effectively zero; singularity or stiff system suspected.

General::stop: Further output of NDSolve::ndsz will be suppressed during this calculation.

What can i do to fix this?

Thank you for reading

$\endgroup$
1
$\begingroup$

It is necessary to correctly calculate the orbital velocity

T = {88.0, 224.7, 365.2, 687.0};(*period in days*)Ecc = {0.20563069, 
  0.00677323, 0.01671022, 
  0.09341233};(*eccentricity of orbit*)\[Mu] = \
0.00029591220828349645438389263285760599`10.122911634152604;\
(*standard gravitational parameter of sun*)M = {174.79439000000002`, 
  0.`, -2.48284000000001`, 
  19.412480000000016`};(*Mean Anomoly*)\[Psi] = 
 M + 2*Ecc*Sin[M Degree] + 5/4*Ecc^2*Sin[2*M Degree] + 
  13/12*Ecc^3*Sin[3*M Degree];(*True Anomoly*)a = {0.38709893, 
  0.7233319899999999, 1.00000011, 1.52366231};
r = Table[
   a[[i]] (1 - Ecc[[i]]^2)/(1 + Ecc[[i]] Cos[\[Psi][[i]] Degree]), {i,
     4}] (*orbital position*);
rx = Table[
   r[[i]] Cos[\[Psi][[i]] Degree], {i, 4}](*x component of position*);
ry = Table[
   r[[i]] Sin[\[Psi][[i]] Degree], {i, 4}](*y component of position*);
v = Table[
   Sqrt[\[Mu] (2/r[[i]] - 1/a[[i]])], {i, 4}](*orbital velocity*);
vx = -v Sin[\[Psi] Degree](*x component of velocity*);
vy = v Cos[\[Psi] Degree](*y component of velocity*);


pl = {Mercury, Venus, Earth, Mars};

planet = Table[
   NDSolve[{x''[t] + (\[Mu]*x[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, 
     y''[t] + (\[Mu]*y[t])/(x[t]^2 + y[t]^2)^(3/2) == 0, 
     x[0] == rx[[i]], y[0] == ry[[i]], x'[0] == vx[[i]], 
     y'[0] == vy[[i]]}, {x, y}, {t, 0, 800}], {i, 4}];

{Table[ParametricPlot[{x[t], y[t]} /. planet[[i]], {t, 0, 800}, 
   PlotLabel -> pl[[i]]], {i, 1, Length[planet]}], 
 ParametricPlot[
  Evaluate[Table[{x[t], y[t]} /. planet[[i]], {i, Length[pl]}]], {t, 
   0, 800}, PlotLegends -> pl]}

Figure 1

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.