последняя версия вопроса - была проблема XY, нужна динамическая оценка E-поля
Вы пытаетесь передать большую интерполированную сетку значений поля E в функцию ODE. Это не то, что вам нужно. И также не возможно, массив параметров не был предназначен для такой цели. (Вот почему это была проблема XY, вы хотите решить X, использовать метод Y и получить проблему, а затем попробовать устранить неполадки Y без сообщения X, но оказывается, что метод Y не является хорошим решением, и вы должны использовать какой-то другой метод Z) * 1004 *
Функция ODE нуждается в значении E-поля в текущих координатах. Просто сделайте интерполятор глобальным объектом и используйте его в функции ODE, согласно документации по функции интерполяции, это должно работать. Использование заданных точек данных для заполнения углов куба с длиной стороны 4, чтение из строки вместо файлов,
griddata = """ 597.8291 0.0 172.9540 -2.0 -2.0 -2.0
561.7756 204.4696 172.9540 -2.0 -2.0 2.0
457.9636 384.2771 172.9540 -2.0 2.0 -2.0
298.9145 517.7352 172.9540 -2.0 2.0 2.0
103.8119 588.7467 172.9540 2.0 -2.0 -2.0
-103.8119 588.7467 172.9540 2.0 -2.0 2.0
-298.9145 517.7352 172.9540 2.0 2.0 -2.0
-457.9636 384.2771 172.9540 2.0 2.0 2.0""";
grid = [ [ float(cc) for cc in line.split()] for line in griddata.split("\n")];
grid = np.asarray(grid);
xyz_grid = grid[:,3:] # xyz_grid = np.array([y4, y5, y6]).T
E_grid = grid[:,:3] # E_grid = np.array([yy1, yy2, yy3]).T
E_field = LinearNDInterpolator( xyz_grid, E_grid )
def trajectory(w, t):
x, vx, y, vy, z, vz = w
Ex, Ey, Ez = E_field([x, y, z])[0] # returns list of arrays
f = [ vx, q*(Ex + vy * Bz - vz * By) / m,
vy, q*(Ey + vz * Bx - vx * Bz) / m,
vz, q*(Ez + vx * By - vy * Bx) / m ]
return f
Обратите внимание, что здесь все константы, используемые в функции, являются глобальными константами, поэтому вызов
wsol = odeint(trajectory, w0, t)
это следует изменить только в том случае, если q и m являются переменными для разных прогонов интеграции.
Вам, вероятно, следует изменить масштаб переменных положения и времени, чтобы координаты и скорости, которые видит odeint, находились в диапазоне звездных величин 0.1..10. В противном случае допуски (по умолчанию) могут привести к странным изменениям в отдельных компонентах.
старая версия вопроса, неправильно построенный вектор параметров
Оболочка lsode odeint пытается преобразовать ваш список параметров в массив. Ожидается, что этот список представляет собой плоский список чисел. Ваш список содержит другие списки, которые дают гетерогенную структуру, которая не помещается в массив numpy.
Необходимо задаться вопросом, для чего предназначены списки fEx и т. Д., Поскольку функция ODE использует эти параметры, как если бы они были числами.
from scipy.integrate import odeint
import numpy as np
m = 9.1e-31
q = 1.6e-19
Bx = 0.1825e-4
By = 0.00942e-4
Bz = 0.46264e-4
fEt = [ 0, 2e-9, 4e-9, 6e-9, 8e-9, 10e-9]
fEx = [0.20507215, 0.20658776, 0.20810338, 0.20961899, 0.21113461, 0.21265022]
fEy = [0.17207596, 0.16972669, 0.16737742, 0.16502815, 0.16267888, 0.1603296]
fEz = [ 3.90810619e-01, 3.60677316e-01, 3.30544013e-01, 3.00410711e-01, 2.70277408e-01, 2.40144105e-01 ]
def trajectory(w, t, p):
q, m = p # not really necessary, global variables work here fine
x1, y1, x2, y2, x3, y3 = w
Ex, Ey, Ez = np.interp(t,fEt, fEx), np.interp(t,fEt, fEy), np.interp(t,fEt, fEz)
f = [y1, q*(Ex + y2 * Bz - y3 * By) / m, y2, q*(Ey - y1 * Bz + y3 * Bx) / m, y3, q*(Ez + y1 * By - y2 * Bx) / m]
return f
x1, y1 = 0.0, 0.0
x2, y2 = 0.0, 0.0
x3, y3 = 0.006, 68999.35
#time
t = np.arange(0, 10, 0.01)*1e-9
p = [q, m]
w0 = [x1, y1, x2, y2, x3, y3]
# Call the ODE solver.
wsol = odeint(trajectory, w0, t, args=(p,))
print wsol
x1, y1, x2, y2, x3, y3 = wsol.T
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig=plt.figure()
ax=fig.gca(projection='3d')
ax.plot(x1,x2,x3,'r',label='charged particle trajectory')
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
ax.legend()
plt.show()
