Python integrate.odeint函数代码示例

def crystal_embed(e,eta):
    """
    Evaluates crystal embedding potential by solving the one-dimensional
    Schroedinger equation through one unit cell in both directions.
    """
    y0=[1.0,0.0,0.0,0.0]
    z1=np.linspace(-z_left,-z_left+alat,101)
    sol1=odeint(schroed,y0,z1,args=(e,eta))
    z2=np.linspace(-z_left+alat,-z_left,101)
    sol2=odeint(schroed,y0,z2,args=(e,eta))
    psi1=complex(sol1[50,0],sol1[50,1])
    psi1_prime=complex(sol1[50,2],sol1[50,3])
    psi2=complex(sol2[50,0],sol2[50,1])
    psi2_prime=complex(sol2[50,2],sol2[50,3])
    wronskian=psi1*psi2_prime-psi2*psi1_prime
    psi1=complex(sol1[100,0],sol1[100,1])
    psi2=complex(sol2[100,0],sol2[100,1])
    cos_ka=0.5*(psi1+psi2)
    ka=arccos(cos_ka)
    if ka.imag<0: ka=np.conj(ka)
    exp_ka=cos_ka+1.0j*sqrt(1.0-cos_ka**2)
    emb_cryst=0.5*wronskian/(exp_ka-psi2)
    if emb_cryst.imag>0.0:
        exp_ka=cos_ka-1.0j*sqrt(1.0-cos_ka**2)
        emb_cryst=0.5*wronskian/(exp_ka-psi2)
    return emb_cryst

def getEta(data,t,xi):
	global samp, ETA, time, agent, XI
	if t < len(data):
		return data[t]
	else:
		XI   = xi # update input function from paramteter
		if len(data) == 0:
			ETA0 = getInitialEta(agent.beta,agent.gamma,XI)
			data.append(ETA0[:])

		for T in range(len(data),t+1):
			# TODO: should this be samp*t so that accuracy is not lost far from 0???
			logging.info('solving ode @ t='+str(T)+', using '+str(samp)+' sub-samples')
			time = linspace(0,T,samp) #(start,end,nSamples)
			etadot_0 = [0,0,0,0,0]	#assumption of 1st order model
			#get arrays of data len=samp*t
			ETA[0] = integrate.odeint(eta1Func,[data[0][0],etadot_0[0]],time)
			ETA[1] = integrate.odeint(eta2Func,[data[0][1],etadot_0[1]],time)
			ETA[2] = integrate.odeint(eta3Func,[data[0][2],etadot_0[2]],time)
			ETA[3] = integrate.odeint(eta4Func,[data[0][3],etadot_0[3]],time)
			ETA[4] = integrate.odeint(eta5Func,[data[0][4],etadot_0[4]],time)
			logging.debug('len(result)='+str(len(ETA[0][:,0])))
			# restructure ETA using [eta#][time , eta_or_dEta] )
			E = [ETA[0][-1,0],\
			     ETA[1][-1,0],\
			     ETA[2][-1,0],\
			     ETA[3][-1,0],\
			     ETA[4][-1,0]]
			data.append(E)
		return data[t]

def ode_wrapper(fun,x0,tsp):
    # array of times at which dynamic uncertainty happens
    tdisc = np.arange(tsp[0],tsp[-1]+DT,DT)

    # determine (ad hoc) whether the called function takes 2 or 3 args
    flag_stoch = False
    try:
        dy = fun(x0,tsp[0],0.0)
        flag_stoch = True
    except TypeError:
        dy = fun(x0,tsp[0])
        flag_stoch = False

    yt = np.zeros((len(tdisc),len(x0)))
    y0 = x0.copy()
    for k in range(len(tdisc)-1):
        if flag_stoch:
            # compute the noise
            v = np.random.normal(scale=q_w)
            yp = sp.odeint(fun,y0,tdisc[k:k+2],args=(v,))
        else:
            yp = sp.odeint(fun,y0,tdisc[k:k+2])
        yt[k,:] = y0.copy()
        y0 = yp[-1,:].copy()
    yt[-1,:] = y0.copy()
    # interpolate in the new history to match the passed-in history tsp
    ysp = np.zeros((len(tsp),len(x0)))
    for k in range(len(x0)):
        ysp[:,k] = np.interp(tsp,tdisc,yt[:,k])
    return ysp

def compute_jacobian(star, differences, surface_guesses, core_guesses, core_masses, surface_masses, mass_step):
    jacobian =(np.zeros((4,4)))

    for i in range(0,4):
        guess_star = copy(star)
        if i == 0:
            step_size = core_guesses[0]*0.01
            guess_star.core_pressure = core_guesses[0] + step_size
        elif i == 1:
            step_size = core_guesses[1]*0.01
            guess_star.core_temp     = core_guesses[1] + step_size
        elif i == 2:
            step_size = surface_guesses[2]*0.01
            guess_star.total_radius  = surface_guesses[2] + step_size
        elif i == 3:
            step_size = surface_guesses[3]*0.01
            guess_star.total_lum     = surface_guesses[3] + step_size

        new_surface = inward_start(guess_star)
        new_core = outward_start(guess_star, mass_step)

        new_differences = difference_is(odeint(derivatives, new_core, core_masses, args=(guess_star,)),
                            odeint(derivatives, new_surface, surface_masses, args=(guess_star,)))
        jacobian[:,i] = np.asarray((new_differences - differences)/step_size)

    return np.linalg.inv(jacobian)

    def __init__(self, f, u0, s, t, dfdu=None):
        self.f = f
        self.t = np.array(t, float).copy()
        self.s = np.array(s, float).copy()

        if self.s.ndim == 0:
            self.s = self.s[np.newaxis]

        if dfdu is None:
            dfdu = ddu(f)
        self.dfdu = dfdu

        u0 = np.array(u0, float)
        if u0.ndim == 1:
            # run up to t[0]
            f = lambda u, t : self.f(u, s)
            assert t[0] >= 0 and t.size > 1
            N0 = int(t[0] / (t[-1] - t[0]) * t.size)
            u0 = odeint(f, u0, np.linspace(0, t[0], N0+1))[-1]

            # compute a trajectory
            self.u = odeint(f, u0, t - t[0])
        else:
            assert (u0.shape[0],) == t.shape
            self.u = u0.copy()

        self.dt = t[1:] - t[:-1]
        self.uMid = 0.5 * (self.u[1:] + self.u[:-1])
        self.dudt = (self.u[1:] - self.u[:-1]) / self.dt[:,np.newaxis]

    def Advect(self, flow_type = None, t0 = 0, dt = 0, extra_args = None):

        if flow_type is not None:
    
            # Number of particles
            num_particles = len(self.Particles);
        
            # Get the positions of all of the particles
            x0, y0, z0 = self.GetCoordinates();
        
            # Time vector
            tf = t0 + dt;
        
            # Time vector
            t = np.array([t0, tf]);
        
            # Initial positions as a numpy vector
            # in the form compatible with the 
            # velocity functions
            xy0 = np.array(x0 + y0);
            # xy0 = np.array([x0, y0]);
    
            # Choose between fields
            if "hama" in flow_type.lower():
                xy = (odeint(HamaVelocity, y0 = xy0, t = t, args = (extra_args,)))[-1, :];
                
            elif "cpipe" in flow_type.lower():
                xy = (odeint(cpipe, y0 = xy0, t = t, args = (extra_args,)))[-1, :];
                
            # Extract the new positions
            x_new, y_new = Parse_Vector_2d(xy);
            
            # Set the new coordinates
            self.SetCoordinates(x = x_new, y = y_new);

def exercise3part1(T=200, Deltat=0.1, epsilon=0.5):
    
    "In this exercise we want to perturb the parameters of the Lorenz system to see how the shape of the Lorenz attractor change"
    
    y0=np.random.randn(3)
    t = np.arange(0.,T, Deltat)
    data = integrate.odeint(vectorfield, y0, t=t,args=(10,28,8/3))
    data_sigma=integrate.odeint(vectorfield, y0, t=t,args=(10+epsilon,28,8/3))
    data_rho=integrate.odeint(vectorfield, y0, t=t,args=(10,28+epsilon,8/3))
    data_beta=integrate.odeint(vectorfield, y0, t=t,args=(10,28,8/3+epsilon))
    fig = plt.figure()
    ax0 = fig.add_subplot(2, 2, 1, projection='3d')
    ax1 = fig.add_subplot(2, 2, 2, projection='3d')
    ax2 = fig.add_subplot(2, 2, 3, projection='3d')
    ax3 = fig.add_subplot(2, 2, 4, projection='3d')
    
    ax0.plot(data[:,0],data[:,1],data[:,2])
    ax1.plot(data_sigma[:,0],data_sigma[:,1],data_sigma[:,2])
    ax2.plot(data_rho[:,0],data_rho[:,1],data_rho[:,2])
    ax3.plot(data_beta[:,0],data_beta[:,1],data_beta[:,2])

    ax0.set_title('Standard parameters')
    ax1.set_title('Sigma perturbed')
    ax2.set_title('Rho perturbed')
    ax3.set_title('Beta perturbed')

    pylab.show()
    print 'From the plot, the system seems to be more sensibile in beta perturbations'

def test_repeated_t_values():
    """Regression test for gh-8217."""

    def func(x, t):
        return -0.25*x

    t = np.zeros(10)
    sol = odeint(func, [1.], t)
    assert_array_equal(sol, np.ones((len(t), 1)))

    tau = 4*np.log(2)
    t = [0]*9 + [tau, 2*tau, 2*tau, 3*tau]
    sol = odeint(func, [1, 2], t, rtol=1e-12, atol=1e-12)
    expected_sol = np.array([[1.0, 2.0]]*9 +
                            [[0.5, 1.0],
                             [0.25, 0.5],
                             [0.25, 0.5],
                             [0.125, 0.25]])
    assert_allclose(sol, expected_sol)

    # Edge case: empty t sequence.
    sol = odeint(func, [1.], [])
    assert_array_equal(sol, np.array([], dtype=np.float64).reshape((0, 1)))

    # t values are not monotonic.
    assert_raises(ValueError, odeint, func, [1.], [0, 1, 0.5, 0])
    assert_raises(ValueError, odeint, func, [1, 2, 3], [0, -1, -2, 3])

def computeTrajectories(func, E0=np.zeros(3), **keywords):
  """Movement of electron and ion under a constant magnetic field.
  
  Positional arguments:
  func -- the name of the function computing dy/dt at time t0
  Keyword arguments:
  E0 -- Constant component of the electric field
  All other keyword arguments are collected in a 'keywords' dictionary 
  and are specific to each func."""
    
  from scipy.integrate import odeint
    
  # Processes the func specific arguments
  re0, rp0 = "ri" in keywords.keys() and keywords["ri"] or [np.zeros(3),np.zeros(3)]
  if "vi" in keywords.keys():   # Initial velocity
    v0 = keywords["vi"]
  else:
    v0 = np.array([0,1,0])
  wce, wcp = "wc" in keywords.keys() and keywords["wc"] or [0,0]

  tf = 350; NPts = 10*tf; t = np.linspace(0,tf,NPts)  # Time points

  # Integration of the equations of movement
  Q0 = np.concatenate((re0,v0))                  # Initial values
  if "wc" in keywords.keys():
    keywords["wc"] = wce
  Qe = odeint(func, Q0, t, args=(-q/me,E0,B0,keywords))  # Trajectory for the "electron"
  Q0 = np.concatenate((rp0,v0))                  # Initial values
  if "wc" in keywords.keys():
    keywords["wc"] = wcp
  Qp = odeint(func, Q0, t, args=(q/Mp,E0,B0,keywords))   # Trajectory for the "ion"
  
  return Qe, Qp

def Rhalf(Rp,Rn,a,b):
    Rn /= Rp
    c = (a*a + b*b)**0.5 / Rp
    r = linspace(0.001,1,100)

    def numer(f,r):
        lo,hi = theta_lohi(Rn,c,0,r)
        return r*r * (sin(lo)-sin(hi))

    def denom(f,r):
        lo,hi = theta_lohi(Rn,c,0,r)
        return r * (hi - lo)

    def arc(f,r):
        lo,hi = theta_lohi(Rn,c,d,r)
        return r * (hi - lo)

    up = odeint(numer,[0],r)
    dn = odeint(denom,[0],r)

    d = up[-1,0]/dn[-1,0]
    area = odeint(arc,[0],r)

    for i in range(1,len(r)):
        if area[i-1] < 0.5*area[-1] and area[i] > 0.5*area[-1]:
            return Rp * r[i]

def integrate(state, N_t=20, T=1., pts=None):
    """
    flow forward integration

    Points pts are carried along the flow without affecting it
    """
    assert(N)
    assert(DIM)
    assert(SIGMA)
    gaussian.N = N
    gaussian.DIM = DIM
    gaussian.SIGMA = SIGMA

    t_span = np.linspace(0. ,T , N_t )
    #print 'forward integration: SIGMA = ' + str(SIGMA) + ', N = ' + str(N) + ', DIM = ' + str(DIM) + ', N_t = ' + str(N_t) + ', T = ' + str(T)

    if parallel:
        odef = ode_function
    else:
        odef = ode_function_single
    
    if pts == None:
        y_span = odeint( odef , state , t_span)
        return (t_span, y_span)
    else:
        y_span = odeint( odef , np.hstack((state,pts.flatten())) , t_span)
        return (t_span, y_span[:,0:get_dim_state()], y_span[:,get_dim_state():])

def poly(n,zi,dz,filename):
    '''
    Integrates the L-E equation and returns and writes to a file the important values.
    n is the index, zi is the interval to use the polynomial solution, dz is the step size, filename is the prefix for the output file (filename+.fits)
    '''
    #set up equation and variables
    LEeq = LEeqn(n)
    zs = np.arange(0.,zi,dz)
    nz = np.size(zs)

    #integrate over initial interval using polynomial
    y = intg.odeint(thPN,np.array([1.,0.]),zs)

    #integrate until th<0
    while y[-1,0]>0.:
        zs2 = np.arange(zs[-1],zs[-1]+200.*dz,dz)
        zs = np.append(zs,zs2)
        y = np.append(y,intg.odeint(LEeq,y[-1],zs2),axis=0)
        #write progress
        sys.stdout.write("\rz={0},th={1}".format(zs[-1],y[-1,0]))
        sys.stdout.flush()
    sys.stdout.write("\n")

    #only keep values where th>0
    ngood= np.size(np.where(y[:,0]>0.))
    zs = zs[:ngood]
    y = y[:ngood]

    #prepare output array
    data = np.array([zs,y[:,0],-(zs**2)*y[:,1],(-3./zs)*y[:,1]])

    #write and return data
    fits.writeto(filename+'.fits',data)
    return data

	def delta_r(self):
		"""
		returns the coeficient delta indicating the RADIAL stability of the movement
		"""
		time = np.linspace(0, 5, 1000)
		def sm1(x, t):
			if x[0] < 0:
				return np.array([x[1], 0]) 
			else:
				return np.array([x[1], self.eta*self.Ez0(x[0])])
		traj = odeint(sm1,[0,1.],time)
		if traj[-1,0] < 0:
			ind = np.where(traj[:,0] > 0)[0][-1]
			tm = [time[ind+1],time[ind]]
			zm = [traj[ind+1,0],traj[ind,0]]
			period = np.interp(0,zm,tm)
			def sm2(x, t):
				t = t - np.floor(t/period)*period
				pos = np.interp(t, time, traj[:,0])
				psi = self.alpha + self.b**2/4. + (self.eta/2.)*self.d2V(pos)
				return np.array([x[1], psi*x[0], x[3], psi*x[2]])
			x0 = [1, 0, 0, 1]
			sol = odeint(sm2, x0, np.linspace(0, period, 1000))
			delta = np.abs((sol[-1,0] + sol[-1,3])/2)
			beta  = np.arccos((sol[-1,0] + sol[-1,3])/2)/np.pi
			if self.b == 0.:
				return delta
			else:
				out = np.log(delta + np.sqrt(np.abs(delta**2-1)))/period - self.b
				return out
		else:
			return 1.0

def runClock(dut, t_start, init_cond):
        global ts, dt, sig_jit
        global p0, z0
        global t_ab, t_ba

        alpha = random.normal(0, sig_jit)
        beta  = random.normal(0, sig_jit)
        print (alpha, beta) 
        t1 = arange(0, t_ab + alpha, dt)
        t2 = arange(0, t_ba + beta,  dt)

        concatenate((t1, [t_ab + alpha])) # 0 ~ t_ab + jitter
        concatenate((t2, [t_ba + beta ])) # 0 ~ t_ba + jitter
        t2 = map(lambda x: x + t1[-1], t2) # 0 ~ t_ab+jitter, t_ab+jitter ~ (t_ab+t_ba+2*jitter)
        
        state_half = odeint(dut, init_cond,      t1, args=(t_start, p0, z0))
        eventAtA(state_half[-1,:]) # comparator / DAC
        
        if(t_ba > dt): # when t_ba = 0 sec (NRZ DAC)
                state_full = odeint(dut, state_half[-1], t2, args=(t_start, p0, z0))
                eventAtB(state_full[-1,:]) # DAC (RZ)
                t     = concatenate((t1,         t2        ))
                state = concatenate((state_half, state_full))
        else:
                t     = t1 
                state = state_half 

        return (t, state)

def Gate(x,t):
  global pars
  xinf = ActivationCurve(V[:,1], pars[0], pars[1])
  tau = Kinetic(V[:,1], pars[2], pars[3], pars[4], pars[5])
  return (xinf-x)/tau

  # the following function is the optimisation part of the algorithm. The two chanels are modeled with initial conditions 
  # and the result is compared to the data, taking account the parameter physiological ranges for identificationdef FullTrace(p,time,y):
  pars = p[0:6]
  rinit1 = ActivationCurve(V[:,0], pars[0], pars[1]) # initial values
  r1 = odeint(Gate,rinit1,time)
  pars = p[7:13]
  rinit2 = ActivationCurve(V[:,0], pars[0], pars[1]) # initial values
  r2 = odeint(Gate,rinit2,time)
  I = p[6]*r1*(V[:,1]-E)+p[13]*r2*(V[:,1]-E)
  
  
  # this plot is to see in real time the reconstructed curves trying to fit the data (very fun)
  plot(time,y,'b'); hold(True)
  plot(time,I,'--r');hold(False)
  draw()  
  
  # this part is added to constrain the algorithm with physiological search ranges
  A = 0
  for a in arange(len(p)):
    if p[a] > HB[a] :
      A = A + (p[a] - HB[a])*1e8
    if p[a] < LB[a] :
      A = A + (LB[a]-p[a])*1e8
  return sum(array(y - I)*array(y - I))+A