# Solution to the rate law of first order reaction followed by second order reaction (consecutive)

I'm trying to derive the rate law of a first order followed by second order consecutive reaction:

$$\frac{\mathrm{d}[\ce{A}]}{\mathrm{d}t}=-k_1[\ce{A}]$$ $$\frac{\mathrm{d}[\ce{B}]}{\mathrm{d}t}=k_1[\ce{A}]-k_2{[\ce{B}]^2}$$ $$[\ce{A}](0)=[\ce{A}]_0,\ce{[B]}(0)=0$$

The solution for A is trivial: $$\ce{[A]}=\ce{[A]}_0\exp(-k_1t)$$

The problem is B, whose rate law ODE is not linear: $$\frac{\mathrm{d}[\ce{B}]}{\mathrm{d}t}=k_1[\ce{A}]_0\exp(-k_1t)-k_2{[\ce{B}]^2}$$

This is a nonlinear first order ODE, but can be rewritten to a linear second order ODE by substituting [\ce{B}] with: $$[\ce{B}]=\frac{\mathrm{d}y}{\mathrm{d}t}\frac{1}{k_2y}$$

Which gives: $$-\frac{\mathrm{d}^2y}{\mathrm{d}t^2}+k_1k_2[\ce{A}]_0\exp(-k_1t)y=0$$

The solution is quite complicate and contains bessel functions of both the first and the second kind:

$$[B] = \frac{C_1\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,J_1\left(\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-\ce{A}_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)\,\sqrt{-A_{0}\,k_{1}\,k_{2}}-\ce{C}_2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,Y_1\left(-\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{2}\,\left(C_1\,J_0\left(\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)+C_2\,Y_0\left(-\frac{2\,{\mathrm{e}}^{-\frac{k_{1}\,t}{2}}\,\sqrt{-A_{0}\,k_{1}\,k_{2}}}{k_{1}}\right)\right)}$$

There are two arbitrary constants but I only have one condition: $$\ce{[B]}(0)=0$$

So my questions are:

• Is my approach solving the rate law correct ?
• If so, how do I get rid of both constants ?
• Is there any approximated solution to this rate law ?

[update]: The ode is solved numerically with Runge-Kutta method. The symbolic solution seems correct.

• Unless you really want an algebraic solution you can numerically solve $dB/dt$ quite easily and you know both initial conditions; at $t=0, [A]=[A_0],[B]=0$. Once you make a 2nd order eqn, you will have 2 initial conditions making 3 in all. I'm not sure what this could be, $[B]=0, t=\infty$ is probably not helpful, so this means knowing $dB/dt$ at some time. Your substitution is a cunning one :), its a pity it does not make a Bernoulli eqn. – porphyrin Jun 25 '18 at 21:57
• @porphyrin I don't see how numerical solution is easy because of the two parameters k1 and k2. To solve that one needs to nest the numerical ode solver inside optimization iterations. The optimization of k1 and k2 will be a global fit of a 2D array. – 7E10FC9A Jun 26 '18 at 15:26
• Of course you have to give k1 and k2 values which you will also have to do for your algebraic solution when you want to plot it. The form of the chemical kinetics is that A will decay exponentially while B will rise and fall until no B is left but only the product P. This forms only from B. The product can be calculated by mass balance since $\ce{A0 + B0 =A + B + P}$ at all times. If you have real data to fit then it is still routine to do nonlinear optimisation with numerical solution of equations. – porphyrin Jun 26 '18 at 16:02

The following code is in Python and uses numpy, scipy and matplotlib that are usually packaged with python. The code will solve the equations $\mathrm{A} \underset{k_3}{\stackrel{k_1}\leftrightharpoons} 2\mathrm{B} \stackrel{k_2} \longrightarrow \mathrm{ P}$ but you can set $k_3$ to zero if you want. (its called km1 in the code). I've used log time spaced points so that a wide time range can be seen. The x axis is log in the plots.

I have used a built in solver, odeint(), buy you could write your own based on the Euler or Runge -Kutta methods.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def dX_dt(X, t, const):     # returns rates dEf/dt, only used in odeint()

k1, km1, k2 = const     # unpack constants needed for calc'n
A = X[0]
B = X[1]
dAdt=    -k1*A  +  km1*B*B                 #  rate equations here
dBdt=     k1*A  - (km1+k2)*B*B

f = np.array( [(dAdt), (dBdt)] )           # evaluate equations;
return f

C0 = np.zeros(2,dtype=float)
C0[0] = 1.0e-3    # C0= [ A(0), B(0)   ]      # initial concentrations
numt = 200
t = np.logspace( -9,-4,  numt )    # time as initial power, final power, num points
k1  = 1e7  # 1st order rate const s^(-1)
km1 = 1e6  # 2nd dm^3mol^(-1)s^(-1)
k2  = 1e9  # 2nd  ....
const = [ k1, km1, k2 ]                            # rate constants for calc'n

soln = odeint(dX_dt, C0, t, args = (const,)  )     # solve equations
Aval, Bval = soln.T     # extract data

plt.plot(t,Aval,color='red',  label='A')
plt.plot(t,Bval,color='green',label='B')
plt.plot(t,C0[0]-Aval-Bval,color='blue',label='P')
plt.xscale('log')
plt.legend()

plt.tight_layout()
plt.show()


The figure shows the populations vs (log) time

• I'm not sure what to do with your answer because I took your advice - numerical solution but I don't use python so I cannot comment on your code. I'm going to accept your answer because you suggested numerical solution. – 7E10FC9A Jun 28 '18 at 0:14
• You will find that it is quite easy to convert it into your preferred language. Most now do have the equivalent of the odeint integrator used above. Alternatively download python (along with numpy, matplotlib & scipy ) as all are free to use. – porphyrin Jul 1 '18 at 7:47