# ParametricNDSolve Value and Manipulate with 3 parameters

I am trying to solve a set of coupled differential equations, using ParametricNDSolveValue, where there are 3 parameters. The parameters are Rast, \[Rho] and one of the initial condition values, a[0]. I want to plot the result of this integration as a Manipulate plot, where each parameter can be independently 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 = QuantityMagnitude[UnitConvert[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\[Rho]a = -((3*L[t]*Qpr*(2 + 3*e[t]^2))/(c^2*(16*Pi*\[Rho]*Rast*a[t]*(1 - e[t]^2)^(3/2))));
RDdedtR\[Rho]a = -((15*L[t]*e[t])/(c^2*(32*Pi*Rast*\[Rho]*a[t]^2*Sqrt[1 - e[t]^2])));

Derivative[1][e][t] == RDdedtR\[Rho]a, a[0] == a0, e[0] == 0.3}, {a[t], e[t]}, {t, 0, 9*Gyr},
{Rast, \[Rho]}];


The parameters are Rast, \[Rho], and the initial condition value, a0. I am unsure how to create a Manipulate plot with all three parameters.

PARAMETER DOMAINS: Rast from 0.001 to 0.01 \[rho] from 1000 to 7000 a0 from 3*au to 20*au and the Plot domain: t=0 to t=9 Gyr

Any help would be greatly appreciated.

• What are the desired ranges of the three parameters and of the plot? Commented Feb 6, 2021 at 14:54
• @BobHanlon parameter domains: Rast from 0.001 to 0.01, \[rho] from 1000 to 7000, a0 from 3*au to 20*au and the Plot domain: t=0 to t=9 Gyr Commented Feb 6, 2021 at 16:03

\$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global*"]


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 = QuantityMagnitude[UnitConvert[Quantity[1, "Years"], "Seconds"]];
Myr = year*10^6;
Gyr = year*10^9;
Mwd = 6/10*Msun;
Cst = 1.27;
U = 10^17;


Functions

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


RDdadtRρa = -((3*L[t]*Qpr*(2 + 3*e[t]^2))/
(c^2*(16*Pi*ρ*Rast*a[t]*(1 - e[t]^2)^(3/2))));
RDdedtRρa = -((15*L[t]*e[t])/
(c^2*(32*Pi*Rast*ρ*a[t]^2*Sqrt[1 - e[t]^2])));


Include a0 as a parameter

RDsolRρa =
ParametricNDSolveValue[{a'[t] == RDdadtRρa, e'[t] == RDdedtRρa,
a[0] == a0, e[0] == 3/10}, {a, e},
{t, 0, 9*Gyr}, {Rast, ρ, a0}];


EDIT: Plot of a rescaled to units of au

Manipulate[
factor = If[func == 1, 1/au, 1];
Plot[factor*RDsolRρa[Rast, ρ, a0*au][[func]][t],
{t, 0, 9 Gyr}],
{{func, 1}, {1 -> "a", 2 -> "e"}},
{{Rast, 0.005}, 0.001, 0.01, 0.001, Appearance -> "Labeled"},
{{ρ, 3000}, 1000, 7000, 50, Appearance -> "Labeled"},
{{a0, 10, "a0 (au)"}, 3, 20, 0.2, Appearance -> "Labeled"}]


• @BonHanlon Amazing! Is there a way of changing just the y axis scale just for a? I've changed the x -axis but find that the changes reflect for both the a and e plots. Commented Feb 6, 2021 at 19:01
• @testing09 - probably; what change do you want? Commented Feb 6, 2021 at 19:04
• I want the y axis of the a plot to be in units of au` , instead of metres. Commented Feb 6, 2021 at 19:06