本文介绍了odeint 中的参数数组的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我正在尝试用 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)

仅当 qm 在集成的不同运行中可变时才应更改此设置.

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 中的参数数组的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!

05-27 16:29
查看更多