Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Many dynamic processes can occur simultaneously, giving more than one ODE. For example, consider two tanks. Tank 1 has concentration c1c_1 of a solute and Tank 2 has concentration c2c_2. Tank 1 flows into Tank 2, which is then passed through a pump and recycled back to Tank 1.

Concentrations in two tanks

Eventually, the two concentrations should equalize, but how do they evolve to over time? c1c_1 and c2c_2 can be modeled using unsteady mass balances. This is a system of first-order ODEs that can be solved simultaneously.

The general, explicit form of a system of first-order ODEs is:

y1=f1(t,y1,,yn)yn=fn(t,y1,,yn)\begin{align} y_1' &= f_1(t, y_1, \cdots, y_n)\\ \vdots \\ y_n' &= f_n(t, y_1, \cdots, y_n) \end{align}

where t is the independent variable, (y1,,yn)(y_1, \cdots, y_n) are the n dependent variables, and (f1,,fn)(f_1, \cdots, f_n) are the n right-hand side functions of the n ODEs. We will notate this system more compactly using vectors

y=f(t,y)\vv{y}' = \vv{f}(t, \vv{y})

where y is the column vector of dependent variables, f is the column vector of right-hand side functions, and y\vv{y}' is a shorthand for the column vector of first derivatives of all dependent variables.

Successive substitution

Some systems of first-order ODEs have a special structure that allows you to solve one equation at a time, then substitute into another that can now be solved. This process can be repeated successively to solve for all dependent variables.

Example: Reaction network

We are analyzing a simple reaction sequence occuring in a batch reactor:

Ak1Bk2C{\rm A} \xrightarrow{k_1} {\rm B} \xrightarrow{k_2} {\rm C}

where k1k_1 and k2k_2 are reaction rate constants. At t=0t = 0, only species A is present at concentration cA,0c_{{\rm A},0}. Determine the concentration of B as a function of time. What happens in the special case k1=k2k_1 = k_2?


Unsteady mole balances give

d ⁣cAd ⁣t=k1cAd ⁣cBd ⁣t=k1cAk2cBd ⁣cCd ⁣t=k2cB\begin{align} \dd{}{c_{\rm A}}{t} &= -k_1 c_{\rm A} \\ \dd{}{c_{\rm B}}{t} &= k_1 c_{\rm A} - k_2 c_{\rm B} \\ \dd{}{c_{\rm C}}{t} &= k_2 c_{\rm B} \end{align}

These ODES can be solved by first finding cAc_{\rm A}, then cBc_{\rm B}, and last cCc_{\rm C}! The ODE for cAc_{\rm A} is separable:

d ⁣cAcA=k1d ⁣tlncA=k1t+K\begin{align} \int \frac{\d{c_{\rm A}}}{c_{\rm A}} &= \int -k_1 \d{t} \\ \ln c_{\rm A} &= -k_1 t + K \end{align}

Applying the initial condition cA(0)=cA,0c_{\rm A}(0) = c_{{\rm A},0} gives:

cA=cA,0ek1tc_{\rm A} = c_{{\rm A}, 0} e^{-k_1 t}

Next, use the solution for cAc_{\rm A} to solve cBc_{\rm B}:

d ⁣cBd ⁣t=k1cA,0ek1tk2cBd ⁣cBd ⁣t+k2cB=k1cA,0ek1t\begin{align} \dd{}{c_{\rm B}}{t} &= k_1 c_{\rm A,0}e^{-k_1 t} -k_2 c_{\rm B} \\ \dd{}{c_{\rm B}}{t} &+ k_2 c_{\rm B} = k_1 c_{\rm A,0}e^{-k_1 t} \end{align}

This is a linear first-order ODE with

p=k2r=k1cA,0ek1tp = k_2 \qquad r =k_1 c_{\rm A,0}e^{-k_1 t}

that can be solved using an integrating factor:

F=ek2d ⁣t=ek2tFrd ⁣t=ek2tk1cA,0ek1td ⁣t=k1cA,0e(k2k1)td ⁣t=cA,0k1k2k1e(k2k1)t\begin{align} F &= e^{\int k_2 \d{t}} = e^{k_2 t} \\ \int F r \d{t} &= \int e^{k_2 t} k_1 c_{\rm A,0} e^{-k_1 t} \d{t} \\ &= k_1 c_{\rm A,0} \int e^{(k_2 -k_1)t} \d{t} \\ &= c_{\rm A,0} \frac{k_1}{k_2 -k_1} e^{(k_2 -k_1)t} \end{align}

Hence,

cB=1F(Frd ⁣t+K)=cA,0k1k2k1ek1t+Kek2t\begin{align} c_{\rm B} &= \frac{1}{F}\left(\int Fr\d{t} + K \right) \\ &= c_{\rm A,0} \frac{k_1}{k_2 - k_1} e^{-k_1 t} + K e^{-k_2 t} \end{align}

Applying the initial condition

cB(0)=cA,0k1k2k1+K=0K=cA,0k1k2k1\begin{align} c_{\rm B}(0) &= c_{\rm A,0} \frac{k_1}{k_2 - k_1} + K = 0 \\ K &= -c_{\rm A,0} \frac{k_1}{k_2 - k_1} \end{align}

Hence,

cB=cA,0k1k2k1(ek1tek2t)c_{\rm B} = c_{\rm A,0} \frac{k_1}{k_2 -k_1} \left(e^{-k_1 t} - e^{-k_2 t} \right)

If k1=k2k_1 = k_2, the solution appears to have an issue because there is a “divide by zero”. However, we can actually evaluate this behavior as a limit! Let’s set k1=kk_1 = k, and take the limit as k2kk_2 \to k:

limk2kcB=cA,0klimk2kektek2tk2k1=cA,0klimk2ktek2t1=cA,0ktekt\begin{align} \lim_{k_2 \to k} c_{\rm B} &= c_{{\rm A},0} k \lim_{k_2 \to k} \frac{e^{-k t} - e^{-k_2 t}}{k_2 -k_1} \\ &= c_{{\rm A},0} k \lim_{k_2 \to k} \frac{t e^{-k_2 t}}{1} \\ &= c_{{\rm A},0} k t e^{-k t} \end{align}

where in the second line, we made use of L’Hopital’s rule to evaluate the limit.

Skill builder problems

Solution to Exercise 1

We begin by solving for y1y_1 in terms of x. The ODE is separable and therefore can be solved using separation of variables.

dy112y1=dx12ln(12y1)=x+c0ln(12y1)=2x+c112y1=c2e2xy1=12(1c2e2x)\begin{align} \int \frac{dy_1}{1-2y_1} &= \int dx \\ -\frac{1}{2} \ln{\left(1-2y_1\right)} &= x + c_0 \\ \ln{\left(1-2y_1\right)} &= -2x + c_1 \\ 1 - 2y_1 &= c_2 e^{-2x} \\ y_1 &= \frac{1}{2} \left(1- c_2 e^{-2x}\right) \\ \end{align}

Now, use the initial condition y1(0)=0y_1(0)= 0 to solve for c2c_2:

y1(0)=12(1c2)=0y_1(0) = \frac{1}{2}\left(1-c_2\right) = 0

so c2=1c_2 = 1 and

y1=12(1e2x)y_1 = \frac{1}{2} \left(1-e^{-2x}\right)

Now, we can solve for y2y_2 knowing y1y_1 and recognizing that y2y_2 can also be solved via separation of variables

y2=y1(2y2)dy22y2=12(1e2x)d ⁣xln(2y2)=12(x+12e2x)+c02y2=c1exp[12(x+12e2x)]y2=2c1exp[12(x+12e2x)]\begin{align} y_2' &= y_1\left(2-y_2\right) \\ \int \frac{dy_2}{2-y_2} &= \int \frac{1}{2}\left(1-e^{-2x}\right) \d{x} \\ -\ln{\left(2-y_2\right)} &= \frac{1}{2}\left(x + \frac{1}{2}e^{-2x}\right) + c_0\\ 2 - y_2 &= c_1 \exp\left[-\frac{1}{2}\left(x+ \frac{1}{2}e^{-2x}\right)\right] \\ y_2 &= 2 - c_1 \exp\left[-\frac{1}{2}\left(x+\frac{1}{2}e^{-2x}\right)\right] \\ \end{align}

Again, use the initial condition y2(0)=0y_2(0)= 0 to solve for c1c_1:

y2(0)=2c1e1/4=0c1=2e1/4\begin{align} y_2(0) &= 2 - c_1 e^{-1/4} = 0 \\ c_1 &= 2e^{1/4} \end{align}

Therefore,

y2=2[1exp(14x214e2x)]y_2 = 2\left[1-\exp\left(\frac{1}{4}-\frac{x}{2}-\frac{1}{4}e^{-2x}\right)\right]