I'm trying to model a stellar atmosphere using Mathematica. Basically I'm working in a physical and chemical model which combined should give me a profile of the chemical abundance in the atmosphere. The physical model is based in two equations:

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

which gives me both the temperature and density as a function of the radius. Given those statements I write the differential equations that will give the rate of formation and destruction of the chemical species: $\displaystyle \frac{d(AB)}{dt}=kn(A)n(B)$, where $k$ is the rate coefficient and it depends on the temperature $t$. Now, I have built the code (with your help) to solve the system of ODE and it is working pretty fine:

AU = 1.49597871*10^13; 
To = 2200; 
ΔT = 440; 
Ro = 1.8;
ΔR = 0.1; 
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[ϕ_]:= 2 + 0.25*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]; 
eqns = {nH'[r] == (r17 + r18 + r62) - (r1 + r2), 
   nOH'[r] == (r2 + r17 + r63) - (r1 + r18 + r62 + r94 + r141), 
   nO'[r] == (r1 + r94) - (r62 + r63 + r138), 
   nH2'[r] == (r1 + r2) - (r17 + r18), 
   nH2O'[r] == (r18 + r94) - (r2 + r143 + r144), 
   nTi'[r] == -(r138 + r141 + r143 + r144), 
   nTiO'[r] == (r138 + r141 + r143) - (r144)};

eqEqn = {nH[r] + nOH[r] + nO[r] + nH2[r] + nH2O[r] + nTi[r] + 
     nTiO[r] + nConst[r] == 0.500029396};

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

sol = Table[
   NDSolve[{eqns, eqEqn, ic} /. ϕ -> ϕi, {nH, nOH, nO, nH2, 
     nH2O, nTi, nTiO}, {r, 2, 9}, MaxSteps -> Infinity, 
    AccuracyGoal -> 4, PrecisionGoal -> 4, 
    Method -> {"TimeIntegration" -> "StateSpace"}], {ϕi, 0, 
    2*Pi, Pi/5}];

solTable = DeleteCases[sol, {}];

But I'm missing one thing. I'm not seeing how could I add the density equation

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

to this code.

It should multiply the species abundances n(X) and then the result to be used to solve the ODE. So for example, when $r=2$, $nOH[2]=nOH[0]*dens[2]$. I tried to multiply the initial conditions to the density, but I get a error:

NDSolve::ndinnt: "Initial condition is not a number or a rectangular array of numbers."

Could anyone help me with that? What should I do to incorporate the density equation to solve this problem? Moreover, after the calculation, the final result should be given divided for the density. So, for example, $nOH[2]=nOH[0]*dens[2]$ is used to calculate the ODE, but the final result should be given by $nOH[2]=nOH[2]/dens$. Got it? Is there some way to do that? Thank you very much.

  • $\begingroup$ Your code does not run as written, because equations such as r1 = k1 nH[r] nOH[r]; should be r1 = k1[r, \[Phi]] nH[r] nOH[r];. Even with this correction, introducing dens in the way you describe overspecifies your dependent variables. $\endgroup$
    – bbgodfrey
    Jun 22, 2015 at 5:09

1 Answer 1


As noted in the comment above, k1 etc must be called with its argument list. With this change the r1 etc equations become

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

I interpret the density discussion in the question to mean that expressions of the form

D[species1[r],r] == c species2[r] species3[r]

should be replaced by expressions of the form

D[species1[r] dens[r], r] == c species2[r] dens[r] species3[r] dens[r]

which can be rewritten as

(D[species1[r], r] dens[r] + species1[r] D[dens[r], r])/dens[r]^2 == 
     c species2[r] species3[r]

The operation which performs this transformation is

Table[eqns[[i, 1]] = eqns[[i, 1]] /. Derivative[1][z_][r] -> 
     ((Derivative[1][z][r] + z[r] dldens[r, ϕ])/dens[r, ϕ], {i, 7}];

which yields

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


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, ϕ];


sol = Table[NDSolveValue[{eqns, ic} /. ϕ -> ϕi, {nH, nOH, nO, nH2, nH2O, nTi, nTiO}, 
      {r, 2, 9}, {ϕi, 0, 2*Pi, Pi/5}];

solves the equations without difficulty, and the first element of sol, for instance, can be plotted as

LogPlot[Evaluate[Table[sol[[1, i]][r], {i, 7}]], {r, 2, 9}, 
   PlotLegends -> {nH, nOH, nO, nH2, nH2O, nTi, nTiO}]

enter image description here

Note that nConst[r] is omitted from the calculation for convenience. It can be computed at the end using Solve on eqEqn, if desired.

  • $\begingroup$ Thank you very much @bbgodfrey! But oddly I get an error message when trying to run it: "NDSolveValue::ndsz: At r == 2.`, step size is effectively zero; singularity or stiff system suspected." I tried some adjustments in the code but it didn't work. $\endgroup$ Jun 23, 2015 at 0:24
  • $\begingroup$ My computation runs quickly and without error. Do the equations you feed to NDSolveValue look like (dldens[r, ϕ] nH[r] + nH]'[r])/ dens[r, ϕ] == -r1 + r17 + r18 - r2 + r62? I do know that replacing dens[r, ϕ] by dens[r, ϕ]^2 in these equations produces the error you describe. $\endgroup$
    – bbgodfrey
    Jun 23, 2015 at 0:33
  • $\begingroup$ Yes, precisely like that. I may have misunderstood the operation Table[eqns[[i, 1]] = eqns[[i, 1]] /. Derivative[1][z_][r] -> ((Derivative[1][z][r] + z[r] dldens[r, ϕ])/dens[r, ϕ], {i, 7}];, to transform the expression because everything else looks exactly as you said. $\endgroup$ Jun 23, 2015 at 0:53
  • $\begingroup$ If your equations look like the one in my last comment after applying the Table operation, I would expect that you would obtain the same results that I do. I shall add all seven equations to my answer. $\endgroup$
    – bbgodfrey
    Jun 23, 2015 at 1:13
  • $\begingroup$ Thank you very much for your help @bbgodfrey! $\endgroup$ Jun 23, 2015 at 1:24

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.