Search code examples
pythonscipyodeodeint

Array of parameter in 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()

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.

Solution

  • recent version of question - was XY problem, need E-field evaluation more dynamical

    You are trying to pass a larger interpolated grid of values of the E field to the ODE function. This is not what you need. And also not possible, the parameter array was not intended for such a purpose. (That's why it was an XY problem, you want to solve X, use method Y and get a problem and then try troubleshooting Y without communicating X, but it turns out that method Y is not a good solution and you should use some other method Z)

    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)

    this should only be changed if q and m are variable over different runs of the integration.

    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.

    old version of question, wrongly constructed parameter vector

    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.

    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()

    3D plot of solution