I am looking for some advice on how to get a simple approximation to an integral that appears in the classical SIR epidemiological model.
One question, that I keep pondering given recent news items, is given different configurations of infectiousness and population size, how long will it take from an initial state to reach the fabled herd immunity . e.g. if all but a small age group are perfectly vaccinated, how long will it take to reach herd immunity in that population with a given transmission rate? Does a population with greater $R_{0}$ always take longer to reach herd? What if the population is smaller? etc.
Following the terminology (and solution strategy) in the wiki page Compartmental models in epidemiology: The SIR model without vital dynamics, then given an initial state of an epidemic as $(S(0),I(0), R(0))=(N(1-\delta),N\delta, 0)$ the time to progress to the herd immunity state where $S(\tau_{herd})/N=1/R_{0}$ (leading to $dI/dt=0$) is given by the integral: $$ \tau_{herd}=\int_{0}^{\ln(R_{0}(1-\delta))/R_{0}} \frac{\mathrm{d}x}{1-x-(1-\delta)\exp(-R_{0}x)} $$ In the limit $\delta\rightarrow 0$ the integral diverges, as is intuitive. However, trying to find an approximate form valid for small $\delta$ doesn't seem obvious to me, due to the mixture of linear and exponential forms in the denominator. Integration by parts does not seem fruitful and methods such as Laplace etc. rely on exponential forms: $\lim_{\lambda\rightarrow \infty} \int \exp(\lambda f(x)) \mathrm{d}x$. It seems hard to force the $\delta$ parameter dependence into such an exponential. A clever substitutions to get the above expression into a more tractable form seems not obvious.
The only approach I could think of was to expand $\exp(-R_{0}x)\simeq 1-R_{0}x +R_{0}^{2}x^{2}/2$ which makes the integral elementary and gives the (awkward) form: $$ \tau_{herd}\simeq \frac{1}{R_{0}(1-\delta)}\frac{1}{ \sqrt{ \frac{ (1+R_{0}(1-\delta))^{2}}{R_{0}^{4}(1-\delta)^{2}} - \frac{2}{R_{0}^{2}}\frac{\delta}{1-\delta} } } \ln \left[ \frac { 1+ \frac{ \ln[R_{0}(1-\delta)]/R_{0} } { \frac{1+R_{0}(1-\delta)}{R_{0}^{2}(1-\delta)} - \sqrt{ \frac{ (1+R_{0}(1-\delta))^{2}}{R_{0}^{4}(1-\delta)^{2}} - \frac{2}{R_{0}^{2}}\frac{\delta}{1-\delta} } } } { 1+ \frac{ \ln[R_{0}(1-\delta)]/R_{0} } { \frac{1+R_{0}(1-\delta)}{R_{0}^{2}(1-\delta)} + \sqrt{ \frac{ (1+R_{0}(1-\delta))^{2}}{R_{0}^{4}(1-\delta)^{2}} - \frac{2}{R_{0}^{2}}\frac{\delta}{1-\delta} } } } \right] $$ (I gather this was the approach in the original 1920s literature). This expression matches numerical integration in an OK-ish manner for small $R_{0}$ and $\delta$ but I feel there is a better approximation here that is valid for larger ranges of $R_{0}$.
Is there a more systematic theory available to find expansions for forms such as $\tau_{herd}(\delta)$?
Thanks!
I think that your approach could be improved using first $x=\frac t {R}$ $$I=\int\frac{dx}{1-x-(1-\delta ) e^{-R x}}=\int\frac{dt}{R-t-R(1-\delta )\, e^{-t}}$$ and, instead of Taylor series for the exponential, use the simplest Padé approximant $$e^{-t}=\frac{2-t}{2+t}+O(t^3)$$ $$I\sim J=-\int \frac {t+2} {t^2+(2-2R+R\delta)t-2R\delta}\,dt$$ Let $$a=\frac{1}{2} \left(2(R-1)-\delta R-\sqrt{((\delta -2) R+2)^2+8 \delta R}\right)$$ $$b=\frac{1}{2} \left(2(R-1)-\delta R+\sqrt{((\delta -2) R+2)^2+8 \delta R}\right)$$ being the roots of the quadratic $$I=\frac{(b+2) \log (t-b)-(a+2) \log (t-a)}{a-b}$$ $$-\int_0^p \frac {t+2} {t^2+(2-2R+R\delta)t-2R\delta}\,dt=\frac{(a+2) \log \left(\frac{a}{a-p}\right)-(b+2) \log\left(\frac{b}{b-p}\right)}{a-b}$$ with $p=\log\big[R(1-\delta)\big]$
Without any idea about what could be $R$, I tried with $R=10$ and a few values of $\delta$. $$\left( \begin{array}{ccc} \delta & \text{approximation} &\text{exact} \\ 0.01 & 0.743557 & 0.774830 \\ 0.02 & 0.665511 & 0.695908 \\ 0.03 & 0.619463 & 0.649007 \\ 0.04 & 0.586514 & 0.615226 \\ 0.05 & 0.560737 & 0.588636 \\ 0.06 & 0.539491 & 0.566595 \\ 0.07 & 0.521367 & 0.547693 \\ 0.08 & 0.505523 & 0.531087 \\ 0.09 & 0.491416 & 0.516233 \\ 0.10 & 0.478676 & 0.502760 \end{array} \right)$$
Edit
If you want to make it better, use $$e^{-t}=\frac{t^2-6 t+12}{t^2+6 t+12}+O(t^5)$$ which would give $$I\sim J=-\int \frac{t^2+6 t+12}{t^3 +(6-\delta R)t^2+6 ((\delta -2) R+2)t-12 \delta R}\,dt$$
Let $(a,b,c)$ be the roots of the cubic. This would give $$J=\frac{\left(a^2+6 a+12\right) \log (t-a)}{(a-b) (a-c)}+\frac{\left(b^2+6 b+12\right) \log (t-b)}{(b-a) (b-c)}+\frac{\left(c^2+6 c+12\right) \log (t-c)}{(c-a) (c-b)}$$ Just apply the bounds as before.
Repeating the same calculations $$\left( \begin{array}{ccc} \delta & \text{approximation} &\text{exact} \\ 0.01 & 0.776267 & 0.774830 \\ 0.02 & 0.697296 & 0.695908 \\ 0.03 & 0.650348 & 0.649007 \\ 0.04 & 0.616520 & 0.615226 \\ 0.05 & 0.589884 & 0.588636 \\ 0.06 & 0.567799 & 0.566595 \\ 0.07 & 0.548853 & 0.547693 \\ 0.08 & 0.532204 & 0.531087 \\ 0.09 & 0.517309 & 0.516233 \\ 0.10 & 0.503795 & 0.502760 \end{array} \right)$$ which looks to be much better.
We could continue for ever using better and better $[n,n]$ Padé approximants for $e^{-t}$. This will make $$J=-\int \frac{P_n(t)}{\prod_{k=1}^{n+1} (t-r_i)}\,dt$$ leading to $(n+1)$ logarithms.