irfpy.util.rungekutta

Runge-Kutta implementations (for particle tracing)

For tutorial of particle tracing, see Tracing particles (RKF45).

In a general term, Runge-Kutta solves

\[dy/dt = f(t, y)\]

The functions

• rk4_progress()

• rk5_progress()

• rkf45_progress()

will integrate the above generic differential equation by 1 step.

Classes

• RK_EqM

• RK4_EqM

• RK5_EqM

• RKF45_EqM

• RKF45_BackEqM

provide the forward integrator of the equation of motion.

\[dv_i/dt = q/m F_i\]

\[dx_i/dt = v_i\]

class irfpy.util.rungekutta.ElectricForce(mass, charge, electric_field)

Bases: object

Parameters:

• mass – The mass of the particle, in kg.

• charge – The charge of the particle, in C.

• vector3dfield – An instance of irfpy.util.fields.Vector3dField representing the electric field (E).

class irfpy.util.rungekutta.MagneticForce(mass, charge, magnetic_field)

Bases: object

Parameters:

• mass – The mass of the particle, in kg.

• charge – The charge of the particle, in C.

• vector3dfield – An instance of irfpy.util.fields.Vector3dField representing the magnetic field (B).

>>> m = 1
>>> q = 1
>>> import irfpy.util.fields
>>> bfield = irfpy.util.fields.UniformVector3dField(0, 0, 1)    #  Bz=1
>>> vel = [1, 0, 0]    # Vx=1
>>> forcefield = MagneticForce(m, q, bfield)
>>> f = forcefield(0, 0, 0, 0, vel[0], vel[1], vel[2])
>>> print(f)
[ 0. -1.  0.]

# A bit more complex example >>> bfield = irfpy.util.fields.UniformVector3dField(-1, 2, 6) >>> vel = [-2, -3, -4] >>> forcefield = MagneticForce(m, q, bfield) >>> f = forcefield(0, 0, 0, 0, vel[0], vel[1], vel[2]) >>> print(f) [-10. 16. -7.]

class irfpy.util.rungekutta.ElectroMagneticForce(mass, charge, electric_field, magnetic_field)

Bases: object

class irfpy.util.rungekutta.GravityForce(center_mass)

Bases: object

Gravity force class to be used for RK.

Parameters:

center_mass – Mass of the central object (or equivalent mass) in kg.

gravity_constant = 6.674e-11

class irfpy.util.rungekutta.RK_EqM(x0, v0, forcefield)

Bases: object

Classical RK.

\[y_{n+1} = y_n + \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4)\]

where

\[\begin{split}k_1 = f(t_n, y_n) \\ k_2 = f(t_n + h/2, y_n + h/2 k_1) \\ k_3 = f(t_n + h/2, y_n + h/2 k_2) \\ k_4 = f(t_n + h, y_n + h k_3)\end{split}\]

Example

Here we try to understand the Io motion around Jupiter.

The mass of Juptier is 1.90e27 kg and Io for 8.93e22.

>>> massj = 1.90e27

Io is at 4.217e5 km with velocity of 17.334 km/s at t=0.

>>> x0 = np.array([4.217e8, 0, 0])
>>> v0 = np.array([0, 17.334e3, 0])

The force per mass (N/kg) is define as a GravityForce class.

>>> gforce = GravityForce(massj)

The RK_EqM class is instanced then.

>>> rk = RK_EqM(x0, v0, gforce)

Let’s put the time resolution as 3600 s.

>>> rk.progress(3600)
>>> print(rk.pos[-1])   
[4.17087687e+08   6.21744701e+07   0.00000000e+00]
>>> print(rk.vel[-1])   
[-2557.7194899  17144.4089548      0.       ]
>>> print(rk.t[-1])
3600

After 21 hours, the Io goes to the other side (-x with -vy).

>>> for i in range(20):
...     rk.progress(3600)
>>> print(rk.pos[-1])   
[-4.20854860e+08   1.22957706e+07   0.00000000e+00]
>>> print(rk.vel[-1])   
[  -506.56403226 -17353.95549486      0.        ]
>>> print(rk.t[-1])
75600

Parameters:

• x0 – Initial position, as numpy array with the shape of (3,)

• v0 – Initial velocity, as numpy array with the shape of (3,)

• forcefield – Force field acting to the particle per mass. This should be a function(-like) that must get 7 floats (t, x, y, z, vx, vy, vz) and return the 3 floats (fx, fy, fz).

progress(dt)

Progress the particle by dt.

irfpy.util.rungekutta.rk4_progress(tn, yk, dt, dfunc)

Integration by RK4 (6 order).

class irfpy.util.rungekutta.RK4_EqM(x0, v0, forcefield)

Bases: object

Runge-Kutta with the order of 4 using 6 k-s.

\[y_{k+1} = y_k + \frac{25}{216}k_1 + \frac{1408}{2565}k_3 + \frac{2197}{4101}k_4 - \frac{1}{5}k_5\]

where

\[\begin{split}k_1 = h f(t_k, y_k) \\ k_2 = h f(t_k + h/4, y_k + k_1/4) \\ k_3 = h f(t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) \\ k_4 = h f(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1 - \frac{7200}{2197}k_2 + \frac{7296}{2197} k_3) \\ k_5 = h f(t_k + h, y_k + \frac{439}{216}k_1 - 8 k_2 + \frac{3680}{513} k_3 - \frac{845}{4104} k_4) \\ k_6 = h f(t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2 k_2 - \frac{3544}{2565} k_3 + \frac{1859}{4104} k_4 - \frac{11}{40} k_5) \\\end{split}\]

Example

Here we try to understand the Io motion around Jupiter.

The mass of Juptier is 1.90e27 kg and Io for 8.93e22.

>>> massj = 1.90e27

Io is at 4.217e5 km with velocity of 17.334 km/s at t=0.

>>> x0 = np.array([4.217e8, 0, 0])
>>> v0 = np.array([0, 17.334e3, 0])

The force per mass (N/kg) is define as a GravityForce class.

>>> gforce = GravityForce(massj)

The RK4_EqM class is instanced then.

>>> rk = RK4_EqM(x0, v0, gforce)

Let’s put the time resolution as 3600 s.

>>> rk.progress(3600)
>>> print(rk.pos[-1])   
[4.17084349e+08   6.21989712e+07   0.00000000e+00]
>>> print(rk.vel[-1])   
[-2558.70808725  17144.27369938      0.        ]
>>> print(rk.t[-1])
3600

After 21 hours, the Io goes to the other side (-x with -vy).

>>> for i in range(20):
...     rk.progress(3600)
>>> print(rk.pos[-1])   
[-4.20899082e+08   1.18280482e+07   0.00000000e+00]
>>> print(rk.vel[-1])
[  -487.75930547 -17353.750109        0.        ]
>>> print(rk.t[-1])
75600

Parameters:

• x0 – Initial position, as numpy array with the shape of (3,)

• v0 – Initial velocity, as numpy array with the shape of (3,)

• forcefield – Force field acting to the particle per mass. This should be a function(-like) that must get 7 floats (t, x, y, z, vx, vy, vz) and return the 3 floats (fx, fy, fz).

progress(dt)

Progress the particle by dt.

irfpy.util.rungekutta.rk5_progress(tn, yk, dt, dfunc)

class irfpy.util.rungekutta.RK5_EqM(x0, v0, forcefield)

Bases: object

Runge-Kutta with the order of 4 using 6 k-s.

\[y_{k+1} = y_k + \frac{16}{135} k_1 + \frac{6656}{12825} k_3 + \frac{28561}{56430} k_4 - \frac{9}{50} k_5 + \frac{2}{55} k_6\]

where

\[\begin{split}k_1 = h f(t_k, y_k) \\ k_2 = h f(t_k + h/4, y_k + k_1/4) \\ k_3 = h f(t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) \\ k_4 = h f(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1 - \frac{7200}{2197}k_2 + \frac{7296}{2197} k_3) \\ k_5 = h f(t_k + h, y_k + \frac{439}{216}k_1 - 8 k_2 + \frac{3680}{513} k_3 - \frac{845}{4104} k_4) \\ k_6 = h f(t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2 k_2 - \frac{3544}{2565} k_3 + \frac{1859}{4104} k_4 - \frac{11}{40} k_5) \\\end{split}\]

Example

Here we try to understand the Io motion around Jupiter.

The mass of Juptier is 1.90e27 kg and Io for 8.93e22.

>>> massj = 1.90e27

Io is at 4.217e5 km with velocity of 17.334 km/s at t=0.

>>> x0 = np.array([4.217e8, 0, 0])
>>> v0 = np.array([0, 17.334e3, 0])

The force per mass (N/kg) is define as a GravityForce class.

>>> gforce = GravityForce(massj)

The RK5_EqM class is instanced then.

>>> rk = RK5_EqM(x0, v0, gforce)

Let’s put the time resolution as 3600 s.

>>> rk.progress(3600)
>>> print(rk.pos[-1])   
[4.17087705e+08   6.21747255e+07   0.00000000e+00]
>>> print(rk.vel[-1])   
[-2557.71081377  17144.41010741      0.        ]
>>> print(rk.t[-1])
3600

After 21 hours, the Io goes to the other side (-x with -vy).

>>> for i in range(20):
...     rk.progress(3600)
>>> print(rk.pos[-1])   
[-4.20862773e+08   1.23126860e+07   0.00000000e+00]
>>> print(rk.vel[-1])
[  -507.29360568 -17353.63323363      0.        ]
>>> print(rk.t[-1])
75600

Parameters:

• x0 – Initial position, as numpy array with the shape of (3,)

• v0 – Initial velocity, as numpy array with the shape of (3,)

• forcefield – Force field acting to the particle per mass. This should be a function(-like) that must get 7 floats (t, x, y, z, vx, vy, vz) and return the 3 floats (fx, fy, fz).

progress(dt)

Progress the particle by dt.

irfpy.util.rungekutta.rk45_progress(tn, yk, h, dfunc)

4- and 5-dimensional Runge-Kutta solver.

The following two calss are equivalent.

rk4, rk5 = rk45_progress(tn, yk, h, dfunc)
rk4, rk5 = rk4_progress(tn, yk, h, dfunc), rk5_progress(tn, yk, h, dfunc)

But for performance reasons, the results of rk4_progress() and rk5_rk5_progress() is calculaated & returned together by this function.

Parameters:

• tn – The time step to be fed into dfunc

• yk – The y coordinate, an (N,) shaped array

• h – Delta to proceed. Sclara.

• dfunc – Differential function. It should get the argument (t, y0, y1, … yN-1), returning the (N,) shaped np-array.

Returns:

Two np-array. 4Dim and 5Dim.

irfpy.util.rungekutta.rkf45_progress(tn, yk, dt, dfunc, dtmin=None, dtmax=None, tol=1e-05, reltol=1e-08)

Single step for RKF.

Parameters:

• tn – The time. To be fed to dfunc.

• yk – The y-value at time tn.

• dt – Time step to evaluate first. The exact time step is defined automatically.

• dfunc – Function to be integrated. This should eat tn and y at tn (array), returning the dy/dt value (array).

• dtmin – The minimum time step for integration. Default is dt/1000.

• dtmax – The maximum time step for integration. Default is dt * 1000.

• tol – Absolute tolerance. Scalar or array with the same shape of yk.

• reltol – Relative tolerance. Scalar or array with the same shape of yk.

Algorithm

• Obtain rk4 = rk4_progress() and rk5 = rk5_progress().

• Error is E=|rk5 - rk4|

• Tolerance is the bigger one out of tol and reltol * (rk5 - yk).

• Compare error and tolerance.

  • If error > tolerance: tn is halven, then recursive call. (Too bad approximation)

  • If error / tolerance is between 0.2 and 1: rk5 is the progressed value.

  • If error < 0.2 * tolerance: tn is increased by 1.5 times, then recursive call. (Too good approximation)

class irfpy.util.rungekutta.RKF45_EqM(x0, v0, forcefield, dt0, dtmin=None, dtmax=None, tolx=1e-05, tolv=1e-05, reltolx=1e-08, reltolv=1e-08)

Bases: object

Solve equation of motion by RKF.

>>> from irfpy.util import fields
>>> efield = fields.UniformVector3dField(0, 1e-6, 0)
>>> eforce = ElectricForce(1.67e-27, 1.6e-19, efield)

>>> x0 = np.array([0, 0, 0])
>>> v0 = np.array([1, 0, 0])
>>> rk = RKF45_EqM(x0, v0, eforce, 30e-3, tolx=0.01, tolv=0.001)
>>> rk.progress()
>>> rk.progress()

Parameters:

• dtmin – Minimum h allowed. Default is 0.01*dt0.

• dtmax – Maximum h allowed. Default is 100*dt0.

progress()

Progress the particle by dt.

class irfpy.util.rungekutta.RKF45_BackEqM(x0, v0, forcefield, dt0, dtmin=None, dtmax=None, tolx=1e-05, tolv=1e-05, reltolx=1e-08, reltolv=1e-08)

Bases: object

Backward integration of the equation of motion.

>>> from irfpy.util import fields
>>> efield = fields.UniformVector3dField(0, 1e-6, 0)
>>> eforce = ElectricForce(1.67e-27, 1.6e-19, efield)

>>> x0 = np.array([0, 0, 0])
>>> v0 = np.array([1, 0, 0])
>>> rk = RKF45_BackEqM(x0, v0, eforce, 30e-3, tolx=0.000001, tolv=0.00001)
>>> rk.retrogress()
>>> rk.retrogress()
>>> rk.retrogress()

Parameters:

• dt0 – The time step, must be positive!

• dtmin – Minimum h allowed. Default is 0.001*dt0.

• dtmax – Maximum h allowed. Default is 1000*dt0.

retrogress()

Progress the particle by dt.