- #1
goldin2008
- 6
- 0
I want to solve a set of equations using Python odeint, but output shows me it is wrong.
Can you help me? Thanks.
Code:
# -*- coding: utf-8 -*-
from scipy.integrate import odeint
import numpy as np
from pylab import *
import math
def func(y, t, k, c, Zr):
#px, py, pz, x, y, z = y & px=y[0]; py=y[1]; pz=y[2]; x=y[3]; y=y[4]; z=y[5]
# apply func formula
return (
-k/((1+y[5]**2/Zr**2)**(1.0/2))*exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))
-((-1.0/2*(1+y[5]**2/Zr**2)**(-3.0/2)*2*y[5]/Zr**2*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+(1+y[5]**2/Zr**2)**(-1.0/2)*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(1+y[5]**2/Zr**2)**(-2.0))*2*y[5]/Zr**2*(y[3]**2+y[4]**2)/(2*Zr/k)*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*(k+(y[3]**2+y[4]**2)*(Zr**2-y[5]**2)/(2.0/k*(Zr**2+y[5]**2))-1.0/(1+y[5]**2/Zr**2)))-1.0/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+
y[3]*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[3])/(2*Zr/k*(1+y[5]**2/Zr**2))*
(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*2*y[3]*y[5]/(2/k*(Zr**2+y[5]**2))+y[5]/Zr*2*y[3]*y[4]/(2.0/k*(Zr**2+y[5]**2))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))))))*y[2]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
0+(y[3]/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[4]/(2*Zr/k*(1+y[5]**2/Zr**2))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+2*y[4]*y[5]/(2/k*(Zr**2+y[5]**2))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+y[5]/Zr*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))))*y[2]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
(y[3]/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(k*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+k*y[5]/Zr*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))))
+(y[0]*((-1.0/2*(1+y[5]**2/Zr**2)**(-3.0/2)*2*y[5]/Zr**2*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+(1+y[5]**2/Zr**2)**(-1.0/2)*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(1+y[5]**2/Zr**2)**(-2.0))*2*y[5]/Zr**2*(y[3]**2+y[4]**2)/(2*Zr/k)*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*(k+(y[3]**2+y[4]**2)*(Zr**2-y[5]**2)/(2.0/k*(Zr**2+y[5]**2))-1.0/(1+y[5]**2/Zr**2)))-1.0/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+
y[3]*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[3])/(2*Zr/k*(1+y[5]**2/Zr**2))*
(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*2*y[3]*y[5]/(2/k*(Zr**2+y[5]**2))+y[5]/Zr*2*y[3]*y[4]/(2.0/k*(Zr**2+y[5]**2))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))))))-
y[1]*((y[3]/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[4]/(2*Zr/k*(1+y[5]**2/Zr**2))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+2*y[4]*y[5]/(2/k*(Zr**2+y[5]**2))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+y[5]/Zr*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))))))/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
y[0]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
y[1]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
y[2]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2))
y0=[0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
k=30.0
Zr=120.0
c=3e8
t = np.arange(0.0, 1000.0, 0.05) # create time point
# use odeint to solve func
track= odeint(func, y0, t, args=(k, c, Zr))
# plot
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track[:,3], track[:,4], track[:,5])
plt.show()
Can you help me? Thanks.
Code:
# -*- coding: utf-8 -*-
from scipy.integrate import odeint
import numpy as np
from pylab import *
import math
def func(y, t, k, c, Zr):
#px, py, pz, x, y, z = y & px=y[0]; py=y[1]; pz=y[2]; x=y[3]; y=y[4]; z=y[5]
# apply func formula
return (
-k/((1+y[5]**2/Zr**2)**(1.0/2))*exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))
-((-1.0/2*(1+y[5]**2/Zr**2)**(-3.0/2)*2*y[5]/Zr**2*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+(1+y[5]**2/Zr**2)**(-1.0/2)*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(1+y[5]**2/Zr**2)**(-2.0))*2*y[5]/Zr**2*(y[3]**2+y[4]**2)/(2*Zr/k)*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*(k+(y[3]**2+y[4]**2)*(Zr**2-y[5]**2)/(2.0/k*(Zr**2+y[5]**2))-1.0/(1+y[5]**2/Zr**2)))-1.0/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+
y[3]*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[3])/(2*Zr/k*(1+y[5]**2/Zr**2))*
(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*2*y[3]*y[5]/(2/k*(Zr**2+y[5]**2))+y[5]/Zr*2*y[3]*y[4]/(2.0/k*(Zr**2+y[5]**2))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))))))*y[2]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
0+(y[3]/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[4]/(2*Zr/k*(1+y[5]**2/Zr**2))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+2*y[4]*y[5]/(2/k*(Zr**2+y[5]**2))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+y[5]/Zr*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))))*y[2]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
(y[3]/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(k*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+k*y[5]/Zr*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))))
+(y[0]*((-1.0/2*(1+y[5]**2/Zr**2)**(-3.0/2)*2*y[5]/Zr**2*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+(1+y[5]**2/Zr**2)**(-1.0/2)*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(1+y[5]**2/Zr**2)**(-2.0))*2*y[5]/Zr**2*(y[3]**2+y[4]**2)/(2*Zr/k)*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*(k+(y[3]**2+y[4]**2)*(Zr**2-y[5]**2)/(2.0/k*(Zr**2+y[5]**2))-1.0/(1+y[5]**2/Zr**2)))-1.0/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+
y[3]*(math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[3])/(2*Zr/k*(1+y[5]**2/Zr**2))*
(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))*2*y[3]*y[5]/(2/k*(Zr**2+y[5]**2))+y[5]/Zr*2*y[3]*y[4]/(2.0/k*(Zr**2+y[5]**2))*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))))))-
y[1]*((y[3]/(Zr*(1+y[5]**2/Zr**2)**(3.0/2))*math.exp(-(y[3]**2+y[4]**2)/(2*Zr/k*(1+y[5]**2/Zr**2)))*(-2*y[4]/(2*Zr/k*(1+y[5]**2/Zr**2))*(math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))-y[5]/Zr*math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))+2*y[4]*y[5]/(2/k*(Zr**2+y[5]**2))*(math.cos(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr))+y[5]/Zr*math.sin(k*(y[5]-c*t)+y[5]*(y[3]**2+y[4]**2)/(2/k*(Zr**2+y[5]**2))-math.atan(y[5]/Zr)))))))/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
y[0]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
y[1]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2),
y[2]/(1+y[0]**2+y[1]**2+y[2]**2)**(1.0/2))
y0=[0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
k=30.0
Zr=120.0
c=3e8
t = np.arange(0.0, 1000.0, 0.05) # create time point
# use odeint to solve func
track= odeint(func, y0, t, args=(k, c, Zr))
# plot
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track[:,3], track[:,4], track[:,5])
plt.show()