I have the following coupled differential equations also known as Guiding Center Approximation. It is used to explain the position- and velocity change of particles (electrons and protons, N = 1000) moving in plasma. I am trying to either solve this using the Runge-Kutta 4th order method or the built-in ODE45 in Matlab.
The equations are as follows:
\begin{align} {d u_{||} \over {dt}} &= {q \over m} E_{||} - {\mu \over {\gamma m B}} (\mathbf{B} \cdot \nabla) B \\ {{d \mathbf{R}} \over dt} &= {u_{||} \over {\gamma B}} \mathbf{B} + {{\mathbf{E} \times \mathbf{B}} \over B^2} + { \mu \over {\gamma q B^2}} \mathbf{B} \times \nabla B + {{ m u^2_{||}} \over {\gamma q B^4}} \mathbf{B} \times (\mathbf{B} \cdot \nabla) \mathbf{B}\\ \mu & = {{m u^2_{\perp}} \over {2B}} = \mathrm{constant} \end{align}
where,
- || is the velocity component parallel to the magnetic field
- $u$ is the particle velocity
- $q$ is the charge
- $m$ is the mass of the particle
- $\mathbf{E}$ is the electric field vector
- $\mu $ is the constant determined in the last equation above
- $\mathbf{B}$ is the magnetic field vector
- $\mathbf{R}$ is the position of the particle
- $\gamma $ is the Lorentz factor given by $\gamma = (1-\frac{v^2}{c^2})^{-\frac{1}{2}}$
I am given the initial position and velocity of each particle. Furthermore I'm given $\mathbf{B}$ & $\mathbf{E}$. I'm also given the initial dt. When $B$ is used without being bold it's either the norm of $\mathbf{B}$ or the component going in the direction of $\mathbf{B}$ and same goes for the electric field.
\begin{align} \mathbf{B} & = {B_0 \over L} (x, y, -2 z)\\ \mathbf{E} & = E_0 L \left( { -y^2 \over {(x^2+y^2)^{1.5}}}, {xy \over {(x^2+y^2)^{1.5}}}, 0 \right) \end{align}
To solve these coupled equations I tried using the RK4 method. What I'm unsure of is how to write my equations in a way of $f(t_n,y_n)$
where,
- $k_1 = f(t_n,y_n) $
- $k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot k_1)$
- $k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot k_2)$
- $k_4 = f(t_n + h, y_n + h \cdot k_3)$
- $y_{n+1} = y_n + \frac{h}{6} \cdot (k_1 + 2\cdot k_2 +2\cdot k_3 + k_4 )$
- $t_{n+1} = t_n + h$
And furthermore, how do I using Matlab, calculate these two equations at the same time (as I assume I must). Lastly, would it be easier or more accurate to use the builtin ODE45?