1
$\begingroup$

I have a set of coupled differential equations, which I am trying to integrate parametrically. Rast is the parameter that is being varied. My code is as follows:

Constants

au = QuantityMagnitude@UnitConvert[Quantity[1, "AstronomicalUnit"], "Meters"]; 
c = QuantityMagnitude@UnitConvert[Quantity[1, "SpeedOfLight"], "MetersPerSecond"]; 
Qpr = 1; 
Lsun = QuantityMagnitude@UnitConvert[Quantity[1, "SolarLuminosity"], "Watts"]; 
Rsun = QuantityMagnitude@UnitConvert[Quantity[1, "SolarRadius"], "Meters"]; 
Msun = QuantityMagnitude@UnitConvert[Quantity[1, "SolarMass"], "Kilograms"]; 
G = QuantityMagnitude@UnitConvert[Quantity[1, "GravitationalConstant"], ("Meters"^2*"Newtons")/"Kilograms"^2]; 
year = QuantityMagnitudeUnitConvert[Quantity[1, "Years"], "Seconds"]; 
Myr = year*10^6; 
Gyr = year*10^9; 
Mwd = 0.6*Msun; 
Cst = 1.27; 
U = 1*10^17; 

Functions

L[t_] := (3.26*Lsun*(Mwd/(0.6*Msun)))/(0.1 + t/Myr)^1.18; 
Roche[dens_] := (0.65*Cst*Rsun*(Mwd/(0.6*Msun))^(1/3))/(dens/3000)^3^(-1); 
Papsis[t_] := a[t]*(1 - e[t]); 

Radiative Drag

RDdadtR = -((1/c^2)*((3*L[t]*Qpr*(2 + 3*e[t]^2))/(16*Pi*2000*Rast*a[t]*(1. - e[t]^2)^(3/2)))); 

RDdedtR = -((1/c^2)*((15*L[t]*e[t])/(32*Pi*Rast*2000*a[t]^2*Sqrt[1. - e[t]^2]))); 

RDsolR = ParametricNDSolveValue[{a'[t] == RDdadtR, e'[t] == RDdedtR, a[0] == 1*au, e[0] == 0.3, 
     WhenEvent[{Evaluate[Papsis[t] <= Roche[2000]]}, {tmax = t, "StopIntegration", Print["Target Accreted"]}]}, {a, e}, 
    {t, 0, 9*Gyr}, {Rast}]; 

RDasolR[Rast_, t_] := RDsolR[Rast][[1]][t]; 

RDesolR[Rast_, t_] := RDsolR[Rast][[2]][t]; 

Plot[Evaluate[Table[RDasolR[Rast, t], {Rast, 0.001, 0.01, 0.001}]], {t, 0, 9*Gyr}, FrameLabel -> {"Time", "Semi-major Axis"}]

Plot[Evaluate[Table[RDesolR[Rast, t], {Rast, 0.001, 0.01, 0.001}]], {t, 0, tmax}, FrameLabel -> {"Time", "Eccentricity"}]

An error is raised, namely Limiting value is not a machine sized real number. I don't know where I'm going wrong. Any help is greatly appreciated.

Improve this question
$\endgroup$
8
  • $\begingroup$ Could you upload Raw InputForm without In[]? $\endgroup$
    – Alex Trounev
    Commented Feb 2, 2021 at 21:06
  • $\begingroup$ Done! Hope that is more readable. $\endgroup$
    – testing09
    Commented Feb 2, 2021 at 21:32
  • 2
    $\begingroup$ These errors are occurring because (1) Gyr has the unit s attached to it, and tmax is not defined until after Plot has been called. In general, I recommend that you avoid using units in your computation. $\endgroup$
    – bbgodfrey
    Commented Feb 3, 2021 at 1:11
  • $\begingroup$ I understand the problem about Gyr having units, but I'm unsure why tmax is defined after Plot has been called, given it should be defined as the time the integration is stopped. Is there a way of getting around this? $\endgroup$
    – testing09
    Commented Feb 3, 2021 at 9:40
  • $\begingroup$ @testing09 there is a typo in your code in a[0] == 6*AU (you define au not AU). Also use non dimensional parameters by adding //QuantityMagnitude, for example year = UnitConvert[Quantity[1, "Years"], "Seconds"] //QuantityMagnitude. Even with this correction there are no events with a[t]*(1 - e[t]) <= Roche[2000] $\endgroup$
    – Alex Trounev
    Commented Feb 3, 2021 at 14:05

1 Answer 1

Reset to default
1
$\begingroup$

We can compute a, e first to show, that there are no such events we are looking for. This code is only working minimal example I been able to derive from the code above:

au = UnitConvert[Quantity[1, "AstronomicalUnit"], "Meters"] // 
   QuantityMagnitude;
c = UnitConvert[Quantity[1, "SpeedOfLight"], "MetersPerSecond"] // 
   QuantityMagnitude;
Qpr = 1;
Lsun = UnitConvert[Quantity[1, "SolarLuminosity"], "Watts"] // 
   QuantityMagnitude;
Rsun = UnitConvert[Quantity[1, "SolarRadius"], "Meters"] // 
   QuantityMagnitude;
Msun = UnitConvert[Quantity[1, "SolarMass"], "Kilograms"] // 
   QuantityMagnitude;
G = UnitConvert[
    Quantity[1, "GravitationalConstant"], ("Meters"^2*"Newtons")/
     "Kilograms"^2] // QuantityMagnitude;
year = UnitConvert[Quantity[1, "Years"], "Seconds"] // 
   QuantityMagnitude;
Myr = year*10^6;
Gyr = year*10^9;
Mwd = 0.6*Msun;
Cst = 1.27;
U = 1*10^17;


(*Functions*)

L[t_] := (3.26*Lsun*(Mwd/(0.6*Msun)))/(0.1 + t/Myr)^1.18;
Roche[dens_] := (0.65*Cst*
     Rsun*(Mwd/(0.6*Msun))^(1/3))/(dens/3000)^3^(-1);
Papsis[t_] := a[t]*(1 - e[t]);


(*Radiative Drag*)

RDdadtR = -((1/
       c^2)*((3*L[t]*Qpr*(2 + 3*e[t]^2))/(16*Pi*2000*Rast*
         a[t]*(1. - e[t]^2)^(3/2))));

RDdedtR = -((1/
       c^2)*((15*L[t]*e[t])/(32*Pi*Rast*2000*a[t]^2*
         Sqrt[1. - e[t]^2])));

RDsolR = ParametricNDSolveValue[{Derivative[1][a][t] == RDdadtR, 
    Derivative[1][e][t] == RDdedtR, a[0] == 6*au, e[0] == 0.3}, {a[t],
     e[t]}, {t, 0, 9*Gyr}, {Rast}];


RDasolR[Rast_] := RDsolR[Rast][[1]];

RDesolR[Rast_] := RDsolR[Rast][[2]];


{Plot[Evaluate[
   Table[RDsolR[Rast][[1]], {Rast, 0.001, 0.01, 0.001}]], {t, 0, 
   9*Gyr}, FrameLabel -> {"Time", "Semi-major Axis"}, Frame -> True], 
 Plot[Evaluate[
   Table[RDsolR[Rast][[2]], {Rast, 0.001, 0.01, 0.001}]], {t, 0, 
   9*Gyr}, FrameLabel -> {"Time", "Eccentricity"}, Frame -> True], 
 Plot[Evaluate[
   Table[RDsolR[Rast][[1]] (1 - RDsolR[Rast][[2]]) - 
     Roche[2000], {Rast, 0.001, 0.01, 0.001}]], {t, 0, 9*Gyr}, 
  FrameLabel -> {"Time", "Eccentricity"}, AxesOrigin -> {0, 0}]}

Figure 1 The last picture is the expression Papsis[t] - Roche[2000]] used for event detection. There is no time where Papsis[t] - Roche[2000]]<=0, and therefore there are no such events.

Update 1. In a case of initial data a[0] == 1*au, e[0] == 0.3 there are several events and we can detect and plot solution with events as follows (here we use a[t] in astronomical unit, and t in Myr)

au = UnitConvert[Quantity[1, "AstronomicalUnit"], "Meters"] // 
    QuantityMagnitude // Rationalize;
c = UnitConvert[Quantity[1, "SpeedOfLight"], "MetersPerSecond"] // 
    QuantityMagnitude // Rationalize;
Qpr = 1;
Lsun = UnitConvert[Quantity[1, "SolarLuminosity"], "Watts"] // 
    QuantityMagnitude // Rationalize;
Rsun = UnitConvert[Quantity[1, "SolarRadius"], "Meters"] // 
    QuantityMagnitude // Rationalize;
Msun = UnitConvert[Quantity[1, "SolarMass"], "Kilograms"] // 
    QuantityMagnitude // Rationalize;
G = UnitConvert[
     Quantity[1, "GravitationalConstant"], ("Meters"^2*"Newtons")/
      "Kilograms"^2] // QuantityMagnitude // Rationalize;
year = UnitConvert[Quantity[1, "Years"], "Seconds"] // 
    QuantityMagnitude // Rationalize;
Myr = year*10^6;
Gyr = year*10^9;
Mwd = 6/10*Msun;
Cst = 127/100;
U = 10^17;


(*Functions*)

L[t_] := (326/100*Lsun*(Mwd/(6/10*Msun)))/(1/10 + t)^(118/100);
Roche[dens_] := (65/100*Cst*
     Rsun*(Mwd/(6/10*Msun))^(1/3))/(dens/3000)^(-1/3);
Papsis[t_] := au a[t]*(1 - e[t]);

(*Radiative Drag*)

RDdadtR = -((1/
       c^2)*((3*L[t]*Qpr*(2 + 3*e[t]^2))/(16*Pi*2000*Rast*
         au a[t]*(1. - e[t]^2)^(3/2))));

RDdedtR = -((1/
       c^2)*((15*L[t]*e[t])/(32*Pi*Rast*2000*au^2 a[t]^2*
         Sqrt[1. - e[t]^2])));


Do[tmax[r] = 9 Gyr/Myr; 
 sol[r] = NDSolve[{Derivative[1][a][t] == (Myr/au) RDdadtR /. {Rast ->
          r}, Derivative[1][e][t] == Myr RDdedtR /. {Rast -> r}, 
      a[0] == 1, e[0] == 3/10, 
      WhenEvent[{Evaluate[Papsis[t] <= Roche[2000]]}, {tmax[r] = t, 
        "StopIntegration"}]}, {a[t], e[t]}, {t, 0, 9 Gyr/Myr}, 
     AccuracyGoal -> 20, PrecisionGoal -> 20, 
     WorkingPrecision -> 35][[1]] // Quiet; 
 tm[r] = If[tmax[r] < 9 Gyr/Myr, tmax[r], 9 Gyr/Myr]; 
 p[r] = {Plot[a[t] /. sol[r], {t, 0, tm[r]}, 
    FrameLabel -> {"Time, Myr", "Semi-major Axis, AU"}, Frame -> True,
     PlotLabel -> Row[{"Rast = ", r 1.}], PlotRange -> All], 
   Plot[e[t] /. sol[r], {t, 0, tm[r]}, 
    FrameLabel -> {"Time, Myr", "Eccentricity"}, Frame -> True, 
    PlotLabel -> Row[{"tmax = ", N[tm[r] 1., 3]}], 
    PlotRange -> All]};, {r, 10^-3, 1/100, 10^-3}]
Table[p[r], {r, 10^-3, 1/100, 10^-3}]

Figure 2

Improve this answer
$\endgroup$
7
  • $\begingroup$ @AlexTrounev- thank you so much! One last question- you've removed the explicit time dependence from the expression. How do I now scale the x axis so it's in Myr? $\endgroup$
    – testing09
    Commented Feb 6, 2021 at 14:19
  • 1
    $\begingroup$ @testing09 Actually in your code Time is in second, not in Myr. But if you wish to plot Time, Gyr or Time, Myr with FrameLabel option we can scaled Plot output. $\endgroup$
    – Alex Trounev
    Commented Feb 6, 2021 at 15:43
  • 2
    $\begingroup$ @testing09 We can also use scaling Plot[Evaluate[ Table[RDsolR[Rast][[1]] /. t -> ts Myr, {Rast, 0.001, 0.01, 0.001}]], {ts, 0, 9*Gyr/Myr}, FrameLabel -> {"Time, Myr", "Semi-major Axis"}, Frame -> True] $\endgroup$
    – Alex Trounev
    Commented Feb 6, 2021 at 16:55
  • $\begingroup$ I have a question about your solution- I understand that you have specified the stop condition in the Plot function. Does Mathematica automatically know that it should stop plotting when this expression goes to zero? The reason I ask is that I am trying to extract the time at which the integration is stopped. Sorry for so many questions but I'm trying to get better at Mathematica :) – $\endgroup$
    – testing09
    Commented Feb 6, 2021 at 17:03
  • 1
    $\begingroup$ Oh I understand- my original problem was trying to get the integration to stop when that condition was met. This way the plots don't break. If the condition is met, plots 1 and 2 don't work. e.g. for a value of a[0]=3 au, the first 2 plots break. Hence, why I was trying to use WhenEvent to stop the integration from proceeding once the condition is met $\endgroup$
    – testing09
    Commented Feb 6, 2021 at 17:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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