Intuition behind variation of parameters method for solving differential equations I have used the variation of parameters method (and have been taught it, although not hugely in depth) and I was wondering if I've understood the intuition behind it. In particular I've been thinking about the method for second order ODEs: $$a_{2}(x)y''(x)+a_{1}(x)y'(x)+a_{0}(x)y(x)=f(x)$$
Does the motivation for considering a particular solution of the form $$y_{p}(x)=u_{1}(x)y_{1}(x)+u_{2}(x)y_{2}(x)$$ (where $y_{1},\;y_{2}$ are solutions to the corresponding homogeneous equation) because $(y_{p}/y_{1})\neq\text{ constant}$ and likewise $(y_{p}/y_{2})\neq\text{ constant}$, as otherwise we would just obtain the complementary solution again. This suggests that both are instead functions of $x$, i.e. $(y_{p}/y_{1})=u_{1}(x)$ and $(y_{p}/y_{2})=u_{2}(x)$, leading to the form of the ansatz I gave above?
Secondly, is the reason why we place a further constraint on the form of $y_{p}$ (other that it be a solution to the original ODE) because, in principal, there will be an infinite number of particular solutions of the form $$y_{p}(x)=u_{1}(x)y_{1}(x)+u_{2}(x)y_{2}(x)$$ but we only require one particular solution, and so by imposing an additional constraint we have two equations for the two unknowns $u_{1}$ and $u_{2}$, thus enabling us to uniquely determine a solution for each of them and subsequently enabling us to find a general solution to the original ODE?!
 A: Here is one way of getting the variation of parameters formulas. Firstly, some motivation mixed with a bit of jargon: 
We think of $a_{2}(x)y''(x)+a_{1}(x)y'(x)+a_{0}(x)y(x)$ as a result of applying to $y(x)$ an operation $D=a_{2}(x)\frac{d^2}{dx^2 }+a_{1}(x)\frac{d}{dx}+a_{0}(x) Id$. This is a linear operation (aka a linear differential operator), meaning that $D(c_1y_1(x)+c_2y_2(x))=c_1 D(y_1(x))+c_2 D((y_2(x))$. Our equation is then $D (y(x))=f(x)$.
Linearity ensures that if $y_a(x)$ and $y_b(x)$ solve $Dy_a=f_a$ and $Dy_b=f_b$ then $y_a+y_b$ solves $Dy=f_a+f_b$. So if we can decompose our inhomogeneous term $f(x)$ as a sum $f=f_a+f_b$ and solve the two "pieces" $Dy=f_a$ and $Dy=f_b$ we  would be done - just add the two solutions and get a solution for $f$.
The key idea now (this is known as Green's function method, or perhaps sometimes Duhamel's principle) is to think of $f$ as being a "sum" of a continuum of delta functions $f(x)=\int f(t)\delta(x-t) dt$. Thus we want to solve $Dy(x)= \delta(x-t)$ for each $t$ (these solutions are the Green's functions for our equation), to get solutions $y_t(x)=y(t,x)$, and then write our solution to $Dy=f$ as the "sum" over t $y(x)=\int f(t) y_t (x) dt = \int f(t) y (t , x) dt $. 
Continuing in this way (without paying much attention to technical details or introducing the theory of distributions needed to provide such details and make the above rigorous), we are now tasked with finding $y_t(x)$ i.e. solving $Dy_t(x)=\delta(x-t)$. The solution is of course not unique, but is only unique up to adding $y_h$ - an arbitrary solution of the homogeneous equation $Dy=0$. Note that in fact for $x>t$  and for $x<t$ we have $\delta(x-t)=0$, and our $y_t$ will coincide with some $y_h$. But it can not be the same $y_h$, or we won't get $\delta(x-t)$, just $0$ everywhere. So we need to jump at $x=t$. Namely, we can start with $y=0$ for $x<t$ and continue as some $y_h$ after. 
Now, any  $y_h$ is a combination of fixed homogeneous solutions $y_1$ and $y_2$, i.e. $y_h(x)=c_1y_1(x)+c_2y_2(x)$. So our Green's function $y_t=0$ for $x<t$ and $y_t(x)=c_1y_1(x)+c_2y_2(x)$ for $x>t$. What should $c_1, c_2$ be? We want $Dy_t(x)= a_{2}(x)y_t''(x)+a_{1}(x)y_t'(x)+a_{0}(x)y_t(x)= \delta(x-t)$. So $y_t$ should be continuous, so that $y_t'$ has only a step discontinuity at $x=t$, and $y_t''$ only a delta (and not $\delta'(x-t)$). Continuity imposes $c_1y_1(t)+c_2y_2(t)=0$.  We also want to have a unit delta jump and not bigger or smaller. The size of the jump is controlled by $a_2 \frac{d^2}{dx^2}$ at $x=t$ and so is $a_2(t)[ c_1y'_1(t)+c_2y'_2(t)]$.
Equating it to 1 imposes $ c_1y'_1(t)+c_2y'_2(t)=\frac{1}{a_2(t)}$
To summarize:
The Green's function at $t$ is [$0$ for $x<t$, then  $c_1(t)y_1(x)+c_2(t)y_2(x)$ for $x>t$]
with $c_1(t), c_2(t)$ subject to 
$c_1(t)y_1(t)+c_2(t)y_2(t)=0$
$ c_1y'_1(t)+c_2y'_2(t)=\frac{1}{a_2(t)}$
Compare with the variation of parameters: $c_1=\frac{u_1'}{f}, c_2=\frac{u_2'}{f}$
And the solution to $Dy=f$ is
$y(x)=\int f(t) y_t (x) dt = \int f(t) y (t , x) dt $
Plugging in, we get
$y(x)=\int_{x>t} f(t) (c_1(t)y_1(x))+c_2(t)y_2(x)) dt=[\int_{-\infty}^{x} f(t)c_1(t)) dt] y_1(x)+[\int_{-\infty}^{x} f(t)c_2(t))dt] y_2(x)=u_1(x)y_1(x)+u_2(x)y_2(x)$.
So we recover $y$ as a combination of $y_1$ and $y_2$ with function coefficients, and we know exactly which relations the derivatives of these coefficients should satisfy (they have to combine to give Green's functions times $f(x)$). Writing out these relations we get the same formulas as for the variation of parameters formulas, which we can then justify independently in an ad hoc manner.
Of course this generalizes in a straightforward manner to higher order linear equations - the $c_i$s will have more linear `matching of derivatives' constraints to satisfy to give a Green's function, and the rest is the same.

A: This is an extremely good and clever question.
Amazingly, the raised issue, the freedom in our choice of the constraint, had been waiting for discussion for more than two centuries, and was addressed and exploited in the literature only recently.
I begin with an elementary example borrowed from this paper. It was proposed by my coauthor William Newman, so I call it Newman's example.
Consider a harmonic oscillator disturbed by a force $\Delta F(t)$:
$$
\ddot{x} + x = \Delta F(t)\quad,\quad\mbox{with $x(0)$ and $\dot x(0)$ known}\;,\qquad(1)
$$
and seek its solution by using the ansatz
$$
x = C_1(t)\,\sin t + C_2(t)\,\cos t\;.\qquad\qquad (2) 
$$
Together, the former and the latter equations will yield:
$$
\dot{x} = [\dot{C}_1(t)\,\sin t + \dot{C}_2(t)\,\cos t] + C_1(t)\,\cos t - C_2(t)\,\sin t\;.\qquad\qquad(3)
$$
It is common to impose at this point an analogue to the Lagrange constrain, i.e., to set $\,[\dot{C}_1(t)\,\sin t + \dot{C}_2(t)\,\cos t]\,$ equal to zero. This step is convenient but not obligatory. We instead may agree on
$$
\dot{C}_1(t)\,\sin t + \dot{C}_2(t)\,\cos t =\phi(t)\;,\qquad\qquad(4)
$$
$\phi(t)$ being an arbitrary function of time. This results in
$$
\ddot{x} = \dot{\phi} + \dot{C}_1(t)\,\cos t -\dot{C}_2(t)\,\sin t - C_1(t)\,\sin t - C_2(t)\,\cos t\;\,,\qquad\qquad(5)
$$
summation whereof with (2) gives us
$$
\ddot{x} + x = \dot{\phi} + \dot{C}_1\,\cos t - \dot{C}_2\,\sin t\;.\qquad\qquad(6)
$$
Plugging this expression into (1) produces, when combined with identity (4), the following system:
$$
\dot{\phi} + \dot{C}_1\,\cos t - \dot{C}_2\,\sin t = \Delta F(t)\qquad\qquad(7)
$$
$$
\dot{C}_1(t)\,\sin t + \dot{C}_2(t)\,\cos t =\phi(t)\;.\qquad\qquad(8)
$$
This leads to
$$
\dot{C}_1 = \Delta F\,\cos t - \frac{d}{dt}(\phi\,\cos t)\;,\qquad\qquad(9)
$$
$$
\dot{C}_2 = - \Delta F\,\sin t + \frac{d}{dt}(\phi\,\sin t)\;,\qquad\qquad(10)
$$
integration of which system entails:
$$
C_1 = \int^t \Delta F(t')\,\cos t'\,dt' - \phi\,\cos t + a_1\;,\qquad\qquad(11)
$$
$$
C_2 = - \int^t \Delta F(t')\,\sin t'\,dt' _ \phi\,\sin t + a_2\;,\qquad\qquad(12)
$$
$a_1$ and $a_2$ being constants.
Substitution of the two last expressions into (9) results in a complete cancellation of the arbitrary $\phi(t)$ term:
$$
x = C_1(t)\,\sin t + C_2(t)\,\cos t = 
$$
$$
- \cos t \int^t \Delta F (t')\,\sin t'\,dt' + \sin t \int^t \Delta F(t')\,\cos t'\,dt' + a_1\,\sin t + a_2\,\cos t\;.\qquad\qquad(13)
$$
Naturally, the physical trajectory $x(t)$ comes out invariant under the choice of the gauge function $\phi(t)$.
Does this imply that gauge freedom is unimportant? No way! Sometimes it helps a lot in analytical calculations. On other occasions, a clever choice of the gauge function(s) helps to mitigate the numerical error in computations, see this work.
Now, some details and a bit of history.
Let us recall the geometric idea underlying the Variation of Parameters (VOP) method in its initial form, as suggested by Euler and Largange. We need to model a perturbed orbit by a sequence of "simple" curves. If the orbit is bound, it is convenient (not obligatory - but convenient) to choose these "simple curves" as ellipses sharing one of their foci. One such ellipse can be parameterised, in the inertial Cartesian frame, as
$$
{\bf r} = {\bf f}(t,\,C_1,\;.\;.\;.\;C_6)\;,\qquad\qquad (14)
$$
where $C_i$ are six parameters. A popular choice of these parameters is the so-called Keplerian elements: the major semiaxis $a$, the eccentricity $e$, the three Euler angles (the longitude of the node $\Omega$, the inclination $i$ of the orbit on the fiducial plane, argument of the pericentre $\omega$), and the initial condition (the value ${\cal{M}}_0$ at the initial time $t_0$). Why does the initial condition enter this set? Because we astronomers are interested not only in the geometric shape of a "simple" orbit, but also in the position of the orbiter on it. So, equation (14) gives not just a "simple" orbit, but a "simple" solution, with the initial position included.
In the above equation, $\bf f$ is the implicit function serving as a solution to the unperturbed equation
$$ 
\ddot{\bf r} = {\bf F}({\bf r})\;,\qquad\qquad(15)
$$
where ${\bf F}({\bf r})$ is the inverse-square force. Euler and Largange were interested in this particular example, and this is why in this case the unperturbed force ${\bf F}$ is a function of ${\bf r}$ solely.
Now, include a perturbation $\Delta \bf F$. We have to solve a more complex equation of motion
$$ 
\ddot{\bf r} = 
{\bf F}({\bf r}) + \Delta {\bf F} ({\bf r},\,{\bf\dot r})\;.\qquad\qquad(16)
$$
Watch my hands: while Euler and Lagrange were interested only in position-dependent perturbation, I have set $\Delta \bf F$ to depend also on velocities, because the modern celestial mechanics has to deal with such perturbations. Such perturbations emerge, e.g., when we have to consider motions in noninertial frames or when relativistic corrections come into play. This detail, however, will not change the central idea of our story.
To solve equation (16) by the VOP method, we endow the "constants" with time-dependencies of their own:
$$
{\bf r} = {\bf f}(t,\,C_1(t),\;.\;.\;.\;C_6(t)\,)\;,\qquad\qquad (17)
$$
while keeping the functional form of $\bf f$ unchanged. In simple words, we model the perturbed orbit with a sequence of "simple solutions" parameterised with $t$. At the risk of sounding repetitive, I again highlight that a "simple solution" is not just an ellipse with one focus fixed, but an ellipse taken together with an initial condition. And why with one focus fixed? -- because it is natural to assemble a perturbed orbit of "simple solutions" that are elliptic orbits about the fixed point, Sun. Each such simple orbit donates one point to the perturbed trajectory, the actual solution.
At this point, Lagrange faced a mathematical difficulty. The first time differentiation rendered
$$
{\bf\dot r} = \frac{\partial \bf r}{\partial t} + {\bf \Phi}\quad,\qquad\Phi \equiv \sum_{j=1}^6\frac{\partial \bf f}{\partial C_j}\,\dot C_j\;,\qquad\qquad(18)
$$
while the second differentiation gave him
$$
{\bf\ddot r} = \frac{\partial^2 \bf r}{\partial t^2} + 
\sum_{j=1}^6\frac{\partial^2\bf f}{\partial t\,\partial C_j}\,\dot{C}_j + {\bf \dot\Phi}\;.
$$
The insertion thereof into (16) left Lagrange with
$$
\sum_{j=1}^6\frac{\partial^2\bf f}{\partial t\,\partial C_j}\,\dot{C}_j + {\bf \dot\Phi} = \Delta {\bf F}\;.\qquad\qquad(19)
$$
This didn't look good, because at his disposal Lagrange had only three equations, the x,y,z projections of (19). These were not enough to obtain solutions for the six functions $C_j(t)$.
To overcome this difficulty, he set, by hand, what we now call the Lagrange constraint:
$$
{\bf \Phi} = 0\;.\quad\quad\quad(20)
$$
The physical meaning of this constraint is simple: it postulates that we model the perturbed trajectory by a sequence of ellipses that are tangent to the trajectory. Each ellipse donates one point to the trajectory -- and is tangent to the trajectory at that point. It does not cross the trajectory. Lagrange termed these tangent ellipses osculating (osculare being Latin for kissing).
Mathematically, the constraint is arbitrary. As was emphasised in several recent publications (this, this, and this), it is possible to make $\bf\Phi$ and arbitrary function of $\bf r$. Its presence will alter the mathematical form of the resulting solution, but will not change the actual physical trajectory. Geometrically, the trajectory will be assembled of points provided by nonosculating ellipses, which will now be permitted to cross the trajectory.
Hence we are facing an example of gauge freedom.
Lagrange surely realised the presence of these three degrees of freedom, the three projections of ${\bf\Phi}({\bf r})\,$. However, neither he nor Poincare not other great celestial mechanicians of the past bothered to look into this issue. My guess is that they did not expect any calculational gains from using alternative gauges.
Gains (big ones) emerged only when it became necessary to solve the problems with velocity-dependent perturbations $\Delta {\bf F}({\bf r},\,{\bf\dot r})$. It has turned out that in such situations the imposition of the Lagrange constraint entails very complex equations for the orbital parameters -- while the choice of a non-Lagrange gauge simplifies the equations drastically.
I would also mention that the attitude mechanics, i.e., the problem of perturbed rotation of an unsupported body is mathematically analogous to the mechanics of orbits. Just as a perturbed orbit gets assembled of points donated by "simple solutions" (usually, conics) -- so a perturbed rotation of a solid body gets assembled of simple "Eulerian cones", each of which is an unperturbed rotation. Like in orbital mechanics, so in attitude mechanics the "Eulerian cones"  can be osculating or nonosculating, see this paper.
