Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python ode problem (need help)

  1. Nov 4, 2011 #1
    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()
     
  2. jcsd
  3. Nov 5, 2011 #2

    Mark44

    Staff: Mentor

    You're kidding, right? You don't really expect that we're going to look at that huge wad of code, do you?
     
  4. Nov 5, 2011 #3
    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()
     
  5. Nov 5, 2011 #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: Nov 5, 2011
  6. Nov 5, 2011 #5

    Pythagorean

    User Avatar
    Gold Member

    And the error message?
     
  7. Nov 6, 2011 #6
    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.
     
  8. Nov 7, 2011 #7

    Pythagorean

    User Avatar
    Gold Member

    no idea, but here's how I've successfully called the ode integrator in one of my python programs:

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

    Code (Text):
    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.
     
  9. Nov 7, 2011 #8

    chiro

    User Avatar
    Science Advisor

    Even assembler code is easier to read than this!
     
  10. Nov 7, 2011 #9
    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()
     
  11. Nov 7, 2011 #10
    thanks for reply.

    what does "info" means? and what is usefulness of "full_output=1"?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Python ode problem (need help)
  1. Help in Python rk4 (Replies: 1)

  2. Python problem (Replies: 3)

  3. Basic python help? (Replies: 3)

Loading...