0
$\begingroup$

I am still learning the intricacies of Mathematica. I have an NDSolve within which I am trying to vary a specific parameter (I understand that ParameticNDSolve is likely a more efficient tool--the results it is producing are physically incorrect, I am hoping NDSolve will allow for the correct solution). To do this, I have a Table on the outside of the NDSolve. The result is a table with interpolating functions. I need help figuring out how to extract these results from the table so that I can graph them. Thank you for the help!

Note: Once I can extract values from this table I plan to build two more tables so as to vary additional variables (d and zH below).


P[r_, z_] := (2. r ((2 - (4 r)/((1 + r)^2 + 
             z^2)) EllipticK[(4 r)/((1 + r)^2 + z^2)] - 
       2 EllipticE[(4 r)/((1 + r)^2 + z^2)]))/(\[Pi] ((4 r)/((1 + 
           r)^2 + z^2)) \[Sqrt]((1 + r)^2 + z^2));
Zmaxreg = 15;
Zminreg = -Zmaxreg;
Rmaxreg = 54;
Rminreg = 0;

d = 0; 
zH = 0; 
pref = 1;
a = 4;

Clear[H2, H]
Quiet[
 Table[sol = NDSolve[
    {(*EQUATIONS *)
     D[H2[r, z], z] == -(zH2 + d)*H2[r, z]*(P[r, z]/pref )^a,
     D[H[r, z], z] == 
      2*d*H2[r, z]*(P[r, z]/pref )^a - zH*H[r, z]*(P[r, z]/pref )^a,
     
     (*INITIAL CONDITIONS*)
     H2[r, Zminreg] == 1,
     H[r, Zminreg] == 0
     },
    (*VARIABLES TO SOLVE FOR*)
    {H2, H},
    (*RANGE TO SOLVE OVER*)
    {r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg}
    ],
  {zH2, 1, 10^3, 10}
  ]
 ]

H2[r_, z_] = H2[r, z] /. sol[[1]];
H[r_, z_] = H[r, z] /. sol[[1]];

~~~~~
$\endgroup$
4
  • 1
    $\begingroup$ Please post a minimal working example. $\endgroup$
    – cvgmt
    Jun 30, 2022 at 0:38
  • $\begingroup$ Well, you already know there's a Part ([[]]), don't you? That's exactly what you need. Please check the document of it more carefully. $\endgroup$
    – xzczd
    Jun 30, 2022 at 2:50
  • $\begingroup$ @xzczd I read through the documentation prior to posting and wasn't able to find the answer I needed. Are you able to elaborate on how exactly I would be able to pull out a single interpolating function from the table? $\endgroup$
    – Bill
    Jul 7, 2022 at 22:53
  • $\begingroup$ It's confusing to me you're having difficulty with [[]]. (Especially after seeing you claimed you know how to do that in MATLAB in comment below. It's similar in Mathematica. ) How do you understand the sol[[1]] in your code? You simply copy it from somewhere without understanding it? Anyway, observe the output of the following code: solutionlist = Table[NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2 x[t]-y[t]^3,x[0]==k, y[0]==1}, {x, y}, {t, 0, 1}], {k, 1, 7}]; solutionlist[[3]], then observe solutionlist[[3, 1]], then solutionlist[[3, 1, 1]], then solutionlist[[3, 1, 1, -1]]. $\endgroup$
    – xzczd
    Jul 8, 2022 at 0:13

1 Answer 1

2
$\begingroup$
$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global`*"]

P[r_, z_] := (2 r ((2 - (4 r)/((1 + r)^2 + z^2)) *
        EllipticK[(4 r)/((1 + r)^2 + z^2)] -
       2 EllipticE[(4 r)/((1 + r)^2 + z^2)]))/
   (π ((4 r)/((1 + r)^2 + z^2)) √((1 + r)^2 + z^2));
Zmaxreg = 15;
Zminreg = -Zmaxreg;
Rmaxreg = 54;
Rminreg = 0;

d = 0;
zH = 0;
pref = 1;
a = 4;

sys = {(*EQUATIONS*)
   D[H2[r, z], z] == -(zH2 + d)*H2[r, z]*(P[r, z]/pref)^a,
   D[H[r, z], z] == 
    2*d*H2[r, z]*(P[r, z]/pref)^a - zH*H[r, z]*(P[r, z]/pref)^a,
   (*INITIAL CONDITIONS*)
   H2[r, Zminreg] == 1, H[r, Zminreg] == 0};

Using NDSolve

Clear[H2, H]
Quiet[
  tab = Table[{zH2,
     NDSolve[sys,(*VARIABLES TO SOLVE FOR*)
      {H2, H},
      (*RANGE TO SOLVE OVER*)
      {r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg}]}, 
  {zH2, 1, 10^3, 10}]];

Manipulate[
 Plot3D[Evaluate[{H2[r, z], H[r, z]} /. tab[[n, 2]]],
  {r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg},
  AxesLabel -> (Style[#, 14] & /@ {r, z}),
  PlotLabel -> Style[StringForm["zH2 = ``", tab[[n, 1]]], 14],
  PlotLegends -> (Style[#, 14] & /@ {H2, H})],
 {{n, 61}, 1, Length[tab], 1, Appearance -> "Labeled"}]

enter image description here

Using ParametricNDSolve

sol = ParametricNDSolve[sys,
   (*VARIABLES TO SOLVE FOR*)
   {H2, H},
   (*RANGE TO SOLVE OVER*)
   {r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg},
   (*Parameter*) 
   {zH2}];

Manipulate[
 Plot3D[Evaluate[{H2[zH2][r, z], H[zH2][r, z]} /. sol],
  {r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg},
  AxesLabel -> (Style[#, 14] & /@ {r, z}),
  PlotLegends -> (Style[#, 14] & /@ {H2, H})],
 {{zH2, 601}, 1, 1000, 10, Appearance -> "Labeled"}]

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you for the response! That is helpful. What I'm still not sure about is how to access a specific cell within the table instead of all of them at once. For example, say that I would like to pull only H2 for zH2 = 200. In MatLAB, that would be a simple matter of understanding the shape of the matrix that the data is stored in and calling it. I'm struggling to find the equivalent operation in Mathematica. $\endgroup$
    – Bill
    Jul 7, 2022 at 20:25

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.