# Adding small correction term to ODE solution

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.

• What kind of orbit do you expect to get? Is this just for a time step of $dt$ Have you solved it numerically? If a certain $r_0$ is important and you are looking to find some sort of deviation from a circle the shouldn't you use a polar coordinate attempt? The approximation in (2) might not hold for long, so if there is a periodic orbit you will be away fro $r_0$ most of the time/ Jun 10, 2013 at 23:49
• Basically, I'm sampling the magnetic field at every $\mathbf{r}_{\circ}$ and using this constant value in the known expression for the trajectory of a particle in such a field to find its location a short time later and repeating. The field is actually changing in this short time, but by taking smaller and smaller steps, it gets better. However, it isn't quite good enough, hence why I want to add a correction term. Jun 10, 2013 at 23:55
• You are neglecting the $q\mathbf{E}$ part of the Lorentz force on the particle, while your comment leads me to believe that you have time-varying magnetic fields, which of course lead to time varying electric fields. What kind of time-variation / frequency are you concerned with? Are you certain that you can neglect the electric field's ability to accelerate the particle? Jun 13, 2013 at 15:24
• Can you specify if you're trying to come up with a first order correction for the time variation of the fields, or a correction for spatial variation of the fields? Your Taylor expansion of the field hints at correction for spatial gradients, but then you try to derive time corrections? Jun 13, 2013 at 15:42
• @rajb245 I'm only concerned with particles in a magnetic field right now. Also, I'm only looking at magnetic fields with spatial variations, not temporal ones. The correction I derived isn't a correction to the magnetic field, it's a correction to the positions and velocities of the particle, which are functions of time. The reason that I need a correction is because the field varies spatially, so that's why I needed to expand $\mathbf{B}$. Jun 13, 2013 at 15:56

$$\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.