Relating the Nonlinear and Linear operators in the Homotopy Analysis Method The question refers to chapter two of the book Liao, Shijun. Homotopy analysis method in nonlinear differential equations. Beijing: Higher Education Press, 2012. 
Link to book pdf from Chinese .edu site (presumably author related).
In section 2.3.2 of the book the homotopy analysis is being applied to a simple nonlinear differential equation.
Assume a differential equation of the form:
$$\mathcal{N}(x)=\gamma x^{''}+\lambda 
x+\varepsilon x^3=0 \tag{2.35}\label{eq235}$$
where $\mathcal{N}$ is the specified nonlinear operator based on $\lambda$ and $\varepsilon$. Assume $x_0\in S_p$ is an initial guess for the true value of $x$ that solves the differential equation and is contained in some periodic function $S_p$ that also contains the true value of $x$. The book then notes that $\mathcal{N}(x)\in S_p$ also.
Next it creates an auxiliary linear operator $\mathcal{L}$ with the property of 
$\mathcal{L}(0)=0$ and assumes that $\mathcal{L}$ is chosen such that $\mathcal{L}[x-x_0]\in S_p$. It then states that the linear and nonlinear operators have the same period and can therefore be deformed into each other before creating the following homotopy:
$$\mathcal{H}(x:q):=(1-q)\mathcal{L}[x-x_0]-\mathcal{N}(x) \tag{2.37}\label{eq237}$$
This appears to me that the homotopy is being created between the linear and nonlinear operators to convert the nonlinear problem into a linear one. Why does this work when the linear and nonlinear operators are different, and in fact the linear operator is somewhat arbitrary? This appears different to me than the example below where a nonlinear equation is solved. The book just performs the steps but does not explain the why or how of the transition.
The basic explanation of Homotopy applied to nonlinear equations is below:
Considering a nonlinear equation $f(x)=0$ with $x\in[a,b]$
Assume that $x_0$ is an initial guess for the value of $x$ that solves the nonlinear equation.
$f(x)-f(x_0)$ is also continuous on $x\in[a,b]$ and can thus be deformed to $f(x)$
A homotopy can then be constructed:
$$\mathscr{H}(x:q)=\left( 1-q \right) \left[f(x)-f(x_0)\right]+qf(x)=0$$
If $x$ is replaced by homotopy $\tilde{x}(q)$ where $\tilde{x}(q)=(1-q)x_0+qx$
and $\tilde{x}(q)$ is assumed to be analytic at $q=0$. Then the original homotopy becomes:
$\mathscr{H}(\tilde{x}(q):q)=(1-q)[f(\tilde{x}(q))-f(x_0)]+qf(\tilde{x}(q))=0 \tag{2.10}\label{eq210}$
Based on assuming that $\tilde{x}$ is analytic at $q=0$ you assume that a Maclaurin series of $\tilde{x}$ using the variable $q$ will converge.
$$\tilde{x}(q)\sim x_0+\sum_{k=1}^{\infty}{x_kq^k}\tag{2.11}\label{eq211}$$
$$x_n=\frac{1}{k!}\frac{\mathcal{d}^kx}{\mathcal{d}q^k}|_{q=0}=\frac{1}{k!}\mathcal{D}_k[\tilde{x}(q)]
 \tag{2.12}\label{eq212}$$
Where $\mathcal{D}_k$ is the Homotopy derivative of order $k$. and can be obtained by taking the $k^{th}$ derivative of Eq. $(2.10)$ and setting $q=0$
For example:
$$\frac{\partial}{\partial q} \left\{ \left(1-q\right)\left[f(\tilde{x}(q))-f(x_0)\right]+qf(x_0)=0\right\} \tag{2.15}\label{eq215}$$
$$f^{'}(x_0) \frac{\partial\tilde{x}}{\partial q}|_{q=0}=-f(x_0) $$
$$x_1=\mathcal{D}_1=-\frac{f(x_0)}{f^{'}(x_0)}$$
Then set q=1 and approximate the true solution of $f(x)=0$ as:
$$x\approx x_0+x_1+x+2+...=x_0-\frac{f(x_0)}{f^{'}(x_0)}-\frac{[f(x_0)]^2f^{''}(x_0)}{2[f^{'}(x_0)]^3}+... \tag{2.18}\label{eq218}$$
 A: With the choice ${\cal L}[x]=γx''+λx$ you get the homotopy equation
$$
0={\cal H}_q[x]=γx''+λx+qεx^3-(1-q)(γx_0''+λx_0).
$$
Set $x_0=R\sin(\frac{2\pi}{P}t)$ for a $P$ periodic function. But what amplitude $R$ to choose?
To get an initially sensible amplitude consider an infinitesimal parameter  $q\gtrapprox0$ and look for a solution in the form $x=x_0+qx_1+...$. The demand here is that $x_1$ is of a normal size. Inserting the perturbation series we get as equation for $x_1$
$$
γx_1''+λx_1=-(γx_0''+λx_0+εx_0^3)
=γω^2R\cos(ωt)-λR\cos(ωt)-ε\frac{R^3}4[\cos(3ωt)+3\cos(ωt)]
$$
To conclude the task of finding $R$, demand that the amplitude and phase of the base frequency $ω=\frac{2\pi}P$ remains unchanged, that is, the contributions of this frequency to the forcing term add to zero, so that
$$
0=γω^2-λ-ε\frac{3R^2}4\implies R=\sqrt{\frac{4(γω^2-λ)}{3ε}}.
$$
This of course assumes that $λ<γω^2$ or $λP^2<4\pi^2γ$. 
Now go from solution to solution, for instance using multiple shooting where the last iterate is the initial guess for the solver in the next step.
