I prefer to have the singular behaviour near $0$ rather than at $\frac{\pi}{2}$, so let's make the substitution $\varphi = \frac{\pi}{2} - \theta$. We obtain
$$K(k) = \int_0^{\frac{\pi}{2}} \frac{d\varphi}{\sqrt{1 - k^2\cos^2\varphi}}.$$
Split the integral at $\frac{\pi}{4}$. The part
$$\int_{\frac{\pi}{4}}^{\frac{\pi}{2}} \frac{d\varphi}{\sqrt{1 - k^2\cos^2\varphi}}$$
remains harmless as $k \to 1$ and tends to
$$\int_{\frac{\pi}{4}}^{\frac{\pi}{2}} \frac{d\varphi}{\sin \varphi}.$$
For the other part, write $1 = \sin^2 \varphi + \cos^2 \varphi$ to obtain $1 - k^2\cos^2\varphi = \sin^2 \varphi + (1-k^2)\cos^2 \varphi$. Let $\varepsilon = \sqrt{1-k^2}$. Then
\begin{align}
\int_0^{\frac{\pi}{4}} \frac{d\varphi}{\sqrt{1 - k^2\cos^2 \varphi}} &= \int_0^{\frac{\pi}{4}} \frac{\cos^2 \varphi + \sin^2 \varphi}{\sqrt{\cos^2 \varphi + \sin^2 \varphi}\cdot \sqrt{\sin^2 \varphi + \varepsilon^2 \cos^2\varphi}}\,d\varphi\\
&= \int_0^{\frac{\pi}{4}} \frac{1 + \tan^2\varphi}{\sqrt{1 + \tan^2 \varphi}\cdot \sqrt{\varepsilon^2 + \tan^2 \varphi}}\,d\varphi \tag{$t = \tan \varphi$}\\
&= \int_0^1 \frac{dt}{\sqrt{1+t^2}\cdot \sqrt{\varepsilon^2 + t^2}}\\
&= \int_0^1 \frac{dt}{\sqrt{\varepsilon^2 + t^2}} - \int_0^1 \Biggl(1 - \frac{1}{\sqrt{1+t^2}}\Biggr) \frac{dt}{\sqrt{\varepsilon^2 + t^2}}.
\end{align}
Since
$$1 - \frac{1}{\sqrt{1+t^2}} = \frac{\sqrt{1+t^2}-1}{\sqrt{1+t^2}} = \frac{t^2}{\sqrt{1+t^2}\cdot (1 + \sqrt{1+t^2})},$$
the last integral remains bounded and tends to
$$\int_0^1 \frac{t \,dt}{1 + t^2 + \sqrt{1+t^2}}$$
as $\varepsilon \to 0$.
And the substitution $t = \varepsilon u$ gives us
\begin{align}
\int_0^1 \frac{dt}{\sqrt{\varepsilon^2 + t^2}} &= \int_0^{\frac{1}{\varepsilon}} \frac{du}{\sqrt{1+u^2}}\\
&= \operatorname{Ar sinh} \frac{1}{\varepsilon}\\
&= \log \biggl(\frac{1}{\varepsilon} + \sqrt{1 + \frac{1}{\varepsilon^2}}\biggr)\\
&= \log \frac{1}{\varepsilon} + \log 2 + \log \frac{1 + \sqrt{1+\varepsilon^2}}{2}.
\end{align}
Thus we have
$$K(k) = \log \frac{1}{\sqrt{1-k^2}} + O(1) = \frac{1}{2}\log \frac{1}{1-k} - \frac{1}{2}\log (1+k) + O(1) = \frac{1}{2}\log \frac{1}{1-k} + O(1).$$
In our specific case, with $k = \sin \bigl(\frac{\pi}{2} - \frac{\delta}{2}\bigr) = \cos \frac{\delta}{2}$, we have $\varepsilon = \sqrt{1-k^2} = \sin \frac{\delta}{2} = \frac{\delta}{2} + O(\delta^3)$, so
$$\log \frac{1}{\varepsilon} = \log \frac{2}{\delta} + O(\delta^2)$$
and overall
$$K\bigl(\cos \tfrac{\delta}{2}\bigr) = \log \frac{1}{\delta} + O(1),$$
where the $O(1)$ term is not only bounded, it in fact converges as $\delta \to 0$. We have the relevant terms:
$$K\bigl(\cos \tfrac{\delta}{2}\bigr) = \log \frac{1}{\delta} + 2\log 2 + \int_0^1 \frac{t\,dt}{1+t^2+\sqrt{1+t^2}} + \int_{\frac{\pi}{4}}^{\frac{\pi}{2}} \frac{d\varphi}{\sin \varphi} + o(1).$$
Since $f(x) = x^2+1$, then $f''(x) = 2$. Also $a=0, b=1$, then
$$|\delta| \leq \frac{\max|f''(x)|(b-a)^3}{12n^2} = \frac{2}{12n^2} = \frac{1}{6n^2}.$$
With an error less than $0.01$ we have $\frac{1}{6n^2} \leq 0.01$, so $n \geq \sqrt{\frac{100}{6}} = 4.08$. So have to take $n = 5$ to acheive the required error.
For three subintervales $$|\delta| \leq \frac{1}{54} \simeq 0.019.$$
Best Answer
Basically you're doing everything right, but several improvements can be suggested.
First, the expansion for $K(r)$ near $r = 1$ is $$ K(r) = \log 4 - \frac{1}{2} \log(1-r^2) + \dots $$ The missing $\log 4$ term introduces $O(\Delta r)$ error in your computations.
I suggest splitting $K(r)$ into regular and singular parts as $$ K(r) = -\frac{1}{2} \log(1-r^2) + R(r). $$ Here $R(r)$ is smooth enough to use in numerical quadrature (though high order quadratures would loose orders of convergence due to $R'(1) = -\infty$)
Then do the splitting $$ \int_0^1 r K(r) f(r) dr = \int_0^1 r R(r) f(r) dr - \int_0^1 \frac{r}{2} \log(1-r^2) f(r) dr. $$ Note that the singular part is used for the whole domain $[0, 1]$ and not only for the $[1-\epsilon, 1]$.
The first integral can be evaluated using the trapezoidal rule. Use $R(1) = \log 4$ at the last point and $R(r) = K(r) + \frac{1}{2} \log (1-r^2)$ for the rest.
The second integral needs to be integrated analytically. To do this we need to split the $[0, 1]$ into segments $[r_m, r_{m+1}]$ in which we treat $f(r)$ as a linear function. In other words we consider $f(r)$ as a piecewise-linear function. Using $a = r_m, b = r_{m+1}$ for simplicity we need to evaluate the following integral $$ \int_a^b \frac{r}{2} \log(1 - r^2) \left[ f(a) + \frac{f(b) - f(a)}{b-a} (r - a) \right] dr = \\ = f(a) J_1(a, b) + \frac{f(b) - f(a)}{b-a} J_2(a, b),\\ J_1(a, b) = \int_a^b \frac{r}{2} \log(1 - r^2) dr\\ J_2(a, b) = \int_a^b \frac{r}{2} \log(1 - r^2) (r-a) dr. $$ The exact expressions for $J_1$ and $J_2$ are $$ J_1(a, b) = \frac{1}{4} \left( a^2 - b^2 + (1-a^2)\log(1-a^2) - (1-b^2)\log(1-b^2) \right)\\ J_2(a, b) = \frac{b-a}{36} (5 a^2+5 a b-4 b^2-12) +{} \\ {} + \frac{ (3 a - 3ab^2 + 2b^3) \log (1-b^2) - a(3-a^2) \log (1-a^2) }{12} + {} \\ {} + \frac{\operatorname{arctanh}(b) - \operatorname{arctanh}(a)}{3} $$ I've obtained them using Wolfram Mathematica with some simplifications by myself.
Below is the implementation of the method in Python: