Trying to model a chemical reaction - NDSolve::litarg error

I am trying to model a chemical reaction with NDSolve.

I have specified the known differential equations and parameters,

With[{kH = 9.63*10^-5, kR = 0., k = .01, cE0 = 0.25*10^-3},
NDSolve[{
cEOH'[t] == kH*(cE0 - cB[0] - cEOH[t] + cB[t]),
cB'[t] == -k*cB[t]*(cE0 - cB[0] - cEOH[t] + cB[t]),
cR'[t] == kR*(cB[0] - (cB[t] + cR[t])),
cEOH[0] == 0.,
cB[0] == 1.624*10^-3,
cR[0] == 0.
},
{cEOH, cB, cR},
{t, 0, 3600*4}
]
]


however, Mathematica can't solve it:

NDSolve::litarg: To avoid possible ambiguity, the arguments of the dependent
variable in {(cEOH^′)[t]==0.0000963 (0.00025 -cB[0]+cB[t]-cEOH[t]),
(cB^′)[t]==-0.01 cB[t] (0.00025 -cB[0]+cB[t]-cEOH[t]),
(cR^′)[t]==0.} should literally match the independent variables. >>


I don't understand the error and the related help article did not clarify it for me. How can I fix the equations?

I replaced cB[0] in the equations with cB0, since I thought the ambiguity comes from that. (Of course I also gave the corresponding value to cB0 in With's first argument.) That seems to work:

 With[{kH = 9.63*10^-5, kR = 0., k = .01, cE0 = 0.25*10^-3, cB0 = 1.624*10^-3},
NDSolve[{
cEOH'[t] == kH*(cE0 - cB0 - cEOH[t] + cB[t]),
cB'[t] == -k*cB[t]*(cE0 - cB0 - cEOH[t] + cB[t]),
cR'[t] == kR*(cB0 - (cB[t] + cR[t])),
cEOH[0] == 0.,
cB[0] == cB0,
cR[0] == 0.}, {cEOH, cB, cR}, {t, 0, 3600*4}]]


A simpler system to get the message in the question is:

NDSolve[{
cR'[t] == (1 + cR[t]),
cR[0] == 1}, {cR}, {t, 0, 1}]


I would conjecture that an explanation for it is related to DSolve finding this solution (note "cR[0]"):

DSolve[{
cR'[t] == (cR[0] + cR[t]),
cR[0] == 1}, {cR}, {t, 0, 1}]

(* Out[37]= {{cR -> Function[{t}, -cR[0] + E^t (1 + cR[0])]}} *)


When you tell NDSolve "{cEOH, cB, cR}" is the list of independent functions, it requires that you always use them as functions, not as constants. I.e., cB[t] is fine (if t is on the list of independent variables), but cB[0] is not a function, it's a constant.

This is largely because of internal representation and parsing: When the expression cB[0] is parsed, the "cB" is seen first, found on the dependent variables list, then its argument "0" is not found on the independent variables list and much hilarity ensues. The usage in specifying initial conditions is caught as a special case and must be of the form {indepVar}[{constExpr}]=={constExpr} and likewise for periodic bounadary conditions.

The documentation says all this (not very clearly). The problem you're bumping into is "covered" by "In ordinary differential equations, the function $u_i$ must depend only on the single variable $t$. In partial differential equations, they may depend on more than one variable."

As an added and related bonus, try "f'[x] == f'[t]". E.g., (and never mind that this system is hopeless)

NDSolve[{
f'[x] == f'[t],
f[0] == 0
}, {f}, {x, 0, 1}, {t, 0, 1}]
(* NDSolve::derlen: "The length of the derivative operator Derivative[1] in f'[x] is not the same as the number of arguments." *)


Or to rephrase in "comprehensible-ese" -- "the function $f$ doesn't depend on all the independent variables, so I'm going to sulk". NDSolve[] is written to assume that the domain of definition is some (multidimensional) rectangle and all dependent functions will be defined on that rectangle (or on its boundary sub-rectanges) and all uses of them must have the right variable lists.

Believe it or not, this is a good thing -- it prevents attacking the halting-problem-equivalent problem of detecting things which are really D[f[x, x], x] no matter how convolutedly they are actually expressed. (Really, partial derivatives are defined in terms of slots -- for example "the derivative of $f(\cdot,\cdot)$ as the first argument is varied" -- but that's not the notation we use, in which notation $\frac{\partial}{\partial x} f(x,x)$ is not syntactically valid.")

Sorry to ramble on... The short version is that constants should look like constants, not dependent function invocations. (This is adequately demonstrated in the other answer(s).)

• This is a good research on the subject, but I think the first two paragraphs are not clear. More specifically, in the code of the question: (1) cB[0] is using cB as a function (as it is advised in the first paragraph), and (2) the initial conditions are following the rule "{indepVar}[{constExpr}]=={constExpr}" (as advised in the second paragraph.) – Anton Antonov Oct 3 '15 at 22:39
• @AntonAntonov : No. "cB[0]" is a constant, it is whatever constant initial value the function has. It is not something that can be meaningfully differentiated and the independent variables do not appear explicitly. – Eric Towers Oct 3 '15 at 22:41
• @AntonAntonov : Or perhaps you don't realize that NDSolve[] is complaining about the usage of, for example, "cB[0]" in "cE0 - cB[0] - cEOH[t] + cB[t]", which is not an initial condition usage. This is quite evident in the OP's quoted error message. – Eric Towers Oct 3 '15 at 22:43
• "No. "cB[0]" is a constant,[...]" -- of course it is. Read the first paragraph of your explanation -- it says: "When you tell NDSolve "{cEOH, cB, cR}" is the list of independent functions, it requires that you always use them as functions, not as constants." What I am saying is that "cB[0]" is using "cB" as a function, so your explanation is not sufficient. – Anton Antonov Oct 3 '15 at 23:25
• @AntonAntonov : If cB[0] were a function, we could differentiate with respect to its nonexistent parameter. It doesn't have one, we can't, and it's not a function invocation. I know it "looks" like one and I "get" that you wish it were. Too bad. – Eric Towers Oct 4 '15 at 4:39

OK, I solved this and it looks so stupid that I am considering deleting the question.

The solution is to replace all instances of cB[0] in the differential equations with the actual value 1.624*10^-3:

...
NDSolve[{
cEOH'[t] == kH*(cE0 - 1.624*10^-3 - cEOH[t] + cB[t]),
cB'[t] == -k*cB[t]*(cE0 - 1.624*10^-3 - cEOH[t] + cB[t]),
cR'[t] == kR*(1.624*10^-3 - (cB[t] + cR[t])),
...


This is very counterintuitive and I don't know why it doesn't work as specified in the OP.

• It does seem strange. I would be inclined to keep the question, since hopefully someone might be able to provide a justification for this behavior. – Oleksandr R. Oct 3 '15 at 15:08
• @shrx Please do not delete the question. I was puzzled by that too... – Anton Antonov Oct 3 '15 at 15:16
• Thanks for your comments, I will keep it. – shrx Oct 3 '15 at 15:26

I just had a similar problem where I wanted to have some output variable $y(t)$ be an index of a stock variable, e.g. $y(t) = \frac{x(t)}{x(0)}$. As with the OP's cB[0] $x(0)$ actually is a constant.

The general solution in this case could be to declare a separate stock variable as a "memory", e.g. make it a discrete variable (DiscreteVariables) and not change it or only change it at certain times using WhenEvent:

With[
{
kH = 9.63*10^-5,
kR = 0.,
k = .01,
cE0 = 0.25*10^-3
},
NDSolve[
{
(* 'flow' equations *)
cEOH'[t] == kH*(cE0 - cB0[t] - cEOH[t] + cB[t]),
cB'[t] == -k*cB[t]*(cE0 - cB0[t] - cEOH[t] + cB[t]),
cR'[t] == kR*(cB0[t] - (cB[t] + cR[t])),
(* initial stock values *)
cEOH[0] == 0.,
cB[0] == 1.624*10^-3,
cR[0] == 0.,
cB0[0] == cB[0]
},
{ cEOH, cB, cR, cB0 }, (* dependent variables *)
{ t, 0, 3600 × 4},        (* independent variables *)

DiscreteVariables -> {cB0}

]
]


This works out fine.

Here we introduced another dependent variable cB0[t] which is initialized as cB0[0] == cB[0]. We have declared it to be a discrete variable and we then never change it. So cB[t] is a valid dependent variable on t which can be used meeting the requirements in the flow equations but it will only refer to cB[0] and thus behave like a constant.

• In this simple case we might of course have declared cB0 simply as a constant, as AntonAntonov has done in his answer, but the solution here is more general and will work if cB[0] is to be some more complex function of other stocks. – gwr May 6 '16 at 13:12