4
$\begingroup$

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[2] == 10^-7, nOH[2] == 4.9*10^-7, nO[2] == 10^-5, nH2[2] == 5*10^-1, nH2O[2] == 4*10^-6, nCONS[2] == 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?
I hope I made myself clear. Thank you in advance!

$\endgroup$
3
  • 2
    $\begingroup$ 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? $\endgroup$
    – MarcoB
    May 26, 2015 at 4:51
  • 1
    $\begingroup$ 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) $\endgroup$ May 26, 2015 at 11:19
  • $\begingroup$ 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? $\endgroup$
    – Phab
    May 26, 2015 at 12:36

1 Answer 1

3
$\begingroup$

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

I took your code:

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[2] == 10^-7, nOH[2] == 4.9*10^-7, 
  nO[2] == 10^-5, nH2[2] == 5*10^-1, nH2O[2] == 4*10^-6, 
  nCONS[2] == 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,ϕ]"]

enter image description here

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ May 26, 2015 at 13:32
  • $\begingroup$ @AndréOliveira In what way should it change? Do you want a table with a solution for each ϕ? Or should ϕ change with r? $\endgroup$
    – Phab
    May 27, 2015 at 6:38
  • $\begingroup$ Just like you did. Thank you very much @Phab! $\endgroup$ May 27, 2015 at 14:04

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.