1
$\begingroup$

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,

  1. $k_1 = f(t_n,y_n) $
  2. $k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot k_1)$
  3. $k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot k_2)$
  4. $k_4 = f(t_n + h, y_n + h \cdot k_3)$
  5. $y_{n+1} = y_n + \frac{h}{6} \cdot (k_1 + 2\cdot k_2 +2\cdot k_3 + k_4 )$
  6. $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?

Improve this question
$\endgroup$
  • 2
    $\begingroup$ In my experience, the amount of work you would have to do to surpass the builtin ODE solvers makes it impractical not to use them. $\endgroup$ – probably_someone Nov 7 '17 at 20:48
  • $\begingroup$ I've never used the ODE45 before but I've read up on and and even watched videos of how to do it. But I'm still not certain how I can translate my equations into something the ODE45 could understand and solve, and I'm afraid that it might not be accurate enough. $\endgroup$ – J.Doe Nov 7 '17 at 20:55
  • 3
    $\begingroup$ Would Computational Science be a better home for this question? (If so, flag for a moderator to migrate it; don't cross-post.) $\endgroup$ – Emilio Pisanty Nov 7 '17 at 21:09
  • 2
    $\begingroup$ Use the existing ODE solver packages -- you will save yourself weeks of work! $\endgroup$ – Wolfgang Bangerth Nov 8 '17 at 0:25
  • $\begingroup$ Having worked on this problem, I can tell you that you probably want to use a variable order solver such as ode113- the guiding center paths are typcially very smooth with abrupt changes where the particles "turn around." ode113() will adjust to this. $\endgroup$ – Brian Borchers Nov 9 '17 at 18:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy