1
$\begingroup$

I have the following code, where I am trying to solve one set of differential equations parametrically. I then want to use this solution to solve another equation parametrically. However, I get the error: ParametricNDSolveValue::pdvar: Dependent variables {Rast,s} cannot depend on parameters {W,Rast,[Rho],s0,a0}.

The code is as follows:

Clear["Global`*"]

SetOptions[Plot, Frame -> True, Axes -> False, LabelStyle -> {FontFamily -> "Baskerville", FontSize -> 12, Bold}, FrameStyle -> Directive[Thick, Black], 
   ImageSize -> Medium]; 

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

scrit[dens_] := (7.48*Sqrt[dens/2000])/10^4; 

Yarkovsky

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]*Qpr)/(c^2*(32*Pi*Rast*\[Rho]*a[t]^2*Sqrt[1 - e[t]^2]))); 

YdadtR\[Rho]a = -((1/c)*((3*L[t])/(16*Pi*Rast*\[Rho]*(1 - e[t]^2)*Sqrt[G*Mwd*a[t]]))); 

YdedtR\[Rho]a = -((1/c)*((3*L[t]*(1 - Sqrt[1 - e[t]^2]))/(32*Pi*Rast*\[Rho]*e[t]*Sqrt[G*Mwd*a[t]^3]))); 

dadtR\[Rho]a = Evaluate[YdadtR\[Rho]a*UnitStep[Rast - 1] + RDdadtR\[Rho]a]; 

dedtR\[Rho]a = Evaluate[YdedtR\[Rho]a*UnitStep[Rast - 1] + RDdedtR\[Rho]a]; 

YORPdsdtWR\[Rho]s = (W/(2*Pi*\[Rho]*Rast^2))*(1/(a[t]^2*Sqrt[1 - e[t]^2]))*(U*(L[t]/Lsun)); 

solR\[Rho]a = ParametricNDSolveValue[{Derivative[1][s][t] == YORPdsdtWR\[Rho]s, Derivative[1][a][t] == dadtR\[Rho]a, Derivative[1][e][t] == dedtR\[Rho]a, s[0] == s0, a[0] == a0, 
     e[0] == 0.3}, {s, a, e}, {t, 0, 9*Gyr}, {W, Rast, \[Rho], s0, a0}]; 

fticks = {{Automatic, Automatic}, {Charting`FindTicks[{0, 1}, {0, 1/Myr}], Automatic}}; 

Manipulate[Column[{Style["YORP working plot", Bold], Plot[fun[func, t], {t, 0, 9*Gyr}, FrameTicks -> fticks, 
     Epilog -> {Red, Dashed, InfiniteLine[{{0, scrit[\[Rho]]}, {10, scrit[\[Rho]]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}], Style["Compiled Plot", Bold], 
    If[comp === {}, Plot[fun[func, t], {t, 0, 9*Gyr}, FrameTicks -> fticks, Epilog -> {Red, Dashed, InfiniteLine[{{0, scrit[\[Rho]]}, {10, scrit[\[Rho]]}}]}, 
      PlotStyle -> {Directive[Blue, Thickness[0.01]]}], Plot[comp, {t, 0, 9*Gyr}, FrameTicks -> fticks, 
      Epilog -> {Red, Dashed, InfiniteLine[{{0, scrit[\[Rho]]}, {10, scrit[\[Rho]]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}]]}], {{func, 1}, {1 -> "s"}}, 
  {{W, 0}, 0, 1, 0.1, Appearance -> "Labeled"}, {{Rast, 50}, 1, 4000, 10, Appearance -> "Labeled"}, {{\[Rho], 2000}, 1000, 7000, 50, Appearance -> "Labeled"}, 
  {{s0, 0, "s0 (\!\(\*SuperscriptBox[\(rads\), \(-1\)]\))"}, 0, 100, 1, Appearance -> "Labeled"}, {{a0, 1, "a0 (au)"}, 1, 20, 1, Appearance -> "Labeled"}, 
  Button["Append", AppendTo[comp, fun[func, t]]], Button["Reset", comp = {}], TrackedSymbols :> {func, W, Rast, \[Rho], s0, a0}, 
  Initialization :> {comp = {}, fun[func_, t_] := solR\[Rho]a[W, Rast, \[Rho], s0, a0][[1]][t]}]

Any help would be appreciated.

$\endgroup$
5
  • $\begingroup$ @Bill I just saw the same problem myself. A simple solution is to replace {a, e} by {a} in the first ParametricNDSolve and remove [[1]] from the definition of asol. $\endgroup$
    – bbgodfrey
    Mar 9 at 2:35
  • $\begingroup$ Another approach would be to combine the two instances of ParametricNDSolve. $\endgroup$
    – bbgodfrey
    Mar 9 at 2:51
  • $\begingroup$ @Bill Rast is a parameter which I am varying. Once I solve the second equations, I can use Manipulate to vary Rast and see the behaviour of the system. $\endgroup$
    – testing09
    Mar 9 at 11:22
  • $\begingroup$ @bbgodfrey I have tried your suggested simple solution, however this doesn't work. How would I combine the two instances into one ParametricNDSolveValue? $\endgroup$
    – testing09
    Mar 9 at 11:43
  • $\begingroup$ I've removed asol with a slight refactor. I have attempted to try to resolve why Rast has a time dependence given it is just a parameter which I vary and thus has no time-dependence. $\endgroup$
    – testing09
    Mar 9 at 12:14
1
$\begingroup$

The code has three difficulties. First, asol appears in YORPdsdtWRρs as asol[Rast, ρ, a0, t], which evaluates to

Clear[asol]
asol[Rast_, ρ_, a0_, t_] := solRρa[Rast, ρ, a0][[1]][t]
asol[Rast, ρ, a0, t]
(* Rast[t] *)

To prevent this, instead use

Clear[asol]
asol[Rast_?NumericQ, ρ_?NumericQ, a0_?NumericQ, t_] := solRρa[Rast, ρ, a0][[1]][t]
asol[Rast, ρ, a0, t]
(* asol[Rast, ρ, a0, t] *)

which does not evaluate prematurely. Next, plot asol with the default parameter values given in Manipulate in the queston.

Plot[Evaluate[asol[50, 2000, 1, t]], {t, 0, 1.726 10^-6}, 
    ImageSize -> Large, AxesLabel -> {t, a}, LabelStyle -> {15, Bold, Black}]

enter image description here

Visibly, the stellar radius goes to zero at about t = 1.7260 10^-6, which is the source of the error message from Evaluate[asol[50, 2000, 1, t]].

Finally, ParametricNDSolve objects with attempting to evaluate solWRρs, because "Dependent variables {s,asol[Rast,ρ,a0,t]} cannot depend on parameters {W,Rast,ρ,s0,a0}". To solve this final problem, merge the computations of solRρa and solWRρs.

YORPdsdtWRρs = (W/(2*Pi*ρ*Rast^2))*(1/(a[t]^2*Sqrt[1 - ecc^2]))*(U*(L[t]/Lsun));
solWRρs = ParametricNDSolveValue[{Derivative[1][s][t] == YORPdsdtWRρs, s[0] == s0, 
    Derivative[1][a][t] == dadtRρa, Derivative[1][e][t] == dedtRρa, a[0] == a0, 
    e[0] == 0.3}, {s[t], a[t], e[t]}, {t, 0, 9*Gyr}, {W, Rast, ρ, s0, a0}];

(Note that I have not corrected the second problem by changing the upper bound on t immediately above, because it causes no substantial harm.)

Plot[Evaluate[First@solWRρs[1, 50, 2000, 0, 1]], {t, 0, 1.726 10^-6}, 
    ImageSize -> Large, AxesLabel -> {t, s}, LabelStyle -> {15, Bold, Black}]

enter image description here

The plot of a, and of e as well, can be recovered by

Plot[Evaluate[Rest@solWRρs[1, 50, 2000, 0, 1]], {t, 0, 1.726 10^-6}, 
    ImageSize -> Large, AxesLabel -> {t, a}, LabelStyle -> {15, Bold, Black}]
$\endgroup$
7
  • $\begingroup$ Can you post the full code? I still get the same error for some reason :( $\endgroup$
    – testing09
    Mar 9 at 16:49
  • $\begingroup$ @testing09 The whole code is everything in the question up to and including dedtRρa = Evaluate[YdedtRρa*UnitStep[Rast - 1] + RDdedtRρa]; plus the code for YORPdsdtWRρs and solWRρs in my answer. It does not include asol. Does this help? $\endgroup$
    – bbgodfrey
    Mar 9 at 16:58
  • $\begingroup$ the default parameter you have used for t is wrong. It should be 9 Gyr. Also, now I get a different error: `"Input value {5.79812*10^12} lies outside the range of data in the \ interpolating function. Extrapolation will be used" However, the graph you have in your answer has the correct shape, but I cannot reproduce it. Do you have access to the full code? Btw really appreciate your help with this! $\endgroup$
    – testing09
    Mar 9 at 18:39
  • $\begingroup$ I have edited the code in the question to reflect the above potential answer. But I now get an interpolation error. Could you check this over plz $\endgroup$
    – testing09
    Mar 9 at 18:53
  • $\begingroup$ @testing09 Do not use 9 Gyr for a plot range, because a goes to zero at about t = 1.7260 10^-6, and ParametricNDSolve stops there. Trying to evaluate the InterpolationFunction for much larger t gives an error. $\endgroup$
    – bbgodfrey
    Mar 9 at 19:12

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.