1
$\begingroup$

I am trying to solve a set of coupled differential equations, where I would create plots for a(t) and e(t), where the parameter RAST can be varied. However, every time I try to plot, I get blank plots. How would I go about doing this?

To emphasise, I am looking to generate 2 separate plots; one for a(t) and one for e(t), which each contain multiple, different coloured graphs showing the variation of the parameter, RAST.

My code:

In[3]:=
LSUN = 3.83*10^26; 
MSUN = 1.989*10^30; 
MWD = MSUN; 
year = QuantityMagnitude[UnitConvert[
     Quantity[1, "Years"], "Seconds"]]; 
km = QuantityMagnitude[UnitConvert[
     Quantity[1, "Kilometers"], 
     "Meters"]]; timelimit = 
   year*2*10^9; Myr = 1*10^6*year; 

In[10]:=
PROPCONST = 1.27; 

In[7]:=
c = QuantityMagnitude[UnitConvert[
     Quantity[1, "SpeedOfLight"], 
     "Meters/second"]]; \[Rho] = 2000; 
MAST[rast_] := (4/3)*Pi*rast^3*\[Rho]; 
RSUN = 696340000; 
AU = QuantityMagnitude[UnitConvert[
     Quantity[1, "AstronomicalUnit"], 
     "meters"]]; 

In[9]:=
L[t_] := (3.26*LSUN*(MWD/(0.6*MSUN)))/
    (0.1 + t/Myr)^1.18; 

In[10]:=
dadt = (-(1/c^2))*((RAST^2*L[t]*
      (2 + 3*e[t]^2))/(4*MAST[RAST]*a[t]*
      (1 - e[t]^2)^(3/2))); 
dedt = (-(5/c^2))*((RAST^2*L[t]*e[t])/
     (8*MAST[RAST]*a[t]^2*
      Sqrt[1 - e[t]^2])); 

In[16]:=
prsol = ParametricNDSolveValue[
   {Derivative[1][a][t] == dadt, 
    Derivative[1][e][t] == dedt, 
    a[0] == 2*AU, e[0] == 0.5}, {a, e}, 
   {t, 0, timelimit}, {RAST}]

Out[16]= ParametricFunction[ <> ]

In[18]:=
asol[RAST_, t_] := prsol[RAST][[1]][t]; 
esol[RAST_, t_] := prsol[RAST][[2]][t]; 
```
$\endgroup$
2
  • 1
    $\begingroup$ tried Plot[Evaluate[Table[prsol[ r][[1]][t], {r, 1, 10}]], {t, 0, timelimit}] and Plot[Evaluate[Table[prsol[ r][[2]][t], {r, 1, 10}]], {t, 0, timelimit}]? $\endgroup$
    – kglr
    Nov 30, 2020 at 16:29
  • $\begingroup$ Brilliant. I think that works. Thanks! $\endgroup$
    – testing09
    Nov 30, 2020 at 17:54

1 Answer 1

1
$\begingroup$
Row[{Plot[Evaluate[Table[prsol[r][[1]][t], {r, 1, 10}]], {t, 0, timelimit}, 
    ImageSize -> Medium], 
  Plot[Evaluate[Table[prsol[r][[2]][t], {r, 1, 10}]], {t, 0, timelimit}, 
    ImageSize -> Medium]}]

enter image description here

Using {t, 0, timelimit / 10^4} gives

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.