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
.
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.
cutoff(a)
as follows.
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$.
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)
.
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.
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 |
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. |
rate_constants(T)
which calculates
the vector $k=(k_0, k_1, k_2)$ of rate constants for a given temperature
$T$.
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.
molarity_to_fraction(x)
that
converts the vector $x$ of molarities back to the vector $p$ of
mole fractions.
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$.
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$.
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.