Python ode problem (need help)

  • Python
  • Thread starter goldin2008
  • Start date
  • #1

Main Question or Discussion Point

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

Answers and Replies

  • #2
33,270
4,966
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()
You're kidding, right? You don't really expect that we're going to look at that huge wad of code, do you?
 
  • #3
You're kidding, right? You don't really expect that we're going to look at that enormous lump of goop, do you?
Sorry. Actually, just the equations are a little long.

The idea is simple.

# -*- 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()
 
  • #4
different colors mean 6 equations for y[0], y[1], y[2], y[3], y[4], y[5]
I just want to integrate them and get px, py, pz, x, y, and z. And then use x, y, and z to
draw trajectory of the particle.




# -*- 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 & 的d(px)/dt=y[0]; d(py)/dt=y[1]; d(pz)/dt=y[2]; d(x)/dt=y[3]; d(y)/dt=y[4]; d(z)/dt=y[5]
# apply func formula

return (

d(px)/dt,

d(py)/dt,

d(pz)/dt,

d(x)/dt,
d(y)/dt,
d(z)/dt)

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()
 
Last edited:
  • #5
Pythagorean
Gold Member
4,191
255
And the error message?
 
  • #6
And the error message?
I want to print track[:,3] (which is x coordinate)
but the output is

lsoda-- at current t (=r1), mxstep (=i1) steps
taken on this call before reaching tout
In above message, I1 = 500
In above message, R1 = 0.4578094279642E-05
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
[ 0.00000000e+00 -9.26095625e-12 0.00000000e+00 ..., 0.00000000e+00
0.00000000e+00 0.00000000e+00]


And the 3D trajectory is obviously wrong, which should not be a straight line.
 
  • #7
Pythagorean
Gold Member
4,191
255
no idea, but here's how I've successfully called the ode integrator in one of my python programs:

Code:
result,info = odeint(func,y0, arange(0,tf,dt), args = (I,D),full_output=1)
compared to your:

Code:
track= odeint(func, y0, t, args=(k, c, Zr))
the only difference I see is full_output argument and I called the second output argument.
 
  • #8
chiro
Science Advisor
4,790
131
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()
Even assembler code is easier to read than this!
 
  • #9
Even assembler code is easier to read than this!
I think this should be easier to read.

# -*- 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 & 的d(px)/dt=y[0]; d(py)/dt=y[1]; d(pz)/dt=y[2]; d(x)/dt=y[3]; d(y)/dt=y[4]; d(z)/dt=y[5]
# apply func formula

return (

d(px)/dt,

d(py)/dt,

d(pz)/dt,

d(x)/dt,
d(y)/dt,
d(z)/dt)

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()
 
  • #10
no idea, but here's how I've successfully called the ode integrator in one of my python programs:

Code:
result,info = odeint(func,y0, arange(0,tf,dt), args = (I,D),full_output=1)
compared to your:

Code:
track= odeint(func, y0, t, args=(k, c, Zr))
the only difference I see is full_output argument and I called the second output argument.
thanks for reply.

what does "info" means? and what is usefulness of "full_output=1"?
 

Related Threads on Python ode problem (need help)

Replies
17
Views
1K
  • Last Post
Replies
3
Views
2K
  • Last Post
Replies
3
Views
2K
  • Last Post
Replies
1
Views
3K
  • Last Post
Replies
8
Views
6K
  • Last Post
Replies
11
Views
703
Replies
3
Views
531
  • Last Post
Replies
11
Views
830
  • Last Post
Replies
10
Views
9K
Replies
6
Views
813
Top