# Solving a chemical system using NDSolve

I have to solve a set of coupled differential equations and am not sure how to proceed. It is a model of a chemical process where I'm interested to find the species concentration as $r$ varies.
For a chemical reaction given by $A + B \rightarrow AB + hv$,

$\displaystyle \frac{d(AB)}{dr}=k*n(A)*n(B)$, is the rate of production of the species $AB$

So far I have built the following code

Te[ϕ_]:= 2200 + 440*Cos[ϕ]
Ra[ϕ_]:= 2 + 0.25*Cos[ϕ + Pi]
t[r_,ϕ_]:= Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4)
Note that, for my case, $k_{i}$ also changes with ϕ (to make it simple to solve, I can set ϕ manually from 0 to 2π)

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

Now I should build the rate equations, based on the differential equantion I pointed in the begining. And there is where I'm not sure how to proceed. The Wolfram Documentation has a good example about how to work with chemical systems, but there the $k_i$ are constants which make it much more simple to solve. There, the rate equations would be write:
r1 = k1 nH[r] nOH[r];
r2 = k2 nH[r] nH2O[r];
r17 = k17 nH2[r] nO[r];
r18 = k18 nH2[r] nOH[r];
r62 = k62 nO[r] nOH[r];
r63 = k63 nO[r] nH2O[r];
r94 = k94 nOH[r] nOH[r];
As the $k's$ in my case are dependent on $r$ and $ϕ$, should I write it like
r1[r_,ϕ_] = k1[r,ϕ]*nH[r]*nOH[r];
for each rate equation?

Just as an example, I wrote the code to solve the system where $k_{i}$ are constants:
r1 = k1 nH[r] nOH[r];
r2 = k2 nH[r] nH2O[r];
r17 = k17 nH2[r] nO[r];
r18 = k18 nH2[r] nOH[r];
r62 = k62 nO[r] nOH[r];
r63 = k63 nO[r] nH2O[r];
r94 = k94 nOH[r] nOH[r];
eqns = {nH'[r] == (r17 + r18 + r62) - (r1 + r2),
nOH'[r] == (r2 + r17 + r63) - (r1 + r18 + r62 + r94),
nO'[r] == (r1 + r94) - (r62 + r63),
nH2'[r] == (r1 + r2) - (r17 + r18), nH2O'[r] == (r18 + r94) - (r2)};

eqEqn = {nH[r] + nOH[r] + nO[r] + nH2[r] + nH2O[r] + nCONS[r] ==
0.50001459};
ic = {nH == 10^-7, nOH == 4.9*10^-7, nO == 10^-5,
nH2 == 5*10^-1, nH2O == 4*10^-6, nCONS == 0};

params = {k1 -> 7.1143*10^-13, k2 -> 2.41635*10^-12,
k17 -> 1.744057*10^-11, k18 -> 2.00061*10^-11,
k62 -> 1.91656*10^-11, k63 -> 2.70681*10^-12,
k94 -> 1.594279*10^-11 };

sol = NDSolve[{eqs, igEqs, ic} /. params, {nH, nOH, nO, nH2,
nH2O}, {r, 2, 8}];
and it worked pretty fine. How could I rewrite this code to work with $k$ changing with $r$? Moreover, is it possible to relate the independent variable $r$ that I used both in the equations and NDSolve? For example, in the NDSolve I've used $r=2$, as the first step. The program will use it to evaluate the previous equations
t[2,ϕ_]:= Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/2)^2])^(1/4);
k1[2,ϕ_]:= (6.99*10^-14)*((t[2,ϕ]/300)^2.8)*Exp[-1950/t[2,ϕ]]
k2[2,ϕ_]:= (1.59*10^-11)*((t[2,ϕ]/300)^1.2)*Exp[-9610/t[2,ϕ]]
k17[2,ϕ_]:= (3.14*10^-13)*((t[2,ϕ]/300)^2.7)*Exp[-3150/t[2,ϕ]]
k18[2,ϕ_]:= (2.05*10^-12)*((t[2,ϕ]/300)^1.52)*Exp[-1736/t[2,ϕ]]
k62[2,ϕ_]:= (1.77*10^-11)*Exp[178/t[2,ϕ]]
k63[2,ϕ_]:= (1.85*10^-11)*((t[2,ϕ]/300)^0.95)*Exp[-8571/t[2,ϕ]]
k94[2,ϕ_] := (1.65*10^-12)*((t[2,ϕ]/300)^1.14)*Exp[-50/t[2,ϕ]]  
and use the results on NDSolve? If don't, is there some other way to appoach this problem?

• Out of curiosity, What chemical process are you modeling? If I am not misunderstanding, your rate constants change in the course of your process. What causes the change? – MarcoB May 26 '15 at 4:51
• I'm trying to model the chemical abundance of a Mira star as a function of its radius. The temperature t[r_,ϕ_] changes as a function of the star radius and so the rate constants k, as it depends exponentially on the temperature k = α*((t/300)^β)*(Exp[^- γ/T) – André Oliveira May 26 '15 at 11:19
• I may not understand your problem, but why don't you just write your r eqns as r1 = k1[r, \[Phi]] nH[r] nOH[r] etc? – Phab May 26 '15 at 12:36

I think you just have to rewrite your "r" equations.

Te[ϕ_] := 2200 + 440*Cos[ϕ]
Ra[ϕ_] := 2 + 0.25*Cos[ϕ + π]

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


Here I only changed your code the way your ks depend on r and ϕ, like you defined them above:

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


With your equations and boundary conditions:

eqns = {nH'[r] == (r17 + r18 + r62) - (r1 + r2),
nOH'[r] == (r2 + r17 + r63) - (r1 + r18 + r62 + r94),
nO'[r] == (r1 + r94) - (r62 + r63),
nH2'[r] == (r1 + r2) - (r17 + r18), nH2O'[r] == (r18 + r94) - (r2)};

eqEqn = {nH[r] + nOH[r] + nO[r] + nH2[r] + nH2O[r] + nCONS[r] ==
0.50001459}; ic = {nH == 10^-7, nOH == 4.9*10^-7,
nO == 10^-5, nH2 == 5*10^-1, nH2O == 4*10^-6,
nCONS == 0};


This gives at least a numerical solution. Note, that I defined ϕ to be 2π (you mentioned you can change it manually):

sol = NDSolve[{eqns, eqEqn, ic} /. ϕ -> 2π, {nH,
nOH, nO, nH2, nH2O}, {r, 2, 8}]


In addition: If you not want to change ϕ manually, just generate a table over your different ϕ (note: I didnt get a solution for all of them):

solTable = Table[
NDSolve[{eqns, eqEqn, ic} //. ϕ -> ϕi, {nH, nOH, nO, nH2,
nH2O}, {r, 2, 8}],
{ϕi, 0, 2π, π/10}];

(* delete empty entrys *)
shortSolTable = DeleteCases[solTable, {}];


To visualise your solutions you can use:

Table[nH[r] /. shortSolTable[[i, 1]], {i, 1, Length[shortSolTable]}];
Plot[%, {r, 2, 8}, PlotLegends -> Automatic,PlotLabel -> "nH[r,ϕ]"] • Thank you Phab! I'm a beginner in the Mathematica programming so am using this problem to learn as much as I can. It didn't occur to me that the solution would be so simple. By the way, what if I wanted to change ϕ automatically? How should I write it? Let's say I'd like ϕ to change from 0 to 2π with a step of π/10. – André Oliveira May 26 '15 at 13:32
• @AndréOliveira In what way should it change? Do you want a table with a solution for each ϕ? Or should ϕ change with r? – Phab May 27 '15 at 6:38
• Just like you did. Thank you very much @Phab! – André Oliveira May 27 '15 at 14:04