# Problem with ParametricNDSolveValue & WhenEvent

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"];
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]);

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.

• Could you upload Raw InputForm without In[]? Feb 2, 2021 at 21:06
• Done! Hope that is more readable. Feb 2, 2021 at 21:32
• 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. Feb 3, 2021 at 1:11
• 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? Feb 3, 2021 at 9:40
• @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] Feb 3, 2021 at 14:05

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]);

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])));

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}]}


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]);

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}]


• @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? Feb 6, 2021 at 14:19
• @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. Feb 6, 2021 at 15:43
• @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] Feb 6, 2021 at 16:55
• 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 :) – Feb 6, 2021 at 17:03
• 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 Feb 6, 2021 at 17:25