# Solving multiple differential equations parametrically

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

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

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

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}, {ChartingFindTicks[{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.

• @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. Mar 9, 2021 at 2:35
• Another approach would be to combine the two instances of ParametricNDSolve. Mar 9, 2021 at 2:51
• @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. Mar 9, 2021 at 11:22
• @bbgodfrey I have tried your suggested simple solution, however this doesn't work. How would I combine the two instances into one ParametricNDSolveValue? Mar 9, 2021 at 11:43
• 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. Mar 9, 2021 at 12:14

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


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


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


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

• Can you post the full code? I still get the same error for some reason :( Mar 9, 2021 at 16:49
• @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? Mar 9, 2021 at 16:58
• 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! Mar 9, 2021 at 18:39
• 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 Mar 9, 2021 at 18:53
• @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. Mar 9, 2021 at 19:12