问题描述
我正在尝试用 odeint 求解微分方程.这里一些常量参数是固定的,一些在列表中.
I am trying to solve the differential equation with odeint. Here some constant parameters are fixed and some are in a list.
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import LinearNDInterpolator
#equation of motion in the direction of x ===== ### d^2x/dt^2 = q[Ex + dy/dt * Bz - dz/dt * By]/m
#equation of motion in the direction of y ===== ### d^2y/dt^2 = q[Ey - dx/dt * Bz + dz/dt * Bx]/m
#equation of motion in the direction of z ===== ### d^2z/dt^2 = q[Ez + dx/dt * By - dy/dt * Bx]/m
m = 9.1 *(10)**(-31)
q = 1.6 *(10)**(-19)
#Electric field from FEMM
with open("Elecric_field_x.txt") as f:
flines = f.readlines()
yy1 = [float(line.split()[0]) for line in flines]
with open("Elecric_field_y.txt") as f:
flines = f.readlines()
yy2 = [float(line.split()[0]) for line in flines]
with open("Elecric_field_z.txt") as f:
flines = f.readlines()
yy3 = [float(line.split()[0]) for line in flines]
#Position x,y,z from FEMM
with open("Electric_position_x.txt") as f:
flines = f.readlines()
y4 = [float(line.split()[0]) for line in flines]
with open("Electric_position_y.txt") as f:
flines = f.readlines()
y5 = [float(line.split()[0]) for line in flines]
with open("Electric_position_z.txt") as f:
flines = f.readlines()
y6 = [float(line.split()[0]) for line in flines]
#data sample from FEMM inside the text file
yy1 yy2 yy3 y4 y5 y6
2.677026732329115255e-01 0.000000000000000000e+00 3.908106187718196067e-01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
1.639206109489374516e-17 2.677026732329115255e-01 3.908106187718196067e-01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
-2.677026732329115255e-01 3.278412218978749032e-17 3.908106187718196067e-01 -0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
4.031888048269389202e+01 0.000000000000000000e+00 -1.452685819581209046e+02 5.000000000000000278e-02 0.000000000000000000e+00 0.000000000000000000e+00
2.468819396416788133e-15 4.031888048269389202e+01 -1.452685819581209046e+02 3.061616997868383172e-18 5.000000000000000278e-02 0.000000000000000000e+00
-4.031888048269389202e+01 4.937638792833576266e-15 -1.452685819581209046e+02 -5.000000000000000278e-02 6.123233995736766344e-18 0.000000000000000000e+00
-2.020413445543617001e+02 -0.000000000000000000e+00 -2.380940300071312777e+03 1.000000000000000056e-01 0.000000000000000000e+00 0.000000000000000000e+00
-1.237146429519632942e-14 -2.020413445543617001e+02 -2.380940300071312777e+03 6.123233995736766344e-18 1.000000000000000056e-01 0.000000000000000000e+00
2.020413445543617001e+02 -2.474292859039265884e-14 -2.380940300071312777e+03 -1.000000000000000056e-01 1.224646799147353269e-17 0.000000000000000000e+00
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.549999999999999989e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.549999999999999989e-01
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -0.000000000000000000e+00 0.000000000000000000e+00 1.549999999999999989e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 5.000000000000000278e-02 0.000000000000000000e+00 1.549999999999999989e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 3.061616997868383172e-18 5.000000000000000278e-02 1.549999999999999989e-01
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -5.000000000000000278e-02 6.123233995736766344e-18 1.549999999999999989e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.000000000000000056e-01 0.000000000000000000e+00 1.549999999999999989e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 6.123233995736766344e-18 1.000000000000000056e-01 1.549999999999999989e-01
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -1.000000000000000056e-01 1.224646799147353269e-17 1.549999999999999989e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 3.099999999999999978e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 3.099999999999999978e-01
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -0.000000000000000000e+00 0.000000000000000000e+00 3.099999999999999978e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 5.000000000000000278e-02 0.000000000000000000e+00 3.099999999999999978e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 3.061616997868383172e-18 5.000000000000000278e-02 3.099999999999999978e-01
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -5.000000000000000278e-02 6.123233995736766344e-18 3.099999999999999978e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.000000000000000056e-01 0.000000000000000000e+00 3.099999999999999978e-01
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 6.123233995736766344e-18 1.000000000000000056e-01 3.099999999999999978e-01
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -1.000000000000000056e-01 1.224646799147353269e-17 3.099999999999999978e-01
#array for electric field components
Ex1 = np.array(yy1, dtype=object)
Ey1 = np.array(yy2, dtype=object)
Ez1 = np.array(yy3, dtype=object)
#array for position
x = np.array(y4, dtype=object)
y = np.array(y5, dtype=object)
z = np.array(y6, dtype=object)
def fE(x,y,z,yy1,yy2,yy3,y4,y5,y6):
#array for electric field components
Ex1 = np.array(yy1, dtype=object)
Ey1 = np.array(yy2, dtype=object)
Ez1 = np.array(yy3, dtype=object)
#array for position
x = np.array(y4, dtype=object)
y = np.array(y5, dtype=object)
z = np.array(y6, dtype=object)
#linear interpolation of electric field
ex = LinearNDInterpolator((x, y, z), Ex1)
ey = LinearNDInterpolator((x, y, z), Ey1)
ez = LinearNDInterpolator((x, y, z), Ez1)
#array of new point
x1 = np.linspace(0, 31, 100)
y1 = np.linspace(0, 10, 100)
z1 = np.linspace(0, 10, 100)
#creating array([x1,y1,z1],[x2,y2,z2],....) for new grids
X = np.dstack((x1,y1,z1))
points = np.array(X)
#Electric field at new grids after linear interpolation
fEx = ex(points)
fEy = ey(points)
fEz = ez(points)
return fEx, fEy, fEz
fEx, fEy, fEz = fE(x,y,z,yy1,yy2,yy3,y4,y5,y6)
#Magnetic field
Bx = 0.1825 *(10)**(-4)
By = 0.00942 *(10)**(-4)
Bz = 0.46264 *(10)**(-4)
def trajectory(w, t, p):
###====Cartesian coordinate system=====#####
#x = x1
#x_prime = y1 #dx/dt
#y = x2
#y_prime = y2 #dy/dt
#z = x3
#z_prime = y3 #dz/dt
x1, y1, x2, y2, x3, y3 = w
q, m, fEx, fEy, fEz, Bx, By, Bz = p
f = [y1, q*(fEx + y2 * Bz - y3 * By) / m, y2, q*(fEy - y1 * Bz + y3 * Bx) / m, y3, q*(fEz + y1 * By - y2 * Bx) / m] #with magnetic field
return f
#Initial conditions
x1 = 0.0
y1 = 0.0
x2 = 0.0
y2 = 0.0
x3 = 0.006
y3 = 68999.35
#time
t = np.linspace(0*(10)**(-9), 10.0*(10)**(-9), 100)
p = [q, m, fEx, fEy, fEz, Bx, By, Bz]
w0 = [x1, y1, x2, y2, x3, y3]
# Call the ODE solver.
wsol = odeint(trajectory, w0, t, args=(p,))
X = wsol[:,0] #for x
Y = wsol[:,2] #for y
Z = wsol[:,4] #for z
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t,X, color= 'b', label=('x'))
ax.plot(t,Y, color= 'r', label=('y'))
ax.plot(t,Z, color= 'c', label=('z'))
ax.set_xlabel('Time(ns)')
ax.set_ylabel('position(m)')
plt.show()
但我收到以下错误:回溯(最近一次调用最后一次):文件trajectory_cartesian.py",第 205 行,在wsol = odeint(trajectory, w0, t, args=(p,))ValueError: 使用序列设置数组元素.
But I am getting following error: Traceback (most recent call last): File "trajectory_cartesian.py", line 205, in wsol = odeint(trajectory, w0, t, args=(p,)) ValueError: setting an array element with a sequence.
推荐答案
最近版本的问题 - 是 XY 问题,需要更动态的 E-field 评估
您正在尝试将 E 字段值的较大内插网格传递给 ODE 函数.这不是你需要的.而且也不可能,参数数组不是为了这样的目的.
ODE 函数需要当前坐标处的电场值.只需将插值器设为全局对象并在 ODE 函数中使用它,根据插值函数的文档,这应该可以工作.使用给定的数据点填充边长为 4 的立方体的角,从字符串而不是文件中读取,
The ODE function needs the value of the E-field at the current coordinates. Just make the interpolator a global object and use it in the ODE function, per the documentation of the interpolation function this should work. Using the given data points to fill the corners of a cube of sidelength 4, reading from a string instead of files,
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
注意这里函数中用到的所有常量都是全局常量,所以调用的是
Note that here all constants used in the function are global constants, so the call is
wsol = odeint(trajectory, w0, t)
仅当 q
和 m
在集成的不同运行中可变时才应更改此设置.
this should only be changed if q
and m
are variable over different runs of the integration.
您可能应该重新调整位置和时间变量,以便 odeint
看到的坐标和速度都在幅度范围 0.1..10
内.否则(默认)公差可能会在单个组件中产生奇怪的变化.
You probably should rescale the position and time variables so that both the coordinates and the velocities that odeint
sees are in the magnitude range 0.1..10
. Otherwise the (default) tolerances might produce strange variations in the single components.
lsode
包装器 odeint
尝试将您的参数列表转换为数组.它期望这个列表是一个平面的数字列表.您的列表包含其他列表,这些列表提供了不适合 numpy
数组的异构结构.
The lsode
wrapper odeint
tries to cast your parameter list into an array. It expects that this list is a flat list of numbers. Your list contains other lists, which gives a heterogeneous structure that does not fit into a numpy
array.
人们不得不质疑 fEx
等列表的目的是什么,因为 ODE 函数使用这些参数就好像它们是数字一样.
One has to question what the lists fEx
etc. are meant to achieve, as the ODE function uses these parameters as if they were numbers.
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()
这篇关于odeint 中的参数数组的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!