Recall, a rate equation looks like this: $$\frac{dy}{dt} = f(y,t).$$
The exact solution is $$C(t) = C_0e^{-t/\tau},$$ where $C_0$ is the initial concentration.
odeint
¶odeint
from scipy.integrate import odeint
f(y,t)
t
C0
odeint
as C=odeint(f, C0, times)
from scipy.integrate import odeint
Solve this problem with $\tau=1$, $C_0=1$, and t_end=5.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint
#------ Step 1, define the function
#------ Step 2, set times to get the solution at
#------ Step 3, set initial condition
#------ Step 4, solve the problem
#------ Step 5 plot the answer
#------ Step 6: does the answer make sense? Analyse it.
Solve the rate equation for $v(t)$ at discrete points $t$ and plot the result:
$$\frac{dv}{dt} = g - cv^2.$$
What if we have multiple rate equations in multiple variables?
If we have multiple rate equations in multiple unknowns, we solve the problem in the same way, but using vector notation.
fsolve
.Solve for the velocity of a falling object: \begin{align} \frac{dv}{dt} &= g - cv^2, \\ \frac{dx}{dt} &= v. \end{align}
Here, g is gravitational acceleration, x is position, v is velocity, and c is a drag coefficient.
We first cast this in the form $$\frac{d\vec{y}}{dt} = \vec{f}(\vec{y},t).$$
This is equivalent to \begin{align} \frac{dy_0}{dt} &= f_0(y_0,y_1,t) \\ \frac{dy_1}{dt} &= f_1(y_0,y_1,t) \end{align}
For our example problem $v\equiv y_0$, $x\equiv y_1$, $g-cv^2=f_0(v,x,t)$, and $v=f_1(v,x,t)$
fsolve
:#------ Step 1, define the function
#------ Step 2, set times to get the solution at
#------ Step 3, set initial condition
#------ Step 4, solve the problem, recover v, x arrays
#------ Step 5 plot the answer
#------ Step 6: does the answer make sense? Analyse it.
Solve the following rate equation for a coal particle in a furnace. $$ \frac{dT}{dt} = \left(\frac{hA}{m c_p}(T_f - T) + \frac{\sigma A}{m c_p}(T_f^4-T^4)\right).$$
The initial particle temperature is $T_0=500$ K.
Use D=100 $\mu$m, tend = 0.05 s.
The following data is given:
Variable | Value | Units |
---|---|---|
$\rho_p$ | 1000 | kg/m$^3$ |
$c_p$ | 1380 | J/kg$\cdot$K |
$k$ | 0.1 | W/m$\cdot$K |
$Nu$ | 2 | -- |
$\sigma$ | 5.67E-8 | W/m$^2$K$^4$ |
$T_f$ | 1500 | K |
Also, the area, mass and Nusselt number $Nu$ are given, respectively, by: $$ A = \pi D^2, $$ $$ m = \frac{\pi}{6}D^3\rho_p,$$ $$ Nu = \frac{hD}{k}.$$
We are performing a chemical reaction as follows. $$\mbox{Rxn 1: }\phantom{xxx} A+B\rightarrow C $$ $$\mbox{Rxn 2: }\phantom{xxx} B+C\rightarrow D $$ Here, symbols $A$, $B$, $C$, $D$ denote species concentrations in mol/L.
The initial concentrations are $A_0=1$, $B_0=1$, $C_0=0$, $D_0=0$. Also, $k_1=1\,L/mol\cdot s$, and $k_2=1.5\,L/mol\cdot s$.
The concentrations obey the following rate equations:
\begin{align} \frac{dA}{dt} &= -k_1AB, \\ \frac{dB}{dt} &= -k_1AB - k_2BC, \\ \frac{dC}{dt} &= k_1AB - k_2BC, \\ \frac{dD}{dt} &= k_2BC. \end{align}Solve for the concentrations of $A$, $B$, $C$, and $D$ as functions of time to $t=3$ s.
Plot the concentrations of A, B, C, and D as functions of time on the same plot. Label the axes as "time (s)" and “concentration (mol/L)”.