I am just learning programming. I solved 2nd order differential equation in python with constant parameters but I want to solve the equation with multiple data for the parameters from text file. I have more than 1000 different values in the text file for all Ex, Ey and Ez. Ex, Ey and Ez are the function of x, y and z. How to solve the differential equation with all the values? I want single plot (trajectory) at the end. Here is my code:
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()
I am not getting the right output from this solution. I am expecting single plot at the end.
