# MAS2008 Scientific Computing: Lab 4 Solving differential equations

## Instructions

The following notebooks and videos are relevant for this lab:

## Task 1: The equation $\dot{x}=x^3$

Use 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

It does not matter whether you use option 1 or option 2, because 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.

This task (which is not included in the online test) is about the combustion of methane. It loosely follows this paper by Andersen et al.. The following chemical compounds are involved:
 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.