python - odeint 中的参数数组

Question

我正在尝试用 odeint 求解微分方程。这里一些常量参数是固定的，一些在列表中。

    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(轨迹, w0, t, args=(p,)) ValueError:使用序列设置数组元素。

Answer 1

问题的最新版本 - 是 XY 问题，需要更加动态的电场评估

您正在尝试将 E 字段值的更大插值网格传递给 ODE 函数。这不是你需要的。而且也是不可能的，参数数组并不是用于这样的目的。 _{(这就是为什么它是一个 XY 问题，你想解决 X，使用方法 Y 并得到一个问题，然后尝试在不与 X 沟通的情况下对 Y 进行故障排除，但事实证明方法 Y 不是一个好的解决方案，你应该使用其他一些方法Z)}

ODE 函数需要当前坐标处的电场值。只需将插值器设为全局对象并在 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()