Блог

Я хочу численно решить уравнения движения релятивистской частицы в поле E/M:

где, , и,
.

Я начал с упрощения с Симпатией, чтобы решить первое уравнение для . С помощью

 """ Define symbols"""
t,m,q, c = smp.symbols('t m q c', Real =True)
x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz, gamma = smp.symbols('x y z vx vy vz  E_x E_y E_z B_x B_y B_z gamma ',cls=smp.Function, real=True)

""" Define functions"""

x    = x(t);                y   = y(t);               z    = z(t);                 # Position   
vx   = vx(t);              vy   = vy(t);              vz   = vz(t);                # Velocity  
Ex   = Ex(x, y, z, t);     Ey   = Ey(x, y, z, t);     Ez   = Ez(x, y, z, t);       # Components of E field 
Bx   = Bx(x, y, z, t);     By   = By(x, y, z, t);     Bz   = Bz(x, y, z, t);       # Components of B field
x_d  = smp.diff(x,t);      y_d  = smp.diff(y,t);      z_d  = smp.diff(z,t);        # dx_{i}/dt ,    x_{i} = x, y, z
vx   = x_d;                vy   = y_d;                vz   = z_d;                  # v_{i} ,           i  = x, y, z 
vx_d = smp.diff(x,t,2);    vy_d = smp.diff(y,t,2);    vz_d = smp.diff(z,t,2);      # d^{2}x_{i}/dt^2 , i  = x, y, z 

""" Define Matrices"""
r      = smp.Matrix([x, y, z])
E      = smp.Matrix([Ex, Ey, Ez])
B      = smp.Matrix([Bx, By, Bz])
v      = smp.Matrix([vx, vy, vz])

""" Define gamma """
v_norm = v.norm()
gamma  = gamma(v_norm)
gamma  = 1/smp.sqrt(1 - (v.norm()/c)**2)

""" Define momentum """
p    = m*gamma*v

""" 2nd order ODE of motion """
de1 = smp.diff(p,t)-q*(E  v.cross(B)/c)

""" Solve for dv_{i}/dt, i = x, y, z """ 
sols   = smp.solve([de1[0],de1[1], de1[2]], (vx_d, vy_d, vz_d), Simplify =False, rational=False)

"""Create numpy functions that we can use with numerical methods"""
list_p = (c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz)


dvx_dt_f  = smp.lambdify(list_p,sols[vx_d], modules=smp)
dx_dt_f   = smp.lambdify(vx,vx, modules=smp)

dvy_dt_f  = smp.lambdify(list_p,sols[vy_d], modules=smp)
dy_dt_f   = smp.lambdify(vy, vy, modules=smp)

dvz_dt_f  = smp.lambdify(list_p,sols[vz_d], modules=smp)
dz_dt_f   = smp.lambdify(vz, vz, modules=smp)
 

Затем я хочу использовать модифицированные выражения и интегрировать их с помощью ivp Scipy в следующей форме.

 def dSdt(S,t):
    vx, x, vy, y, vz, z = S
    return[dvx_dt_f(c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz),
           dx_dt_f(vx),
           dvy_dt_f(c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz),
           dy_dt_f(vy),
           dvz_dt_f(c, m, q, x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz),
           dz_dt_f(vz)
          ]
 

для начальных условий v0 = (x, y, z, vx, vy, vz, Ex, Ey, Ez, Bx, By, Bz)

Однако вместо того, чтобы получить решение в терминах Vx, Vy, Vz, необходимое для численного интегрирования, я получаю его в терминах производных по времени $frac{dx}{dt}.

Что я могу сделать, чтобы изменить производные с помощью $v_x, v_y, v_z$?