Learning Scientific Programming with Python (2nd edition)
P8.2.9: Modelling the evolution of a match flame
Question P8.2.9
A simple model of the evolution of a match flame considers the flame radius, $y$, to change in time as $$ \frac{\mathrm{d}y}{\mathrm{d}t} = \alpha y^2 - \beta y^3, $$ where $\alpha$ and $\beta$ are some constants relating to the transport of oxygen through the surface of the flame and the rate of its consumption inside it. The flame is initally small, $y(0) \ll \alpha/\beta$, but at some point rapidly grows until it reaches a steady state of constant radius (assuming an unlimited supply of fuel).
Taking $\alpha=\beta=1$, solve this ODE numerically using scipy.integrate.solve_ivp using a suitable integration method over a time interval $(0, 5/y(0))$ for (a) $y(0) = 0.01$, (b) $y(0) = 0.0001$. How many time steps must be taken in each case?
The exact solution may be written as
$$
y(t) = \frac{\alpha}{\beta\left[1+W\left(a\mathrm{e}^{a-\alpha^2t/\beta}\right)\right]},
$$
where $a=\alpha/(\beta y(0))$ and $W(x)$ is the Lambert $W$ function, implemented in SciPy as scipy.special.lambertw. Compare the accuracy of the various numerical solutions with this expression.
