# Series Approximation

The basic idea of series approximation is to express the delta quantity $\Delta_n$ in form of an infinite series that can be used to quickly approximate $\Delta_n$ by evaluating a limited number of initial terms. If the approximation turns out to be sufficiently accurate for large $n$, a considerable speed up can be achieved, because it is then no longer necessary to compute $\Delta_n$ by iteration.

In the example of K. I. Martin the series is composed out of three initial terms and a residual term. Writing $\delta$ for $\Delta_0$ Martin's equation looks like this: 
\begin{align*}
\Delta_n = A_n \delta + B_n \delta^2 + C_n \delta^3 + \text{o}(\delta^4)
\end{align*}
In generalized terms, this series can be written down like so: 
\begin{align*}
\Delta_n = a_{1,n} \delta + a_{2,n} \delta^2 + a_{3,n} \delta^3 + \ldots + a_{k,n} \delta^k + \text{o}(\delta^{k+1}) \tag{1} 
\end{align*}
If the residual term $\text{o}(\delta^{k+1})$ is sufficiently small, this formula can be utilized to approximate $\Delta_n$:
\begin{align*}
\Delta_n \approx a_{1,n} \delta + a_{2,n} \delta^2 + a_{3,n} \delta^3 + \ldots + a_{k,n} \delta^k \tag{2}
\end{align*}
In practice, this approach works very well in many cases: The residual term often remains small for substantial time until it exceeds a given error threshold. In such a case, we can speed up the computation of a delta orbit by skipping a large initial segment. This means that instead of computing 
\begin{align*}
\Delta_0, \Delta_1, \Delta_2, \ldots
\end{align*}
the series 
\begin{align*}
\Delta_m, \Delta_{m+1}, \Delta_{m+2}, \ldots \tag{3}
\end{align*}
is computed where $m$ denotes the largest value of $n$ for which $(2)$ is sufficiently accurate. In many cases, a large enough initial piece can be skipped that only a small fraction of the iterations need be performed that would have been necessary without series approximation.

For applying series approximation in practice, we need determine up to which value equation ($2$) is sufficiently accurate. The common solution is to use a certain set of *probe points*. For each of these points, the delta orbit is generated by iteration and compared to the values obtained by series approximation. If the relative error exceeds a certain threshold, calculation is aborted. The test point which exceeds the error threshold first, defines the value of $m$ in equation $(3)$.

We are missing one last piece of the puzzle. So far, we have dealt with $(1)$ as an abstract formula, assuming that the coefficients $a_{1,n}, a_{2,n}, a_{3,n} \ldots$ have already been calculated. We still need to find a way to calculate these coefficients. $n = 0$ is the easy case. Substituting $n$ by $0$ in $(1)$ reveals:
\begin{align*}
a_{1,n} &= 1 \nonumber \\
a_{j,n} &= 0 \;\;\text{for}\;\; j \geq 2 \nonumber
\end{align*}
For $n \geq 1$ we find the answer by applying the delta iteration scheme to:
\begin{align*}
\Delta_n = \sum_{j = 1}^{k} a_{j,n} \delta^j
\end{align*}
This result in:
\begin{align*}
\sum_{j = 1}^{k} a_{j,n+1} \delta^j 
&= \Delta_{n+1} \\
&= 2 x_n \Delta_n +  \Delta_n^2 + \Delta_0 \\
&= 2 x_n \Bigl(\; \sum_{j = 1}^{k} a_{j,n} \delta^j\; \Bigr) + \Bigl(\; \sum_{j = 1}^{k} a_{j,n} \delta^j \; \Bigr)^2 + \delta \\
&= \sum_{j = 1}^{k} 2 x_n a_{j,n} \delta^j + \sum_{i = 1}^{k}\sum_{j = 1}^{k} a_{i,n} a_{j,n} \delta^{i+j} + \delta 
\end{align*}
Comparing of coefficients reveals the relationships we are looking for: 
\begin{align*}
a_{1,n+1} &=  2 x_n a_{1,n} + 1 \nonumber \\
a_{j,n+1} &=  2 x_n a_{j,n} + \sum_{i = 1}^{k} a_{i,n} a_{k - i,n} \;\;\text{for}\;\; j \geq 2 \nonumber
\end{align*}