Adding small correction term to ODE solution

Question

Let $\mathbf{r}(t) = [x(t), y(t), z(t)]$ and $\mathbf{v}(t) = \frac{d}{dt}\mathbf{r}(t)$. I'm trying to solve $$ \frac{d}{dt}\mathbf{v}=\frac{q}{m}(\mathbf{v}\times\mathbf{B}) \tag{1} $$ where $q$ and $m$ are real constants and $\mathbf{B}(\mathbf{r})$ is an arbitrary vector field. Actually, $\mathbf{B}$ is a magnetic field and $(1)$ is the equation for the motion of a particle in such a field, but that isn't vitally important for the purposes of my question.

If $\mathbf{B}$ is a constant,$\mathbf{B}_{0}$, then we can find the exact solution of $(1)$ with initial conditions $\mathbf{r}_{\circ}$ and $\mathbf{v}_{\circ}$; call these solutions $\mathbf{r}_{0}$ and $\mathbf{v}_{0}$. However, $\mathbf{B}$ usually isn't a constant, it's usally a total mess, and $\mathbf{r}_{0}$ and $\mathbf{v}_{0}$ are poor approximations. In order to make the solutions a little bit better, one could consider expanding $\mathbf{B}$ about $\mathbf{r}_{\circ}$ to first order in $x, y$ and $z$, introducing the term $\mathbf{B}_{1}$, and then attempting to solve $(1)$. Even this is usually impractical, so I really only want the lowest order term of the solution. In short, my idea is to sub $\mathbf{B} = \mathbf{B}_{0} + \mathbf{B}_{1}$ and $\mathbf{v}=\mathbf{v}_{0}+\mathbf{v}_{1}$ into $(1)$, use the fact that $\mathbf{v}_{0}$ is known, and discard as much as possible to find a really low-order correction $\mathbf{v}_{1}$. The issue is that I'm not totally sure how to actually get that correction term, ie. what I'm allowed to throw away. The details are below; I would really appreciate it if someone could read it and tell me if what I did is valid and where I should go next.

Derivation

Consider a particle of mass $m$ and charge $q$ in an arbitrary magnetic field $\mathbf{B}$ at position $\mathbf{r}_{\circ}(t_{\circ})=[x_{\circ}, y_{\circ}, z_{\circ}]$ with velocity $\mathbf{v}_{\circ}(t_{\circ})=[v_{x{\circ}}, v_{y{\circ}}, v_{z{\circ}}]$. For $\mathbf{r}(t) = [x(t), y(t), z(t)]$ and $\mathbf{v}(t) = \dot{\mathbf{r}}(t)$, the particle's trajectory satisfies $(1)$. We may Taylor expand the magnetic field about $\mathbf{r}_{\circ}$ as follows: \begin{equation} \mathbf{B}(\mathbf{r})\approx \mathbf{B}(\mathbf{r}_{\circ}) + (\mathbf{r}- \mathbf{r}_{\circ})\cdot\nabla \mathbf{B}(\mathbf{r}_{\circ}) \tag{2} \end{equation} Where $\mathbf{r}\cdot\nabla\mathbf{B}(\mathbf{r}_{\circ}) = (x\partial_{x} + y\partial_{y} + z\partial_{z})\mathbf{B}(\mathbf{r})|_{\mathbf{r}=\mathbf{r}_{\circ}}$. Define $\mathbf{B}_{0}=\mathbf{B}(\mathbf{r}_{\circ})$ along with $\mathbf{B}_{1} = \mathbf{B}_{10} + \mathbf{B}_{11}= -\mathbf{r}_{\circ}\cdot\nabla \mathbf{B}(\mathbf{r}_{\circ}) + \mathbf{r}\cdot\nabla \mathbf{B}(\mathbf{r}_{\circ})$. Further separate $\mathbf{r}=\mathbf{r}_{0}+\mathbf{r}_{1}$ and $\mathbf{v}= \mathbf{v}_{0} + \mathbf{v}_{1}$, where $\mathbf{v}_{0}$ satisfies $(1)$ with field $\mathbf{B}_{0}$, the explicit solution of which is known since the field is simply a constant vector. By subbing $\mathbf{v}$ and $\mathbf{B}$ into $(1)$ and canceling the known solution, we have \begin{equation} \frac{d}{dt} \mathbf{v}_{1} = \frac{q}{m}(\mathbf{v}_{0}\times\mathbf{B}_{1}) + \frac{q}{m}\left(\mathbf{v}_{1}\times(\mathbf{B}_{0}+\mathbf{B}_{1})\right) \tag{3} \end{equation} Then note that we must have $\mathbf{r}_{1}(0)=\mathbf{v}_{1}(0) = \mathbf{0}$ since $\mathbf{r}_{0}(0) = \mathbf{r}_{\circ}$ and $\mathbf{v}_{0}(0) = \mathbf{v}_{\circ}$, ie. the initial conditions are already taken care of. Therefore, the lowest order correction in time must be, for $\mathbf{y}=[\alpha, \beta, \gamma]$, $\mathbf{r}_{1}=\mathbf{y}t^2 \implies \mathbf{v}_{1}=2\mathbf{y}t\implies \frac{d}{dt}\mathbf{v}_{1}=2\mathbf{y}$; by subbing these into $(3)$, we can solve for $\mathbf{y}$. The first-order solution discards all terms with $t^2$ dependence, and therefore we should set $\mathbf{r}_{1}\approx \mathbf{0}$. Also, applying this idea to $\mathbf{B}$, we have $\mathbf{B}_{11}=\mathbf{r}\cdot\nabla\mathbf{B}(\mathbf{r}_{\circ}) = (\mathbf{r}_{0} + \mathbf{r}_{1})\cdot\nabla\mathbf{B}(\mathbf{r}_{\circ}) \approx \mathbf{r}_{0} \cdot\nabla\mathbf{B}(\mathbf{r}_{\circ})\approx (\mathbf{r}_{\circ}+\mathbf{v}_{\circ}t) \cdot\nabla\mathbf{B}(\mathbf{r}_{\circ})$, so that $\mathbf{B}_{1} \approx t\mathbf{v}_{\circ}\cdot\nabla\mathbf{B}(\mathbf{r}_{\circ})$.

Here's the issue. If we now sub everything into $(3)$, we have $\mathbf{B}_{1}\propto t$ and $\mathbf{v}_{1}\propto t$, so the left side is a constant and the ride side is $\propto t$ after discarding quadratic terms. How can you get a condition on $\mathbf{y}$ from that?

Edit:

There might be a simple fix. It seems like if I postulate a solution of the form $r_{1}=\mathbf{y}t^3$ rather than $r_{1}=\mathbf{y}t^2$, then this can work. In this case, the term in $(3)$ with $\mathbf{v}_{1}$ is still discarded, but now we have that the left side is $6\mathbf{y} t$. Writing $\mathbf{v}_{0}\approx\mathbf{v}_{\circ}$, we have a linear term on the right side that can be matched with this, giving $$ \mathbf{y}= \frac{q}{6m}\mathbf{v}_{\circ}\times(\mathbf{v}_{\circ}\cdot\nabla)\mathbf{B}(\mathbf{r}_{\circ}) $$ I'd still like someone to confirm this logic though.

Answer 1

This is a partial answer that should help you a bit. Your unperturbed equation is

$$ \dot{\mathbf{v}}_0 = \frac{q}{m}\mathbf{v}_0\times\mathbf{B}_0\tag{1} $$

When $\mathbf{B}_0=(B_{0x},B_{0y},B_{0z})^T$ is not a function of space, this reduces to

$$ \dot{\mathbf{v}}_0= \frac{q}{m} \left( \begin{array}{ccc} 0 & B_{0z} & -B_{0y} \\ -B_{0z} & 0 & B_{0x} \\ B_{0y} & -B_{0x} & 0 \\ \end{array} \right) \mathbf{v}_0 $$

Which can be solved by matrix exponentials, $\mathbf{v}_0(t)=\exp \left(\frac{q}{m}\underline{\mathbf{B}}t\right)\mathbf{v}_0(0)$, where $\underline{\mathbf{B}}$ is the matrix above. You get velocity orbits that are periodic in a plane, which corresponds periodic motion in two directions superimposed on a constant velocity motion in the perpendicular direction. This is a simple helical motion.

Now the perturbation is to let $\mathbf{v}_0\rightarrow \mathbf{v}_0+\mathbf{v}_1$ and $\mathbf{B}_0\rightarrow \mathbf{B}_0+\mathbf{B}_1$. Inserting into $(1)$ I get: $$ \dot{\mathbf{v}}_0 + \dot{\mathbf{v}}_1= \frac{q}{m}\left( \mathbf{v}_0\times\mathbf{B}_0+ \mathbf{v}_0\times\mathbf{B}_1+ \mathbf{v}_1\times\mathbf{B}_0+ \mathbf{v}_1\times\mathbf{B}_1 \right) $$ You can neglect the product of the two perturbation terms ($\mathbf{v}_1\times\mathbf{B}_1$), since this would be a "second order" kind of term. Usually you carry around a small parameter $\varepsilon$ in these kinds of calculations, in which case that term would be $O(\varepsilon^2)$ and neglectable, and then the $\varepsilon$ cancels out everywhere anyway. Neglecting that term and cancelling the solution that satisfies $(1)$, I get

$$ \dot{\mathbf{v}}_1= \frac{q}{m}\left( \mathbf{v}_0\times\mathbf{B}_1+ \mathbf{v}_1\times\mathbf{B}_0 \right) $$

Now try applying your ideas to this equation. It seems to me that this is simply a forced version of $(1)$ (with the known $\mathbf{v}_0\times\mathbf{B}_1$ playing the role of the forcing) and could be solved exactly in terms of a convolution integral. The correction wouldn't be a simple low order function of $t$, but some pretty complicated function involving matrix exponentials. Nonetheless, you could probably work out the first order effect of field gradients using this technique.

Edit

So upon thinking on it, I see why you go with a low order temporal approximation, so that you can march in short time steps. This is technically another kind of asymptotic approximation, different from the perturbation method outlined above. But, if you're interested in numerical solutions, I must say a standard ODE solver would be a better choice. MATLAB's ODE routines can handle this system directly.