Python integrate.odeint方法代码示例

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def integrate(self,s0,t):
        """Integrates motion in the CRTBP given initial conditions
        
        This method returns a star's position vector in the rotating frame of 
        the Circular Restricted Three Body Problem.  
        
        Args:
            s0 (integer 1x6 array):
                Initial state vector consisting of stacked position and velocity vectors
                in normalized units
            t (integer):
                Times in normalized units

        Returns:
            s (integer nx6 array):
                State vector consisting of stacked position and velocity vectors
                in normalized units
        """
        
        EoM = lambda s,t: self.equationsOfMotion_CRTBP(t,s)
             
        s = itg.odeint(EoM, s0, t, full_output = 0,rtol=2.5e-14,atol=1e-22)
        
        return s 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def integrate_nonlinear_piecewise(self, X_l, U, sigma):
        """
        Piecewise integration to verfify accuracy of linearization.
        :param X_l: Linear state evolution
        :param U: Linear input evolution
        :param sigma: Total time
        :return: The piecewise integrated dynamics
        """
        X_nl = np.zeros_like(X_l)
        X_nl[:, 0] = X_l[:, 0]

        for k in range(self.K - 1):
            X_nl[:, k + 1] = odeint(self._dx, X_l[:, k],
                                    (0, self.dt * sigma),
                                    args=(U[:, k], U[:, k + 1], sigma))[1, :]

        return X_nl 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def coffee_cup_dependent(kappa_hat, T_env, alpha):
    # Initial temperature and time
    time = np.linspace(0, 200, 150)            # Minutes
    T_0 = 95                                   # Celsius

    # The equation describing the model
    def f(T, time, alpha, kappa_hat, T_env):
        return -alpha*kappa_hat*(T - T_env)

    # Solving the equation by integration.
    temperature = odeint(f, T_0, time, args=(alpha, kappa_hat, T_env))[:, 0]

    # Return time and model results
    return time, temperature


# Create a model from the coffee_cup_dependent function and add labels 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def run(self, kappa_hat, T_env, alpha):
        # Initial temperature and time
        time = np.linspace(0, 200, 150)            # Minutes
        T_0 = 95                                   # Celsius

        # The equation describing the model
        def f(T, time, alpha, kappa_hat, T_env):
            return -alpha*kappa_hat*(T - T_env)

        # Solving the equation by integration.
        values = odeint(f, T_0, time, args=(alpha, kappa_hat, T_env))[:, 0]

        # Return time and model results
        return time, values


# Initialize the model 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def compute_depth_gravity_profiles(pressures, densities, surface_gravity, outer_radius):
    gravity = [surface_gravity] * len(pressures) # starting guess
    n_gravity_iterations = 5
    for i in range(n_gravity_iterations):    
        # Integrate the hydrostatic equation
        # Make a spline fit of densities as a function of pressures
        rhofunc = UnivariateSpline(pressures, densities)
        # Make a spline fit of gravity as a function of depth
        gfunc = UnivariateSpline(pressures, gravity)
        
        # integrate the hydrostatic equation
        depths = np.ravel(odeint((lambda p, x: 1./(gfunc(x) * rhofunc(x))), 0.0, pressures))
        
        radii = outer_radius - depths
        
        rhofunc = UnivariateSpline(radii[::-1], densities[::-1])
        poisson = lambda p, x: 4.0 * np.pi * burnman.constants.G * rhofunc(x) * x * x
        gravity = np.ravel(odeint(poisson, surface_gravity*radii[0]*radii[0], radii))
        gravity = gravity / radii / radii
    return depths, gravity


# BEGIN USER INPUTS

# Declare the rock we want to use 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def calculate_discretization(self, X, U, sigma):
        """
        Calculate discretization for given states, inputs and total time.

        :param X: Matrix of states for all time points
        :param U: Matrix of inputs for all time points
        :param sigma: Total time
        :return: The discretization matrices
        """
        for k in range(self.K - 1):
            self.V0[self.x_ind] = X[:, k]
            V = np.array(odeint(self._ode_dVdt, self.V0, (0, self.dt), args=(U[:, k], U[:, k + 1], sigma))[1, :])

            # using \Phi_A(\tau_{k+1},\xi) = \Phi_A(\tau_{k+1},\tau_k)\Phi_A(\xi,\tau_k)^{-1}
            # flatten matrices in column-major (Fortran) order for CVXPY
            Phi = V[self.A_bar_ind].reshape((self.n_x, self.n_x))
            self.A_bar[:, k] = Phi.flatten(order='F')
            self.B_bar[:, k] = np.matmul(Phi, V[self.B_bar_ind].reshape((self.n_x, self.n_u))).flatten(order='F')
            self.C_bar[:, k] = np.matmul(Phi, V[self.C_bar_ind].reshape((self.n_x, self.n_u))).flatten(order='F')
            self.S_bar[:, k] = np.matmul(Phi, V[self.S_bar_ind])
            self.z_bar[:, k] = np.matmul(Phi, V[self.z_bar_ind])

        return self.A_bar, self.B_bar, self.C_bar, self.S_bar, self.z_bar 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def integrate_nonlinear_full(self, x0, U, sigma):
        """
        Simulate nonlinear behavior given an initial state and an input over time.
        :param x0: Initial state
        :param U: Linear input evolution
        :param sigma: Total time
        :return: The full integrated dynamics
        """
        X_nl = np.zeros([x0.size, self.K])
        X_nl[:, 0] = x0

        for k in range(self.K - 1):
            X_nl[:, k + 1] = odeint(self._dx, X_nl[:, k],
                                    (0, self.dt * sigma),
                                    args=(U[:, k], U[:, k + 1], sigma))[1, :]

        return X_nl 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def calculate_discretization(self, X, U):
        """
        Calculate discretization for given states, inputs and total time.

        :param X: Matrix of states for all time points
        :param U: Matrix of inputs for all time points
        :return: The discretization matrices
        """
        for k in range(self.K - 1):
            self.V0[self.x_ind] = X[:, k]
            V = np.array(odeint(self._ode_dVdt, self.V0, (0, self.dt), args=(U[:, k], U[:, k + 1]))[1, :])

            # flatten matrices in column-major (Fortran) order for CVXPY
            Phi = V[self.A_bar_ind].reshape((self.n_x, self.n_x))
            self.A_bar[:, k] = Phi.flatten(order='F')
            self.B_bar[:, k] = np.matmul(Phi, V[self.B_bar_ind].reshape((self.n_x, self.n_u))).flatten(order='F')
            self.C_bar[:, k] = np.matmul(Phi, V[self.C_bar_ind].reshape((self.n_x, self.n_u))).flatten(order='F')
            self.z_bar[:, k] = np.matmul(Phi, V[self.z_bar_ind])

        return self.A_bar, self.B_bar, self.C_bar, self.z_bar 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def _do_problem(self, problem, integrator, method='adams'):

        # ode has callback arguments in different order than odeint
        f = lambda t, z: problem.f(z, t)
        jac = None
        if hasattr(problem, 'jac'):
            jac = lambda t, z: problem.jac(z, t)

        ig = ode(f, jac)
        ig.set_integrator(integrator,
                          atol=problem.atol/10,
                          rtol=problem.rtol/10,
                          method=method)
        ig.set_initial_value(problem.z0, t=0.0)
        z = ig.integrate(problem.stop_t)

        assert_(ig.successful(), (problem, method))
        assert_(problem.verify(array([z]), problem.stop_t), (problem, method)) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def likelihood(parameter_vector):

    parameter_vector = 10**np.array(parameter_vector)

    #Solve ODE system given parameter vector
    yout = odeint(odefunc, y0, tspan, args=(parameter_vector,))

    cout = yout[:, 2]

    #Calculate log probability contribution given simulated experimental values.
    
    logp_ctotal = np.sum(like_ctot.logpdf(cout))
    
    #If simulation failed due to integrator errors, return a log probability of -inf.
    if np.isnan(logp_ctotal):
        logp_ctotal = -np.inf
      
    return logp_ctotal


# Add vector of rate parameters to be sampled as unobserved random variables in DREAM with uniform priors. 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def _do_problem(self, problem):
        t = arange(0.0, problem.stop_t, 0.05)

        # Basic case
        z, infodict = odeint(problem.f, problem.z0, t, full_output=True)
        assert_(problem.verify(z, t))

        # Use tfirst=True
        z, infodict = odeint(lambda t, y: problem.f(y, t), problem.z0, t,
                             full_output=True, tfirst=True)
        assert_(problem.verify(z, t))

        if hasattr(problem, 'jac'):
            # Use Dfun
            z, infodict = odeint(problem.f, problem.z0, t, Dfun=problem.jac,
                                 full_output=True)
            assert_(problem.verify(z, t))

            # Use Dfun and tfirst=True
            z, infodict = odeint(lambda t, y: problem.f(y, t), problem.z0, t,
                                 Dfun=lambda t, y: problem.jac(y, t),
                                 full_output=True, tfirst=True)
            assert_(problem.verify(z, t)) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def _solve(self, y0, a_normalize):
        # solve with critical point at a=1.0, lna=0.
        y = odeint(self.ode, y0, self.lna, tcrit=[0.], atol=0)

        v1 = []
        v2 = []
        for yi, lnai in zip(y, self.lna):
            D1, F1, D2, F2 = yi
            D1p, F1p, D2p, F2p = self.ode(yi, lnai)
            v1.append((D1, F1, F1p))
            v2.append((D2, F2, F2p))

        v1 = np.array(v1)
        v2 = np.array(v2)

        ind = abs(self.lna - np.log(a_normalize)).argmin()
        # normalization to 1 at a=a_normalize
        v1 /= v1[ind][0]
        v2 /= v2[ind][0]
        return v1, v2 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def integrate_nonlinear_full(self, x0, U, sigma):
        """
        Simulate nonlinear behavior given an initial state and an input over time.
        :param x0: Initial state
        :param U: Linear input evolution
        :param sigma: Total time
        :return: The full integrated dynamics
        """
        X_nl = np.zeros([x0.size, self.K])
        X_nl[:, 0] = x0

        for k in range(self.K - 1):
            X_nl[:, k + 1] = odeint(self._dx, X_nl[:, k], (0, self.dt * sigma), args=(U[:, k], U[:, k + 1], sigma))[1,
                             :]

        return X_nl 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def derivs(M, t, Meq, w, w1, T1, T2):
    '''Bloch equations in rotating frame.

    Args:
        w:              Larmor frequency :math:`2\\pi\\gamma B_0` [kRad / s].
	    w1 (complex):   B1 rotation frequency :math:`2\\pi\\gamma B_1`  [kRad / s].
        T1:             longitudinal relaxation time.
        T2:             transverse relaxation time.
        M:              magnetization vector.
        Meq:            equilibrium magnetization.
        t:              time vector (needed for scipy.integrate.odeint).

    Returns:
        integrand :math:`\\frac{dM}{dt}`
    
    '''
    
    dMdt = np.zeros_like(M)
    dMdt[0] = -M[0]/T2+M[1]*w+M[2]*w1.real
    dMdt[1] = -M[0]*w-M[1]/T2+M[2]*w1.imag
    dMdt[2] = -M[0]*w1.real-M[1]*w1.imag+(Meq-M[2])/T1
    return dMdt 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import odeint [as 别名]
def follow(self, Rstart, Zstart, angles, rtol=None, backward=False):
        """Follow magnetic field lines from (Rstart, Zstart) locations
        to given toroidal angles. 
        If backward is True then the field lines are followed in the
        -B direction"""

        Rstart = np.array(Rstart)
        Zstart = np.array(Zstart)
        
        array_shape = Rstart.shape
        assert Zstart.shape == array_shape

        evolving = np.ones(array_shape)
        
        # (R,Z,length) with length=0 initially
        position = np.column_stack((Rstart, Zstart, np.zeros(array_shape))).flatten()
        result = odeint(self.fieldDirection, position, angles,
                        args=(evolving, backward), rtol=rtol)

        return result.reshape(angles.shape + array_shape + (3,))