python - 如何在python中使用文本文件求解具有多个参数数据的微分方程？

Question

我只是在学习编程。我用常数参数在python中求解了二阶微分方程，但我想用来自文本文件的参数的多个数据来求解方程。对于所有 Ex、Ey 和 Ez，我在文本文件中有 1000 多个不同的值。Ex、Ey 和 Ez 是 x、y 和 z 的函数。如何求解所有值的微分方程？我想要最后的单一情节（轨迹）。这是我的代码：

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate

m = 9.1 *(10)**(-31)    
q = 1.6 *(10)**(-19)    
Bx = 0.00942 *(10)**(-4)
By = 0.46264 *(10)**(-4)
Bz = 0.1825 *(10)**(-4)





with open("field_x.txt") as f:
    flines = f.readlines()
    Ex = [float(line.split()[0]) for line in flines] #Ex = 1.579, 3.456, 16.25, 68.99,.....

with open("field_y.txt") as f:
    flines = f.readlines()
    Ey = [float(line.split()[0]) for line in flines]  #Ey = 25.793, 23.14, 14.25, 16.25,......

with open("field_z.txt") as f:
    flines = f.readlines()
    Ez = [float(line.split()[0]) for line in flines]#Ez = -102457.243, 2003.56, 21134.89,210.35,...

#Position x,y,z
with open("position_x.txt") as f:
    flines = f.readlines()
    x = [float(line.split()[0]) for line in flines] #x=0.1, 0.3, 0.5,...

with open("position_y.txt") as f:
    flines = f.readlines()
    y = [float(line.split()[0]) for line in flines] #y= 0.0, 0.5, 0.6,...

with open("position_z.txt") as f:
    flines = f.readlines()
    z = [float(line.split()[0]) for line in flines] #z=0.1, 0.5, 0.64


#array for field components
Ex1 = np.array(Ex)
Ey1 = np.array(Ey)
Ez1 = np.array(Ez)



#array for position
xs = np.array(x) 
ys = np.array(y)
zs = np.array(z)

#linear interpolation of Ex1, Ey1, Ez1
Exf = LinearNDInterpolator((xs, ys, zs), Ex1)
Eyf = LinearNDInterpolator((xs, ys, zs), Ey1)
Ezf = LinearNDInterpolator((xs, ys, zs), Ez1)

#array of new point
x1 = np.linspace(0, 5, 10)
y1 = np.linspace(0, 7, 10)
z1 = np.linspace(0, 10, 10)

#creating array([x1,y1,z1],[x2,y2,z2],....) for new grids
X = np.dstack((x1,y1,z1))
points = np.array(X)

#Field at new grids after linear interpolation
fEx = Exf(points)
fEy = Eyf(points)
fEz = Ezf(points)

def trajectory(w, t, p):

    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] 

#initial condition
x1 = 0.0
y1 = 0.0
x2 = 0.0
y2 = 0.0
x3 = 0.006
y3 = 691762.096

#time
t = np.arange(0*(10)**(-9), 10.0*(10)**(-9), 0.01*(10)**(-9))


p = [q, m, fEx, fEy, fEz, Bx, By, Bz]
w0 = [x1, y1, x2, y2, x3, y3]

#SOlving the equation
wsol = odeint(trajectory, w0, t, args=(p,))

#output data
X = wsol[:,0]      
Y = wsol[:,2]      
Z = wsol[:,4]      

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

我没有从这个解决方案中得到正确的输出。我期待最后的单一情节。

标签： python