Solving a Set of Equations with Python Odeint | Code and Example

  • Thread starter goldin2008
  • Start date
  • Tags
    Ode Python
In summary, Odeint is a function in the SciPy library that is used for numerical integration of ordinary differential equations. It uses an adaptive algorithm to solve the equations and provides accurate results. To use Odeint in your code, you first need to import it from the SciPy library. The inputs for Odeint are the differential equations, initial conditions, and the time points at which the solutions are to be evaluated. To solve a set of equations with Odeint, you first need to define the equations, initial conditions, and time points as mentioned above. Odeint can be used for stiff equations and automatically detects stiffness in the equations, but for highly stiff equations, specialized methods should be used.
  • #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()
 
Technology news on Phys.org
  • #2
goldin2008 said:
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
Mark44 said:
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 mathdef 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
And the error message?
 
  • #6
Pythagorean said:
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
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
goldin2008 said:
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
chiro said:
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
Pythagorean said:
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"?
 

What is Odeint in Python?

Odeint is a function in the SciPy library that is used for numerical integration of ordinary differential equations. It uses an adaptive algorithm to solve the equations and provides accurate results.

How do I import Odeint in Python?

To use Odeint in your code, you first need to import it from the SciPy library. This can be done by using the following code:
from scipy.integrate import odeint

What are the inputs for Odeint?

The inputs for Odeint are the differential equations, initial conditions, and the time points at which the solutions are to be evaluated. The equations can be defined as a function, and the initial conditions can be given as a list or array. The time points are usually defined using the np.linspace function.

How do I solve a set of equations with Odeint?

To solve a set of equations with Odeint, you first need to define the equations, initial conditions, and time points as mentioned above. Then, you can use the odeint function and pass in the necessary inputs. The function will return an array of solutions at the specified time points.

Can Odeint be used for stiff equations?

Yes, Odeint can be used for stiff equations. It automatically detects stiffness in the equations and adapts the integration algorithm accordingly. However, for highly stiff equations, it is recommended to use specialized methods like ode or solve_ivp.

Similar threads

  • Programming and Computer Science
Replies
2
Views
1K
  • Programming and Computer Science
Replies
1
Views
1K
  • Programming and Computer Science
Replies
1
Views
2K
  • Programming and Computer Science
Replies
16
Views
1K
  • Programming and Computer Science
Replies
5
Views
2K
  • Programming and Computer Science
Replies
15
Views
1K
  • Programming and Computer Science
Replies
8
Views
758
  • Programming and Computer Science
Replies
1
Views
1K
Replies
1
Views
1K
  • Programming and Computer Science
Replies
4
Views
499
Back
Top