Функция odeint Python, похоже, не работает

Я хочу изучить движение заряженной частицы во время путешествия через магнитное поле, моделируя его с помощью Python. Я попытался использовать функцию odeint из scipy.integrate, и, похоже, она работает не так, как я ожидал. Вот что я ожидал, учитывая начальное условие:

Движение заряженной частицы в магнитном поле

Но вот что я получил с моей симуляцией:

Кривая не такая, какой должна быть

Проблема связана с моей реализацией функции odeint. Любая помощь приветствуется.

Вот мой код:

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import odeint

from mpl_toolkits.mplot3d import Axes3D


def vect_prod(u, v):
    return np.array([u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]])


def dy(y, t):
    x1, Vx, y1, Vy, z1, Vz = y

    F = q * (E + vect_prod(np.array([Vx, Vy, Vz]), B))

    dy = [Vx, Vx - F[0] / m, Vy, Vy - F[1] / m, Vz, Vz - F[2] / m]


    return dy


E = np.array([0, 0, 0])

B = np.array([0, 0, 1])

q = 1

m = 1

a = 0.04

cond = [0, 1, 0, 1, 0, 1]

t = np.arange(0, 100, 0.1)

sol = odeint(dy, cond, t)

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

plt.plot(sol[:, 0], sol[:, 2], sol[:, 4])

plt.show()

Любая помощь будет оценена!


person Community    schedule 03.03.2018    source источник


Ответы (1)


Я думаю, что масса слишком велика:

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import odeint

from mpl_toolkits.mplot3d import Axes3D


def vect_prod(u, v):
    return np.array([u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]])


def dy(y, t):
    x1, Vx, y1, Vy, z1, Vz = y

    F = q * (E + vect_prod(np.array([Vx, Vy, Vz]), B))

    dy = [Vx, Vx - F[0] / m, Vy, Vy - F[1] / m, Vz, Vz - F[2] / m]


    return dy


E = np.array([0, 0, 0])

B = np.array([0, 0, 1])

q = 1

m = 0.001

a = 0.04

cond = [0, 1, 0, 1, 0, 1]

t = np.arange(0, 0.05, 0.0001)

sol = odeint(dy, cond, t)

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

plt.plot(sol[:, 0], sol[:, 2], sol[:, 4])

plt.show()

выход:

введите здесь описание изображения

person HYRY    schedule 03.03.2018
comment
Я поставил ту же массу, что и вы, и, похоже, это не решило проблему; вот что я получаю: i.imgur.com/v6MauXq.png - person ; 03.03.2018
comment
Вам также нужно изменить время: t = np.arange(0, 0.05, 0.0001) - person HYRY; 03.03.2018
comment
Я проверил, и это сработало, но это было неожиданно, потому что ускорения, которые я задал в функции dy, были неправильными. - person ; 03.03.2018
comment
На самом деле, проведя несколько тестов, я обнаружил, что ваше решение работает очень недолго. Но когда я устанавливаю большую продолжительность, мы видим, что она закручивается по оси Z, что не соответствует движению в магнитном поле. imgur.com/Om7vEGM - person ; 03.03.2018