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