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.

$$d{v}_i/dt=q/m{F}_i$$

$$d{x}_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+2k_2+2k_3+k_4)$$

where

$$\begin{array}{r}{k}_1=f(t_n,y_n)\\ k_2=f(t_n+h/2,y_n+h/2k_1)\\ k_3=f(t_n+h/2,y_n+h/2k_2)\\ k_4=f(t_n+h,y_n+h{k}_3)\end{array}$$

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{array}{r}{k}_1=hf(t_k,y_k)\\ k_2=hf(t_k+h/4,y_k+k_1/4)\\ k_3=hf(t_k+\frac{3}{8}h,y_k+\frac{3}{32}{k}_1+\frac{9}{32}{k}_2)\\ k_4=hf(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=hf(t_k+h,y_k+\frac{439}{216}{k}_1-8k_2+\frac{3680}{513}{k}_3-\frac{845}{4104}{k}_4)\\ k_6=hf(t_k+\frac{1}{2}h,y_k-\frac{8}{27}{k}_1+2k_2-\frac{3544}{2565}{k}_3+\frac{1859}{4104}{k}_4-\frac{11}{40}{k}_5)\end{array}$$

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{array}{r}{k}_1=hf(t_k,y_k)\\ k_2=hf(t_k+h/4,y_k+k_1/4)\\ k_3=hf(t_k+\frac{3}{8}h,y_k+\frac{3}{32}{k}_1+\frac{9}{32}{k}_2)\\ k_4=hf(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=hf(t_k+h,y_k+\frac{439}{216}{k}_1-8k_2+\frac{3680}{513}{k}_3-\frac{845}{4104}{k}_4)\\ k_6=hf(t_k+\frac{1}{2}h,y_k-\frac{8}{27}{k}_1+2k_2-\frac{3544}{2565}{k}_3+\frac{1859}{4104}{k}_4-\frac{11}{40}{k}_5)\end{array}$$

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.