# 4. Integrating ODEs using `odeint`

First, import some stuff we'll need. We give some of them nicknames to save typing later


```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
```

The function `odeint`, which is found in the `integrate` submodule of the `scipy` package, integrates ordinary differential equations of the form
$$
\newcommand{\yy}{\mathbf{y}}
\newcommand{\ff}{\mathbf{f}}
\frac{\mathrm{d} \yy}{\mathrm{d}t} = \ff(\yy, t),
$$where $\yy$ is a vector in $\mathbb{R}^n$, or
$$
\dot{\yy} =  \ff(\yy, t).
$$
That is, it provides a numerical solution of the *initial value problem* given by this differential equation, together with an initial condition
$$
\yy(t = t_0) = \yy_0,
$$
returning the solution curve (or trajectory) $\yy(t)$ as a function of time $t$, within a certain (finite!) range of values of $t$.

Writing $\yy = (y_1, \ldots, y_n)$, we have a system of the form
\begin{align}
\dot{y_1} &= f_1(y_1, \ldots, y_n, t) \\
\dot{y_2} &= f_2(y_1, \ldots, y_n, t) \\
&\ldots \\
\dot{y_n} &= f_n(y_1, \ldots, y_n, t)
\end{align}

Note that in fact, *any* ordinary (system of) ordinary differential equation(s) can be written in this form by defining new variables to represent higher-order derivatives,  
and incorporating them into the vector $\yy$ of variables.

To use `odeint`, we must define at least the Python function `f` which implements the mathematical function $\ff$. It must take a vector `y` and a time `t`,
and return the new vector `f(y,t)`.

## 1D equations

The simplest case is when the spatial dimension is $n=1$. For example, we could wish to solve the ODE
$
\dot{y} = ay,
$
whose exact solution is, of course, $y(t) = y_0 \exp(a t)$.


```python
a = 2.0
def f(y, t):
    return a*y
times = np.linspace(0,6)
y0 = 1.0
y_output = odeint(f, y0, times)
y_exact = y0*np.exp(a*times)
plt.plot(times, y_output, 'r.', label='integrated')
plt.plot(times, y_exact, label='analytical')
plt.legend()
plt.show()
```

# Example E5-6, Calculating X in a reactor with pressure drop
See page [183-188 in the textbook "Essentials of Chemical Reaction Engineering" by Fogler](.

This time we have two variables to integrate (conversion $X$ and pressure ratio $y$), so our function must accept and return a vector. We integrate over weight of catalyst $W$, not time $t$.
The differential equations are:

$$\frac{\mathrm{d} X}{\mathrm{d} W} = \frac{-r_A'}{F_{A0}}$$
$$\frac{\mathrm{d} y}{\mathrm{d} W} = -\alpha \frac{(1+\epsilon X)}{2y}$$



```python
def dYvecdW(Yvec,W):
    """
    Returns the right hand side of the differential equation dY/dW
    where Y is a vector containing two variables, (X, y).
    It thus returns a vector of two variables (dX/dW, dy/dW).
    
    X is conversion (what fraction of the limiting reactant has reacted)
    y is pressure ratio. (P/P_0)
    
    The "vector" Yvec could be an array, list, or tuple. Anything "array-like".
    """
    X, y = Yvec # unpack values from Y vector into two variables
    # Explicit equations
    epsilon = -0.15
    kprime = 0.0074
    Fa0 = 0.1362
    alpha = 0.0367
    raprime = -1*kprime*(1-X)/(1+epsilon*X) * y
    f = (1+epsilon*X)/y # ratio of v / v0
    # Differential equations
    dXdW = -raprime / Fa0
    dydW = -1*alpha * (1 + epsilon * X)/ (2*y) 
    return (dXdW, dydW)  # pack up three differential equations into the dYvecdV vector
```

Set up the steps you want to integrate over. In this case it's Weight not Time


```python
W_output = np.linspace(0, 27., num=50) # Range for the weight of catalyst of the reactor
```

Set up the initial conditions, as a vector of the two values


```python
Y0 = [0., 1.] # Initial values for the dependent variables i.e X and y
```

Perform the integration, and extract the results


```python
Y_result = odeint(dYvecdW, Y0, W_output) # perform the integration
X, y = Y_result.T # extract columns by transposing matrix and splitting into two variables
```

Now we can plot the results


```python
epsilon = -0.15
f = (1+epsilon*X)/y # volumetric flow rate ratio

plt.plot(W_output, X, label='X (conversion)')
plt.plot(W_output, y, label='y ($P/P_0$)')
plt.plot(W_output, f, label='f ($v/v_0$)')
plt.legend(loc='upper left')
plt.xlabel('Weight of catalyst (kg)')

plt.show()
```

# Integrating other variables
Sometimes it's convenient to integrate things like $N_{tot}$ or Volume as well as the variables that you think you care most about. Just add them to your vector.


```python

```

d92395987f47218e23849124841db5a11663c212
ipynb
Jupyter Notebook
4-ODE-integration.ipynb
CHME5137/helpful-examples
b188abd739b6d30d1c326043350e565df998ed7d
MIT
2020-08-23T23:17:59.000Z
2020-08-23T23:17:59.000Z
4-ODE-integration.ipynb
CHME5137/helpful-examples
b188abd739b6d30d1c326043350e565df998ed7d
MIT
4-ODE-integration.ipynb
CHME5137/helpful-examples
b188abd739b6d30d1c326043350e565df998ed7d
MIT
2016-10-17T19:30:26.000Z
2020-08-23T23:18:00.000Z
Qwen/Qwen-72B
1. YES
2. YES