1
$\begingroup$

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

Radiative Drag

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

RDsolR\[Rho]a = ParametricNDSolveValue[{Derivative[1][a][t] == RDdadtR\[Rho]a, 
     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.

$\endgroup$
2
  • $\begingroup$ What are the desired ranges of the three parameters and of the plot? $\endgroup$
    – Bob Hanlon
    Feb 6, 2021 at 14:54
  • $\begingroup$ @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 $\endgroup$
    – testing09
    Feb 6, 2021 at 16:03

1 Answer 1

1
$\begingroup$
$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"]];
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 = 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]);

Radiative Drag

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

enter image description here

$\endgroup$
3
  • $\begingroup$ @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. $\endgroup$
    – testing09
    Feb 6, 2021 at 19:01
  • 1
    $\begingroup$ @testing09 - probably; what change do you want? $\endgroup$
    – Bob Hanlon
    Feb 6, 2021 at 19:04
  • $\begingroup$ I want the y axis of the a plot to be in units of au , instead of metres. $\endgroup$
    – testing09
    Feb 6, 2021 at 19:06

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.