Motion in a time-dependent uniform magnetic field Suppose you have an infinite solenoid generating an uniform magnetic field inside.  The field is oriented along the solenoid axis : unit vector $\vec{\bf n}$.  The field intensity varies linearly with time between $t_1$ and $t_2 = t_1 + \Delta t$, like this (we neglect all electromagnetic waves here) :
\begin{equation}\tag{1}
B(t) = B_1 + \lambda \, (B_2 - B_1)(t - t_1),
\end{equation}
where $B_1$ is the constant magnetic field for time $t < t_1$, $B_2$ is the constant magnetic field for time $t > t_2 = t_1 + \Delta t$, and $\lambda = 1/\Delta t$.  The time variation of this field generates an induced electric field inside the solenoid, for the same time interval (from $t_1$ to $t_2$) :
\begin{equation}\tag{2}
\vec{\bf E} = -\, \frac{1}{2} \; \lambda \; \Delta B \; \vec{\bf n} \times \vec{\bf r},
\end{equation}
where $\Delta B = B_2 - B_1 > 0$ is a simple constant (the magnetic field is increasing in the solenoid).  Take note that the unit vector $\vec{\bf n}$ is also a constant (the solenoid's axis).
Now, you drop a positive charge particle inside the solenoid : $q > 0$, with any position and initial velocity.  The equation of motion is this :
\begin{equation}\tag{3}
\frac{d \vec{\bf p}}{d t} = q \, \vec{\bf E} + q \, \vec{\bf v} \times \vec{\bf B},
\end{equation}
where $\vec{\bf p} = \gamma \, m \, \vec{\bf v}$ is the relativistic linear momentum of the particle.  I'm not interested in solving analytically this equation (I've done it numerically using Mathematica.  The 3D curves are pretty !).  Now, the problem is this :
How can we find analytically the final energy at time $t > t_2$, as a function of the field intensities $B_1$, $B_2$ and the initial velocity (or energy) at time $t < t_1$ ?
I know that there is at least one exact constant of motion for this problem :
\begin{align}
\mathcal{J} &= \vec{\bf n} \cdot \big( \vec{\bf r} \times (\vec{\bf p} + q \, \vec{\bf A}) \big), \\[18pt]
&= \vec{\bf n} \cdot \big( \vec{\bf r} \times \gamma \, m \, \vec{\bf v} + \frac{q}{2} \; B(t) \, \vec{\bf r} \times (\vec{\bf n} \times \vec{\bf r})\big), \tag{4}
\end{align}
where $\vec{\bf A}$ is the potential-vector :
\begin{equation}\tag{5}
\vec{\bf A} = \frac{1}{2} \; B(t) \, \vec{\bf n} \times \vec{\bf r}.
\end{equation}
We could also try to use the kinetic energy theorem (the magnetic field doesn't do any work) :
\begin{align}
\Delta K = W_{\text{em}} &= \int_{t_1}^{t_2} q \, \vec{\bf E} \cdot \vec{\bf v} \; dt \\[18pt]
&\equiv -\, \frac{q}{2} \; \lambda \, \Delta B \int_{t_1}^{t_2} \vec{\bf n} \cdot (\vec{\bf r} \times \vec{\bf v}) \, dt, \tag{6}
\end{align}
but unfortunately it's not helping since I don't know how to evaluate this integral (take note that the vector $\vec{\bf r} \times \vec{\bf v}$ isn't conserved here, and it's not exactly the particle's angular momentum since the relativistic $\gamma$ factor is missing).  However, we recognize the time integral of the particle's magnetic moment $\vec{\boldsymbol{\mu}}(t)$ :
\begin{equation}\tag{7}
\vec{\boldsymbol{\mu}}(t) = \frac{q}{2} \; \vec{\bf r}(t) \times \vec{\bf v}(t),
\end{equation}
so we could write the following kinetic energy variation, but it doesn't help much :
\begin{equation}\tag{8}
\Delta K = -\, \langle \, \vec{\boldsymbol{\mu}} \, \rangle \cdot \Delta \vec{\bf B}.
\end{equation}
The constant of motion $\mathcal{J}$ isn't of any help in this case, even if the motion is constrained to the plane orthogonal to $\vec{\bf n}$ (i.e. motion in the solenoid's cross section).
Any suggestion to find the kinetic energy variation $\Delta K$ ?
I also suspect that there may be another exact conserved quantity for this problem (total energy ?  total magnetic flux on the particle's path ?).  What may be the other conserved quantity ?

Here's a typical trajectory in the solenoid's orthogonal plane :
A picture in the plane http://s10.postimg.org/87h0m0849/motion.jpg[/img]
The large circle is the initial motion at time $t < t_1$ (classical circular motion, of radius $r_1 = \gamma_1 \, m \, v_1/ q \, B_1$).  The small circle inside is the final motion at time $t > t_2$ (another classical circular motion around the final magnetic field lines, of radius $r_2 = \gamma_2 \, m \, v_2/ q \, B_2$).  The path between both circles is the effect of the time varying magnetic field and of the induced electric field (which accelerates the particle : $v_2 > v_1$).  I need to analytically find the energy variation from the large circle to the smaller one, to get the final radius $r_2$ (since we don't know the final linear momentum $p_2 = \gamma_2 \, m \, v_2$).
Here's another picture to show some typical trajectories in 3D.  The drift to the center occurs during the transition $B_1 \Rightarrow B_2 > B_1$ :
A picture in 3D http://s11.postimg.org/yarbwa0ab/induction2.jpg[/img]
The drift is caused by the induced electric field, which is kicking the particles with a local drift velocity $\vec{\bf v}_d = \vec{\bf E} \times \vec{\bf B}/B^2$.

Complement :
This may be interesting.  If we consider non-relativistic motion in a plane only (orthogonal to the magnetic field lines), using polar coordinates give the following radial differential equation :
\begin{equation}
\ddot{\rho} + \omega^2 \, \rho = \frac{\mathcal{J}^2}{\rho^3},
\end{equation}
where $\mathcal{J}$ is the constant of motion defined above (per unit of mass) and $\omega = q B(t)/2m$ is the Larmor angular frequency.  This diff. equation is hard to solve, especially since $\omega$ depends on $t$.  The angular part is given by this equation :
\begin{equation}
\dot{\vartheta} = \frac{\mathcal{J}}{\rho^2} - \omega.
\end{equation}
 A: The problem simply is not integrable and thus we cannot generally trace the evolution of all phase-space variables analytically. The easiest way to describe it is via a Hamiltonian in cylindrical coordinates
$$H = \frac{(p_\phi - c A_{\phi})^2}{2m \rho^2} + \frac{1}{2m} (p_\rho^2 + p_z ^2)$$
where $A_{\phi}=B(t) \rho /2$ (you can easily see that there is no need for electric potential $\Phi$). The obvious symmetries are rotational leading to the conservation of $p_\phi= m\dot{\phi}\rho^2 + c A_\phi$, and translational, leading to the consevation of $p_z=\dot{z}$. Unfortunately, there are still two degrees of freedom, $\rho$ and $t$ which means that general initial conditions may lead even to chaotic scattering. 
If the system truly exhibits chaotic scattering, it is a proof of the fact that you cannot find a general analytical formula. However, sometimes it happens that a system has a "hidden" additional integral. There is no easy way to discern between two cases. I think the easiest thing you can do is resort to some kind of approximation such as assuming $\Delta t$ is small, or on the other hand, that $\Delta t$ is large and you can thus integrate the energy loss as adiabatically evolving through the orbits in the time-independent system.
