Q: Approximations for a spring pendulum's equations of motion in 2D I'm working on Exercise 24 in Classical Mechanics, 3rd ed by Goldstein, Poole, and Safko. It concerns the spring pendulum and approximations to its equations of motion. I'm trained more in pure mathematics than physics and I'm not as good at physics-style approximations as I'd like to be. I have questions on parts (b) and (c), so I'll quote everything up to part (c):

A spring of rest length $L_a$ (no tension) is connected to a support at one end and has a mass $M$ attached at the other. Neglect the mass of the spring, the dimension of the mass $M$, and assume that the motion is confined to a vertical plane. Also, assume that the spring only stretches without bending but it can swing in the plane.
  
  
*
  
*(a) Using the angular displacement of the mass from the vertical and the length that the spring has stretched from its rest length (hanging with the mass $M$), find Lagrange's equations.
  
*(b) Solve these equations for small stretching and angular displacements.
  
*(c) Solve the equations in part (a) to the next order in both stretching and angular displacement. This part is amenable to hand calculations. Using some reasonable assumptions about the spring constant, the mass, and the rest length, discuss the motion. Is a resonance likely under the assumptions stated in the problem?
  

Part (a): The equations I found are
\begin{align*}
(l+a)\ddot\theta + 2\dot\theta\dot a + g\sin\theta &= 0 \\
\ddot a - (l+a)\dot\theta^2 + g(1-\cos\theta) + \frac{k}{M}a &= 0.
\end{align*}
Here $\theta$ is the angle from the vertical and $a$ is the displacement downward from the rest length $l$ the spring has when supporting the mass $M$ vertically at rest (so $l=L_a+Mg/k$). Of course $k$ is the spring constant. Having checked these against other solutions online, I am fairly confident these are correct.
Part (b): I'm supposed to make some approximations, and I assume I'm supposed to obtain the equations
\begin{align*}
\ddot\theta + \frac{g}{l}\theta &= 0 \\
\ddot a + \frac{k}{M} a &= 0
\end{align*}
which are decoupled and have sinusoidal solutions. The obvious approximations to make are $\sin\theta\approx\theta$, $\cos\theta\approx 1$, $l+a\approx l$. But I also need to throw away the terms $2\dot\theta\dot a$ and $-(l+a)\dot\theta^2$. Is there a good reason this is legitimate? The best thing I can come up with is that if I substitute in the sinusoidal solutions to the exact equations, the terms $2\dot\theta\dot a$ and $-(l+a)\dot\theta^2$ have amplitudes squared whereas the other terms have just amplitudes. And since the amplitudes are small...is this reasoning correct? Is there better reasoning?
Part (c): I need a hint as to how to "solve the equations in part (a) to the next order". What substitutions or approximations should I make? Also: Maybe I can't understand this until I get an answer to the question I just asked, but: why is the question about resonance in this part? My first instinct is to compare the two angular frequencies from part (b) and see if they are close. If, say, they were equal, would we expect this to produce resonance behavior or not?
The bolded question is the one I would most like help on, but I would appreciate any thoughts on the other questions too.

A: The simplest way to keep track of all this is to introduce a "dummy" parameter $\epsilon$ which just "counts how small things are".  Thus start with your initial set of full equations
\begin{align*}
(l+a)\ddot\theta + 2\dot\theta\dot a + g\sin\theta &= 0 \tag{1} \\
\ddot a - (l+a)\dot\theta^2 + g(1-\cos\theta) + \frac{k}{M}a &= 0. \tag{2}
\end{align*}
and, under the approximations where $\theta$ and $a$ are small, replace
$\theta\to \epsilon\theta$ and $a\to \epsilon a$.  With this
$\dot \theta \to \epsilon \dot\theta$ etc so your equations become
\begin{align*}
(l+\epsilon a)\epsilon \ddot\theta + 
2\epsilon^2\dot\theta\dot a + g\sin\epsilon\theta &= 0 \tag{3}\\
\epsilon\ddot a - (l+\epsilon a)\epsilon^2 \dot\theta^2 + g(1-\cos\epsilon\theta) + \epsilon\frac{k}{M}a &= 0. \tag{4}
\end{align*}
You can then linearize your equations of motion, i.e. expand everything to terms linear in $\epsilon$, meaning you dump anything $\epsilon^2$ or above.
This immediately gives
\begin{align}
\epsilon l\ddot{\theta}+g\epsilon\theta&=0\\
\epsilon \ddot{a}+\epsilon\frac{k}{M}a&=0\, .
\end{align}
If you have done the job right, the counter $\epsilon$ just drops out (as it does here).  This makes it clear that terms in your first equation like $a\ddot{\theta}$ and $\dot{\theta}\dot{a}$ are of size $\epsilon^2$ and can be ignored.  Likewise this cleanly kills in your second equation the entire $(l+\epsilon a)\epsilon^2\dot{\theta}^2$ as it contains terms in $\epsilon^2$ and $\epsilon^3$, and also kills the $(1-\cos\epsilon \theta)$ term.
In problems such as the one you have, nothing replaces being systematic and carefully keeping track of what's expected to be small.

Edit 1: resonance requires an external driving force and the response of the system to resonant or near resonant input is not linear: the amplitude of the oscillation is dominated by a non-linear factor of $1/\sqrt{(\omega^2-\omega_0^2)^2+(\omega_0^2\omega^2/Q^2)}$ which is nearly singular when the quality factor $Q$ is small.  As a result the assumption of small amplitude of oscillation near a minimum of the effective potential is not valid and linearizing will produce nonsense (terrible approximate solution).

Edit 2: If you want "the next order", you need to solve by successive approximation.  I will sketch this only for the $\theta$ degree of freedom.  You assume
\begin{align}
\theta(t)&=A\cos\omega_\theta t + \epsilon F_\theta(t) \tag{5}\\
a(t)&=B\cos\omega_a t +\epsilon G_a(t) \tag{6}
\end{align}
(I trust the notation is obvious.)  Inserting Eq.(5) and Eq.(6) into the equation of motion for $\theta$, Eq.(3), and setting separate powers of $\epsilon$ to $0$; if my algebra is right, you get
\begin{align}
0&=\epsilon\,A \cos(\omega_\theta t)
\left(g - l\omega^2_\theta\right) \tag{7} \\
0&=\epsilon^2\left(-A B \omega^2_\theta \cos(\omega_a t)\cos(\omega_\theta t)+g F_\theta(t)+2 A B \omega_a\omega_\theta \sin(\omega_\theta t)\sin(\omega_a t)+ l \ddot{F_\theta}\right)\, . \tag{8}
\end{align}
Amazingly (and still provided my algebra is right!), this equation is still decoupled in $\epsilon^2$, in the sense it involves only $F_\theta$ but does not involve $G_a$ so you are left with a standard 2nd order differential equation equivalent to a driven oscillator.
The terms in $\cos(\omega_a t)\cos(\omega_\theta t)$ and $\sin(\omega_\theta t)\sin(\omega_a t)$ indicate that you must use
as ansatz 
$$
F_\theta(t)=C \cos\left((\omega_\theta-\omega_a)t\right)+ D \cos\left((\omega_\theta+\omega_a)t\right)\, ,
$$
(there might also be terms in 
$\sin\left((\omega_\theta\pm\omega_a)t\right)$ also depending on your initial conditions). 
You need to watch carefully for some resonance conditions, which would invalidate the assumption that the extra term $\epsilon F_\theta$ is small.  In general, one must also watch for the appearance of so-called secular terms. I don't think this happens here, but when this occurs, one must then write the first order frequencies $\omega_a$ and $\omega_\theta$ as a series of the form
$$
\Omega_\theta^2=\omega_\theta^2+\epsilon (\omega^{(1)}_\theta)^2+\ldots
$$
in a procedure called the Lindstedt-Poincare method.

A: b)It is legitimate. Because for linear motion, any second order terms or higher order terms will be canceled, such as $o(a^2),o(\theta^2),o(a\theta),o(\dot{a}\dot{\theta})$ etc. 
c) The next order means the any third order or higher order terms will be canceled. Just expand $sin\theta$ and $cos\theta$ to $o(\theta^3)$, you can obtain a coupled equtions. Resonance occurs only when the nature frequency of the system equals the driven frequency or approaches it. For coupled motion, stretching and angular motion drives each other, so you can compare the two frequency.
