1
$\begingroup$

Normally I wouldn't ask such question before even read the material from the Wolfram Documentation about what I'm trying to do but this time I have a pretty tigh dead end. I apologize for that. So, what I need to do is to extract a data from a NDSolve integration to use it in different ways. I will show my whole code and make a precise description for what I need. The code is the following:

AU = 1.49597871*10^13;
n = 0.5*10^15; 
cg = 6.67384/10^8; 
kb = 1.3806488/10^16; 
mh = 1.0078250321*(1.660538921/10^24); 
mf = 1.98892*10^33; 
Te[ϕ_] := 2200 + 440*Cos[ϕ]; 
Ra[ϕ_] := 1.9 + 0.1*Cos[ϕ + Pi]; 
t[r_, ϕ_] := Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4);

k1[r_, ϕ_] := (6.99/10^14)*(t[r, ϕ]/300)^2.8*Exp[-1950/t[r, ϕ]]; 
k2[r_, ϕ_] := (1.59/10^11)*(t[r, ϕ]/300)^1.2*Exp[-9610/t[r, ϕ]]; 
k17[r_, ϕ_] := (3.14/10^13)*(t[r, ϕ]/300)^2.7*Exp[-3150/t[r, ϕ]]; 
k18[r_, ϕ_] := (2.05/10^12)*(t[r, ϕ]/300)^1.52*Exp[-1736/t[r, ϕ]]; 
k62[r_, ϕ_] := (1.77/10^11)*Exp[178/t[r, ϕ]]; 
k63[r_, ϕ_] := (1.85/10^11)*(t[r, ϕ]/300)^0.95*Exp[-8571/t[r, ϕ]]; 
k94[r_, ϕ_] := (1.65/10^12)*(t[r, ϕ]/300)^1.14*Exp[-50/t[r, ϕ]]; 
k138[r_, ϕ_] := (5.94/10^17)*(t[r, ϕ]/300)^0.17*Exp[65.9/t[r, ϕ]]; 
k141[r_, ϕ_] := (1/10^11)*Exp[-4800/t[r, ϕ]]; 
k143[r_, ϕ_] := (1/10^11)*Exp[-4800/t[r, ϕ]]; 
k144[r_, ϕ_] := 1/10^9; 

r1 = k1[r, ϕ]*nH[r]*nOH[r]; 
r2 = k2[r, ϕ]*nH[r]*nH2O[r]; 
r17 = k17[r, ϕ]*nH2[r]*nO[r]; 
r18 = k18[r, ϕ]*nH2[r]*nOH[r]; 
r62 = k62[r, ϕ]*nO[r]*nOH[r]; 
r63 = k63[r, ϕ]*nO[r]*nH2O[r]; 
r94 = k94[r, ϕ]*nOH[r]*nOH[r]; 
r138 = k138[r, ϕ]*nTi[r]*nO[r]; 
r141 = k141[r, ϕ]*nTi[r]*nOH[r]; 
r143 = k143[r, ϕ]*nTi[r]*nH2O[r]; 
r144 = k144[r, ϕ]*nTiO[r]*nH2O[r];

dens[r_, ϕ_] := n*Exp[(-((2*cg*mf*mh)/(kb*Te[ϕ]*Ra[ϕ]*AU)))*
           (1 - (Ra[ϕ]/r)^(1/2))]; 
dldens[r_, ϕ_] := D[dens[r, ϕ], r]/dens[r, ϕ];

eqns = {
   (dldens[r, ϕ]*nH[r] + Derivative[1][nH][r])/dens[r, ϕ] == 
     -r1 + r17 + r18 - r2 + r62, 
   (dldens[r, ϕ]*nOH[r] + Derivative[1][nOH][r])/dens[r, ϕ] == 
     -r1 - r141 + r17 - r18 + r2 - r62 + r63 - r94, 
   (dldens[r, ϕ]*nO[r] + Derivative[1][nO][r])/dens[r, ϕ] == 
      r1 - r138 - r62 - r63 + r94, 
   (dldens[r, ϕ]*nH2[r] + Derivative[1][nH2][r])/dens[r, ϕ] == 
      r1 - r17 - r18 + r2, 
   (dldens[r, ϕ]*nH2O[r] + Derivative[1][nH2O][r])/dens[r, ϕ] == 
     -r143 - r144 + r18 - r2 + r94, 
   (dldens[r, ϕ]*nTi[r] + Derivative[1][nTi][r])/dens[r, ϕ] == 
     -r138 - r141 - r143 - r144, 
   (dldens[r, ϕ]*nTiO[r] + Derivative[1][nTiO][r])/dens[r, ϕ] == 
      r138 + r141 + r143 - r144}; 

ic = {nH[2] == 2/10^7, nOH[2] == 9.8/10^7, nO[2] == 2/10^5, nH2[2] == 5/10, 
      nH2O[2] == 8/10^6, nTi[2] == 1.96/10^7, nTiO[2] == 2/10^8}; 

sol = Table[
   NDSolveValue[{eqns, ic} /. ϕ -> ϕi, {nH, nOH, nO, nH2, nH2O, nTi, nTiO},            
      {r, 2, 6}, MaxStepSize -> 0.001, MaxSteps -> Infinity,
                 AccuracyGoal -> 0, PrecisionGoal -> 0, Method -> {"BDF"}], 
      {ϕi, 0, 2*Pi, Pi/5}];

Question: how could I extract the data from the solution to a especific value of r and ϕ?

Let's say, for example, that I want to know nOH[6] for all values of ϕ, so

nOH[6] for ϕ=0
nOH[6] for ϕ=Pi/5
nOH[6] for ϕ=Pi

and so on. Is there some way to do it?

Moreover, is there some way to export the data from the solution where I can obtain an arrangement like the following?

Table

Thank you very much.

$\endgroup$
  • $\begingroup$ I get several NDSolveValue::ndnum errors when I execute your code. Perhaps you forgot to post something? $\endgroup$ – Michael E2 Dec 21 '15 at 20:46
  • $\begingroup$ maybe it should be cg*mf*mh ..? $\endgroup$ – george2079 Dec 21 '15 at 21:16
  • $\begingroup$ I think this is what you want: TableForm[ Table[{r, sol[[1, 2]][r]} , {r, {2, 2.1, 2.2, 2.3}}] , TableHeadings -> {None, {r, y}}] $\endgroup$ – george2079 Dec 21 '15 at 21:29
  • $\begingroup$ @garej yes, it look reasonable! And yes, cgmfmh= cg * mf * mh. I'm sorry, didn't see it got the messy when I wrote it here. $\endgroup$ – André Oliveira Dec 21 '15 at 21:30
  • 2
    $\begingroup$ I get sol[[1, 2]][6] -> 258767. Andre, you should edit the question to fix that cgmfmh typo. $\endgroup$ – george2079 Dec 21 '15 at 21:33
6
$\begingroup$

In your case you get the matrix of solutions. You can extract them by parts like sol[[2,1]][6] but it is not convenient since you have to think which number stands for what combination of variables: sol[[phi, element]][r].

It is comfortable to set up a Manipulate for the task. There are of course many ways to do that.

Manipulate[
 sol;
 expr = Thread[List[Range[up - 1, up, 0.1], Table[sol[[phi, element]][r], 
        {r, Range[up - 1, up, 0.1]}]]];

 export = Grid[{
      (* for Export first row is better to make as {"\[Phi]=", phi/.phirange} *)
      {Item["\[CapitalPhi]=" Dynamic[phi /. phirange], 
        Alignment -> {Left, Center}, Frame -> {{True, False}}], Nothing},
      {"r", element /. list},
      ##
      }, Frame -> All, ItemSize -> {{10, 10}, {3, {1.6}}}] & @@ expr,

 {{up, 3, "Range r"}, 3, 6, 0.1, Appearance -> "Labeled"},
 Delimiter,
 {{element, 2}, list, SetterBar},
 Delimiter,
 {{phi, 1}, phirange, SetterBar},

 Delimiter,
 (* change Print[export] for Export["file__.xls",export]*)
 Button["Export", Print[export]],

 ControlPlacement -> Left,

 Initialization :> (
   list = 
    Thread[Range[7] -> {"nH", "nOH", "nO", "nH2", "nH2O", "nTi", "nTiO"}];
   phirange = 
    Thread[Range[Length@Range[0, 2 Pi, Pi/5]] -> Range[0, 2 Pi, Pi/5]];
   )
]

Result: enter image description here


Addenum: With solution as sol some, I would recommend:

Manipulate[
 sol;(*Dimensions[sol] gives {11,7}; elements of sol are function heads*)
 expr = 
  TableForm[
   Table[sol[[phi, element]][r], 
   {r, Range[2, up, 0.1]}, 
   {phi, Range[Length@Range[0, p, Pi/5]]}], 
   TableHeadings -> {Range[2, 6, 0.1], Range[0, 2 Pi, Pi/5]}
  ];
 Framed[expr]
 ,
 Button["Print", Print[expr]], 
 (* you may change Print[expr] for Export["file.ext", expr]*)

 {{element, 1}, {1 -> "nH", 2 -> "nOH", 3 -> "nO", 4 -> "nH2", 
   5 -> "nH2O", 6 -> "nTi", 7 -> "nTiO"}}, 
 {{up, 3, "Range r"}, 3, 6, 0.1}, 
 {{p, Pi/5, "Range phi"}, Pi/5, 2 Pi, Pi/5}]

Special thanks to @george2079 comments.

enter image description here

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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