Solving differential equations

Euler's method and the midpoint method for solving ODEs | View | Download | |

Numerical solution of the Lotka-Volterra predator-prey equation | View | Download | |

The Runge-Kutta and Runge-Kutta-Fehlberg methods | View | Download |

`solve_ivp()`

to solve the differential equation $\dot{x}=x^3$
for $0\leq t\leq 0.1$ with $x=1$ at $t=0$, and then plot the result. Your code
should look something like this:
```
t = np.linspace(0, 0.1, 1000)
sol = solve_ivp(..., t_eval=t)
plt.plot(t, ...) # option 1
plt.plot(sol.t, ...) # option 2
```

`sol.t`

will be the same as `t`

.
Now do the same thing but for $0\leq t\leq 1$. You should find that option 1 gives an error but option 2 still works. If you look at the plot for option 2, you will see that the solution goes to infinity at $t=0.5$ and the solver gives up at that point. Thus

`sol.t`

only contains times before $0.5$, but
`t`

has times after $0.5$ which have no corresponding $x$ values,
which is why option 1 gives an error.
Now define a function

`cutoff(a)`

as follows.
- You can assume that $a$ is a real number with $a\geq 1$.
- Ask
`solve_ivp()`

to solve $\dot{x}=x^3$ for $0\leq t\leq 1$ with $x=a$ at $t=0$. Make sure that your original array $t$ of times has length at least $10^4$. - The solution will always go to infinity at some time $t_0\lt 1$.
The last entry in the array
`sol.t`

returned by`solve_ivp()`

will be slightly less than $t_0$. This last entry should be the return value of your function`cutoff(a)`

.

There is in fact a simple formula for

`cutoff(a)`

, which you should
try to discover. You can calculate `cutoff(a)`

for $a=1,2,3,4,\dotsc$
and try to spot the pattern, or you can plot `cutoff(a)`

against $a$
for a range of values of $a$. If the pattern is not apparent, you can plot
`2*cutoff(a)`

or `1/cutoff(a)`

or `np.sqrt(cutoff(a))`

or various other things that might make the pattern clearer.
Alternatively, you can try to solve the differential equation analytically and deduce a formula for the cutoff from that. You are allowed to search the web or ask Wolfram Alpha or an AI assistant for help with this.

You should enter your function

`cutoff(a)`

into the online test system.
If you have found a formula for `cutoff(a)`

, you can enter a function
that just uses the formula and does not call `solve_ivp()`

. If you have
not found a formula, you should enter a function that calls `solve_ivp()`

as described above.
0 | $CH_4$ | Methane |

1 | $O_2$ | Oxygen |

2 | $CO$ | Carbon monoxide |

3 | $H_2O$ | Water |

4 | $CO_2$ | Carbon dioxide |

We will need the following numbers, vectors and matrices. The units are ones typically used in combustion chemistry, based on calories (for energy), centimeters (for length), seconds (for time), Kelvin (for temperature) and moles (for amount of substance). You can click here to see all the numbers in a form that you can copy and paste into your code.

R | 1.987 | The universal gas constant, in $\text{cal}/\text{mol K}$ |

$P_0$ | $2.422\times 10^{-2}$ | Atmospheric pressure, in $\text{cal}/\text{cm}^3$ |

A | $[1.59\times 10^{13}, 3.98\times 10^8, 6.16\times 10^{13}]$ | Arrhenius pre-exponential factors |

$\beta$ | $[0, 0, -0.97]$ | Temperature exponents |

$E$ | $[4.78\times 10^4, 1.00\times 10^4, 7.84\times 10^4]$ | Activation energies, in $\text{cal}/\text{mol}$ |

$N$ | $\begin{bmatrix} -1& 0& 0 \\ -1.5 &-0.5& 0.5\\ 1&-1& 1\\ 2& 0& 0\\ 0& 1&-1 \end{bmatrix}$ | Stoichiometric matrix |

$U$ | $\begin{bmatrix} 1 & 0 & 1 & 0 & 1 \\ 0 & 2 & 1 & 1 & 2 \\ 4 & 0 & 0 & 2 & 0 \end{bmatrix} $ | Numbers of atoms in each compound. |

The modified Westbrook-Dryer model for methane combustion is as follows. Everything happens at a temperature $T$ (in Kelvin) and at atmospheric pressure $P_0$. We first need to calculate rate constants $k_i=A_i T^{\beta_i} \exp(-E_i/(RT))$ for $i=0,1,2$ (where the numbers $A_i, \beta_i, E_i, R$ are given above). As initial data we are also given the molar fractions of the five compounds at time $t=0$. These are numbers $p_0,\dotsc,p_4\geq 0$ with $p_0+\dotsb+p_4=1$. We convert these to molarities $x_i$ by the rule $x_i=p_iP_0/(RT)$. The numbers $x_i$ then change according to a differential equation. To describe this, we first put $y_i=\max(x_i, 10^{-9})$. (This is only necessary to avoid some problems arising from inaccuracies in the solution process.) We then put \begin{align*} r_0 &= k_0 y_0^{0.70} y_1^{0.80} \\ r_1 &= k_1 y_1^{0.25} y_2 y_3^{0.50} \\ r_2 &= k_2 y_1^{-0.25} y_3^{0.50} y_4. \end{align*} This gives a vector $r$ of length $3$, so we can multiply by the matrix $N$ to get a vector $N r$ of length $5$. The differential equation is then given by $\dot{x}=N r$. After solving this, we can convert back to molar fractions $p_i=x_i/(x_0+\dotsb+x_4)$.

- Write a one-line function
`rate_constants(T)`

which calculates the vector $k=(k_0, k_1, k_2)$ of rate constants for a given temperature $T$. -
Write a one-line function
`fraction_to_molarity(p, T, P=P0)`

that converts the vector $p$ of mole fractions to the vector $x$ of molarities. This depends on the temperature $T$ and the pressure $P$, but the pressure should be taken to be atmospheric pressure by default. -
Write a one-line function
`molarity_to_fraction(x)`

that converts the vector $x$ of molarities back to the vector $p$ of mole fractions. -
Write a function
`x_dot(t, x, k)`

that calculates the derivative $\dot{x}$ according to the rule described above. The first argument $t$ is the time, which is required by`scipy.integrate.solve_ivp`

even though it is not used in this case. The second argument $x$ is the vector of molarities. The third argument $k$ is the vector of rate constants, which will be calculated by the function`rate_constants()`

. You should tell`solve_ivp()`

to use the`BDF`

method, with an absolute tolerance of $10^{-9}$. Note that the function`x_dot()`

has an additional argument of $k$ as well as the usual arguments $t$ and $x$. You may need to search for information about how to deal with this in`solve_ivp()`

. You will find that`solve_ivp()`

requires a tuple of additional arguments; remember that a tuple with one element must be written with notation like $(k,)$ rather than $(k)$ or $k$. -
Write a function
`solve_combustion(T, p0, t_max, n)`

that carries out the steps described above: convert the initial data $p0$ to molarities $x0$, calculate the rate constants $k$, solve the differential equation for $x$ using`scipy.integrate.solve_ivp`

, convert the result back to molar fractions. The function`solve_combustion()`

should call the other functions that you have written. It should return a pair $(t, p)$, where $t$ is an array of $n$ equally spaced times from $0$ to $t_\text{max}$ and $p$ is an array of shape $(n,5)$ whose $i$th row is the vector of molar fractions at time $t_i$. -
Write a function
`show_combustion(T, p0, t_max, n)`

that calls`solve_combustion()`

and then plots the molar fractions against time. The plot should have a legend showing the chemical formulae of the compounds, and labels on the axes. You should use a logarithmic scale on the time axis. The number $n$ will need to be large. If you take $n=1000$, you will see that the curve looks jagged; you should increase $n$ until the curve becomes smooth. - Test your function with $T=1500$, $p0=[0.1, 0.25, 0.01, 0.01, 0.63]$ and $t_\text{max}=0.1$. (These initial molar fractions correspond to a situation in which the nitrogen has been removed from the air and replaced by carbon dioxide, which Andersen et al. were studying for applications to carbon capture and storage.)
- The matrix $U$ is described as "numbers of atoms in each compound". That is not a very detailed explanation, but there is only one reasonable interpretation that is consistent with the other information that you have in front of you. What is it? Show that $U\dot{x}=0$, and explain the significance of this equation.