Separating the Variables

A differential equation is an equation in the form:

\[ \frac{dy}{dx} = f(x, y) \]

where \(f(x, y)\) is a function of \(x\) and \(y\). The method of separating the variables works by getting our differential equation into the form:

\[ f_1(y)dy = f_2(x)dx \]

So we have all the ys on the left hand side of the equation and all of the xs on the right hand side, then we can integrate as below:

\[ \int f_1(y) \;dy = \int f_2(x) \;dx \]

---

Achieving the rearrangement of an equation for that all the ys and xs are on opposite sides, is hard to acheive with Python. However, sympy can still be used to find the individual integrals.

from sympy import symbols, integrate

x = symbols('x')

integrate(4 * x ** 3)

\[\displaystyle x^{4}\]

---

Example

If \(\frac{dy}{dx} = \frac{\sin(-x) + e^{4x}}{y^2}\), then use the method of separation of variables to find an expression for \(y\) in terms of \(x\).

Solution: So first we need to get the differential equation into the correect form, with all the xs and y on the right and left sides respectively. So we have:

\[ y^2\;dy = \left(\sin(-x) + e^{4x}\right)\;dx \]

We now integrate both sides, the left with respect to \(y\) and the right with respect to \(x\):

\[\begin{split} \begin{aligned} \int y^2\;dy & = \int \left(\sin(-x) + e^{4x}\right)\;dx \\ \int y^2\;dy & = \int \left(\sin(-x) + \int e^{4x}\right)\;dx \\ \frac{y^3}{3} + C_y & = \cos(-x) + \frac{e^{4x}}{4} + C_x \\ y^3 + 3C_y & = \frac{3}{4}e^{4x} + 3\cos(-x) + 3C_x \\ y & = \left(\frac{3}{4}e^{4x} + 3\cos(-x) + C'\right)^{\frac{1}{3}} \end{aligned} \end{split}\]

where \(C\ = 3C_x - 3C_y\).

---

Again, we start from the reorganised form.

from sympy import sin, exp, Eq, solve

y = symbols('y')

y_side = integrate(y ** 2)
x_side = integrate(sin(-1 * x) + exp(4 * x))
solutions = solve(Eq(y_side, x_side), y)
solutions

[-(3*exp(4*x)/4 + 3*cos(x))**(1/3)/2 - sqrt(3)*I*(3*exp(4*x)/4 + 3*cos(x))**(1/3)/2,
 -(3*exp(4*x)/4 + 3*cos(x))**(1/3)/2 + sqrt(3)*I*(3*exp(4*x)/4 + 3*cos(x))**(1/3)/2,
 (3*exp(4*x)/4 + 3*cos(x))**(1/3)]

There are three possible solutions, accounting for complex results. The final one matches our result.

solutions[-1]

\[\displaystyle \sqrt[3]{\frac{3 e^{4 x}}{4} + 3 \cos{\left(x \right)}}\]

---

Integrated Rate Reaction

In a reaction mixture, at a fixed temperature, the concentration of a reactant, \([A]\), varies with time, \(t\), according to the differential equation:

\[ \frac{d[A]}{dt} = -2k[A]^2 \]

Integrate the equation with the boundary limits of when \(t=0\) and \([A] = [A]_0\) to get the following equation:

\[ \frac{1}{[A]} - \frac{1}{[A]_0} = 2k \]

Solution: First we need to put all the \([A]\)s on the left and all the \(t\)s on the right to give:

\[ \frac{1}{[A]^2} \;d[A] = -2k\;dt \]

We can now apply the integral to both sides:

\[\begin{split} \begin{aligned} \int [A]^{-2} \;d[A] & = \int -2k\;dt \\ -[A]^{-1} & = -2kt + C \end{aligned} \end{split}\]

Now we can substitute in our limits to find \(C\). When \(t=0\) and \([A] = [A]_0\), we have:

\[ -[A]_0^{-1} = -2k \times 0 + C \Rightarrow C = -[A]_0^{-1} \]

Substituting this \(C\) back into our equation, obtained by the integral gives:

\[\begin{split} \begin{aligned} -[A]^{-1} & = -2kt - [A]_0^{-1} \\ \frac{1}{[A]} - \frac{1}{[A]_0} & = 2kt \end{aligned} \end{split}\]

---

Similar to the previous examples, we start from the rearranged equation.

A, t, k, C = symbols('[A] t, k, C')

integrals = Eq(integrate(1 / (A ** 2)), integrate(-2 * k, t) + C)
C = solve(integrals, C)[0].evalf(subs={t: 0})

Eq(integrate(1 / (A ** 2)), integrate(-2 * k, t) + C)

\[\displaystyle - \frac{1}{[A]} = - 2 k t - \frac{1}{[A]}\]

---

Cotrell Equation

In electrochemistry, the Cottrell equation,

\[ I = nFAc\sqrt{\frac{D}{\pi}}t_{-\frac{1}{2}}, \]

describes the current, \(I\), at a time, \(t\), after an electrode is immersed in solution. Current is defined as the rate of change of charge, \(Q\), so:

\[ I = \frac{dQ}{dt}. \]

Using this, find an expression for charge.

Solution: Since \(I = \frac{dQ}{dt}\), we can substitute this into the equation involving \(t\) to form a differential equation:

\[\begin{split} \begin{aligned} \frac{dQ}{dt} & = nFAc\sqrt{\frac{D}{\pi}}t^{-\frac{1}{2}} \\ dQ & = nFAc\sqrt{\frac{D}{\pi}}t^{-\frac{1}{2}}\;dt \\ \int dQ & = \int nFAc\sqrt{\frac{D}{\pi}}t^{-\frac{1}{2}}\;dt \\ \int dQ & = nFAc\sqrt{\frac{D}{\pi}}\int t^{-\frac{1}{2}}\;dt \\ Q & = nFAc\sqrt{\frac{D}{\pi}}\times \frac{t^{-\frac{1}{2} + 1}}{-\frac{1}{2} + 1} + C \\ Q & = nFAc\sqrt{\frac{D}{\pi}}2t^{\frac{1}{2}} + C \end{aligned} \end{split}\]

---

The above is just the integral of the Cottrell equation with respect to \(t\).

from sympy import sqrt

n, F, A, c, D, pi, t = symbols('n F A c D \pi t')

integrate(n * F * A * c * sqrt(D / pi) * t ** -0.5, t)

\[\displaystyle 2.0 A F c n t^{0.5} \sqrt{\frac{D}{\pi}}\]

---

Note

This is one form of the Lennard-Jones potential, with \(a_0\) being the equilibrium separation and \(\varepsilon\) being the energy needed to move the particles apart to infinity.

Lennard-Jones Potential

The force between two particles is modelled to be:

\[ F = \frac{12\varepsilon}{a_0}\left[\left(\frac{a_0}{r}\right)^{13} - \left(\frac{a_0}{r}\right)^7\right] \]

Given that force is the negative derivative of potential, i.e., \(F = -\frac{d}{dr}U\), and at \(r=a_0\), \(U=-\varepsilon\), calculate the potential between the two particles.

Solution: This is a problem of separating variables, so first we set up the problem:

\[\begin{split} \begin{aligned} -\frac{d}{dr}U & = \frac{12\varepsilon}{a_0}\left[a_0^{13}r^{-13} - a_0^{7}r^{-7}\right] \\ -\int dU & = \int \frac{12\varepsilon}{a_0}\left[a_0^{13}r^{-13} - a_0^{7}r^{-7}\right]\;dr \\ -U & = \frac{12\varepsilon}{a_0}\int\left[a_0^{13}r^{-13} - a_0^{7}r^{-7}\right]\;dr \\ -U & = \frac{12\varepsilon}{a_0}\left(-\frac{a_0^{13}}{12}r^{-12} + \frac{a_0^{7}}{6}r^{-6}\right) + C \\ U & = \varepsilon\left(a_0^{12}r^{-12} - 2a_0^6r^{-6}\right) + C \end{aligned} \end{split}\]

Now we find \(C\) by substitution:

\[\begin{split} \begin{aligned} C & = -\varepsilon - \varepsilon\left(a_0^{12}a_0^{-12} - 2a_0^6a_0^{-6}\right) \\ & = -\varepsilon - \varepsilon(1-2) \\ & = 0 \end{aligned} \end{split}\]

Back into our equation, we get \(U = \varepsilon\left(a_0^{12}r^{-12} - 2a_0^6r^{-6}\right)\), which is the Lennard-Jones potential.

---

sympy is particularly helpful for the integration of complex equations, as in this example.

a0, r, epsilon = symbols('a0 r \epsilon')

integral = -1 * integrate(12 * epsilon / a0 * ((a0 / r) ** 13 - (a0 / r) ** 7), r)
integral

\[\displaystyle \frac{\epsilon a_{0}^{12} - 2 \epsilon a_{0}^{6} r^{6}}{r^{12}}\]

And we can solve for \(C\).

C = symbols('C')

solve(Eq(-epsilon, integral.evalf(subs={r: a0}) + C), C)

[0]

---