MPS224 Home

MPS224 Scientific Computing: Lab 4
Solving differential equations

Instructions

The following notebooks and videos are relevant for this lab:
       
The Duffing oscillator View Download
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

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. 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.

Task 2: Simple harmonic motion and error

(This task is not part of the online test).

Use solve_ivp() to solve the differential equation for Simple Harmonic Motion: $$ \ddot{x} + x = 0 $$ for $0\leq t\leq 20$ with $x(0)=1$, $\dot{x}(0)=0$ at $t=0$, and then plot the result in two ways: $x(t)$ as a function of $t$ and as a phase portrait ($x$ against $\dot{x}$).

Before you can use solve_ivp() you will need to make a first-order reduction (i.e. write the system as a pair of first-order ODEs in two variables, with $y \equiv \dot{x}$).

For simple harmonic motion, the energy $E = \tfrac{1}{2}( x^2 + y^2 )$ should remain constant. However, in your data this may not be precisely so, due to error in the numerical solution. Let's investigate this.

Write a function to compute the energy from a vector $u = [x, y]$. Plot $E - E_0$ against $t$ on the domain $t \in [0, 2000]$ (where $E_0$ is the initial energy, $E_0 = 1/2$ in this case). What do you find?

Now read the help file for solve_ivp(), and experiment with changing the optional arguments atol and rtol to smaller values, and try changing the method. Note the effect of changing these parameters on the 'error' (the difference between $E$ and its initial value $E_0$), and also on how long it takes for solve_ivp() to run.

You should be able to reduce the 'drift' in the energy to less than 1 part in $10^{10}$ over this domain. It may help to read about 'stiff' ODEs and suitable numerical methods for stiff systems.

Task 3: Methane combustion

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)$.