Python scipy.integrate.odeint() Examples

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 

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 

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. 

def simulate(self, steps, time, solution = solve_type.RK4):
        """!
        @brief Performs static simulation of oscillatory network based on Hodgkin-Huxley neuron model.
        @details Output dynamic is sensible to amount of steps of simulation and solver of differential equation.
                  Python implementation uses 'odeint' from 'scipy', CCORE uses classical RK4 and RFK45 methods,
                  therefore in case of CCORE HHN (Hodgkin-Huxley network) amount of steps should be greater than in
                  case of Python HHN.

        @param[in] steps (uint): Number steps of simulations during simulation.
        @param[in] time (double): Time of simulation.
        @param[in] solution (solve_type): Type of solver for differential equations.
        
        @return (tuple) Dynamic of oscillatory network represented by (time, peripheral neurons dynamic, central elements
                dynamic), where types are (list, list, list).
        
        """
        
        return self.simulate_static(steps, time, solution) 

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 

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 

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 

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 

def step(self, action):
        if action == 0:
            u = -self._max_u
        elif action == 1:
            u = 0.
        else:
            u = self._max_u

        self._last_u = u

        u += np.random.uniform(-self._noise_u, self._noise_u)
        new_state = odeint(self._dynamics, self._state, [0, self._dt],
                           (u,))

        self._state = np.array(new_state[-1])
        self._state[0] = normalize_angle(self._state[0])

        if np.abs(self._state[0]) > np.pi * .5:
            reward = -1.
            absorbing = True
        else:
            reward = 0.
            absorbing = False

        return self._state, reward, absorbing, {} 

def step(self, action):
        u = self._bound(action[0], -self._max_u, self._max_u)
        new_state = odeint(self._dynamics, self._state, [0, self._dt],
                           (u,))

        self._state = np.array(new_state[-1])
        self._state[0] = normalize_angle(self._state[0])

        if abs(self._state[0]) > np.pi / 2:
            absorbing = True
            reward = -10000
        else:
            absorbing = False
            Q = np.diag([3.0, 0.1, 0.1])

            x = self._state

            J = x.dot(Q).dot(x)

            reward = -J

        return self._state, reward, absorbing, {} 

def step(self, action):
        action = self._discrete_actions[action[0]]
        sa = np.append(self._state, action)
        new_state = odeint(self._dpds, sa, [0, self._dt])

        self._state = new_state[-1, :-1]

        if self._state[0] < -self.max_pos or \
                np.abs(self._state[1]) > self.max_velocity:
            reward = -1
            absorbing = True
        elif self._state[0] > self.max_pos and \
                np.abs(self._state[1]) <= self.max_velocity:
            reward = 1
            absorbing = True
        else:
            reward = 0
            absorbing = False

        return self._state, reward, absorbing, {} 

def div_dep_ext_dt(div, t, d12, d21, mu1, mu2, k_d1, k_d2, covar_mu1, covar_mu2):
	div1 = div[0]
	div2 = div[1]
	div3 = div[2]
	div13 = div[0] + div[2]
	div23 = div[1] + div[2]
	lim_d2 = max(0, 1 - (div1 + div3)/k_d1) # Limit dispersal into area 2
	lim_d1 = max(0, 1 - (div2 + div3)/k_d2) # Limit dispersal into area 1
	dS = np.zeros(5)
	dS[3] = d21 * div2 * lim_d1 # Gain area 1
	dS[4] = d12 * div1 * lim_d2 # Gain area 2
	mu1 = mu1 + covar_mu1 * dS[3] / (div13 + 1.)
	mu2 = mu2 + covar_mu2 * dS[4] / (div23 + 1.)
	dS[0] = -mu1 * div1 + mu2 * div3 - dS[4]
	dS[1] = -mu2 * div2 + mu1 * div3 - dS[3]
	dS[2] = -(mu1 + mu2) * div3 + dS[3] + dS[4]
	return dS

#div_int = odeint(div_dep_ext_dt, np.array([1., 1., 0., 0., 0.]), [0, 1], args = (0.2, 0.2, 0.1, 0.1, np.inf, np.inf, 0.2, 0.2))
#div_int 

def __calculate(self, t, step, int_step):
        """!
        @brief Calculates new amplitudes for oscillators in the network in line with current step.
        
        @param[in] t (double): Time of simulation.
        @param[in] step (double): Step of solution at the end of which states of oscillators should be calculated.
        @param[in] int_step (double): Step differentiation that is used for solving differential equation.
        
        @return (list) New states (phases) for oscillators.
        
        """
        
        next_amplitudes = [0.0] * self._num_osc;
        
        for index in range (0, self._num_osc, 1):
            z = numpy.array(self.__amplitude[index], dtype = numpy.complex128, ndmin = 1);
            result = odeint(self.__calculate_amplitude, z.view(numpy.float64), numpy.arange(t - step, t, int_step), (index , ));
            next_amplitudes[index] = (result[len(result) - 1]).view(numpy.complex128);
        
        return next_amplitudes; 

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

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 = complex_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)
        z2 = ig.y

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

def run(self, **parameters):
        self.set_parameters(**parameters)

        self.h0 = self.h_inf(self.V_rest)
        self.m0 = self.m_inf(self.V_rest)
        self.n0 = self.n_inf(self.V_rest)

        initial_conditions = [self.V_rest, self.h0, self.m0, self.n0]

        X = odeint(self.dXdt, initial_conditions, self.time)
        values = X[:, 0]

         # Add info needed by certain spiking features and efel features
        info = {"stimulus_start": self.time[0], "stimulus_end": self.time[-1]}

        return self.time, values, info 

def compute_separatrices(Xss, Js, func, x_range, y_range, t=50, n_sample=500, eps=1e-6):
    ret = []
    for i, x in enumerate(Xss):
        print(x)
        J = Js[i]
        w, v = eig(J)
        I_stable = np.where(np.real(w) < 0)[0]
        print(I_stable)
        for j in I_stable:  # I_unstable
            u = np.real(v[j])
            u = u / np.linalg.norm(u)
            print("u=%f, %f" % (u[0], u[1]))

            # Parameters for building separatrix
            T = np.linspace(0, t, n_sample)
            # all_sep_a, all_sep_b = None, None
            # Build upper right branch of separatrix
            ab_upper = odeint(lambda x, _: -func(x), x + eps * u, T)
            # Build lower left branch of separatrix
            ab_lower = odeint(lambda x, _: -func(x), x - eps * u, T)

            sep = np.vstack((ab_lower[::-1], ab_upper))
            ret.append(sep)
    ret = clip_curves(ret, [x_range, y_range])
    return ret 

def solve_coupled_constecc_solution(F0, e0, phase0, mc, t):
    """
    Compute the solution to the coupled system of equations
    from from Peters (1964) and Barack & Cutler (2004) at
    a given time.

    :param F0: Initial orbital frequency [Hz]
    :param mc: Chirp mass of binary [Solar Mass]
    :param t: Time at which to evaluate solution [s]

    :returns: (F(t), phase(t))
    """

    y0 = np.array([F0, phase0])

    y, infodict = odeint(get_coupled_constecc_eqns, y0, t, args=(mc, e0), full_output=True)

    if infodict["message"] == "Integration successful.":
        ret = y
    else:
        ret = 0

    return ret 

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 

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 

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 

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 

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

def simulate(self, simulation_time, sample_time):
        n_samp = int(np.round(simulation_time/sample_time, 6))+1
        time0 = self.signals['time'][-1]
        time_axis = np.linspace(time0, (n_samp-1)*sample_time+time0, n_samp)
        state0 = np.r_[self.signals['position'][:, -1],
                       self.signals['velocity'][:, -1],
                       self.signals['acceleration'][:, -1]].T
        if time0 != 0.0:
            state0 -= self.state_incr_interp(time0)
        state = odeint(self._ode, state0, time_axis).T
        state += self.state_incr_interp(time_axis)
        self.signals['position'] = np.c_[self.signals['position'],
                                         state[:self.n_dim, 1:n_samp+1]]
        self.signals['velocity'] = np.c_[self.signals['velocity'],
                                         state[self.n_dim:2*self.n_dim, 1:n_samp+1]]
        self.signals['acceleration'] = np.c_[self.signals['acceleration'],
                                             state[2*self.n_dim:3*self.n_dim, 1:n_samp+1]]
        self.signals['time'] = np.r_[
            self.signals['time'], time_axis[1:n_samp+1]] 

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 

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 

def run(self, **parameters):
        self.set_parameters(**parameters)

        self.h0 = self.h_inf(self.V_rest)
        self.m0 = self.m_inf(self.V_rest)
        self.n0 = self.n_inf(self.V_rest)

        initial_conditions = [self.V_rest, self.h0, self.m0, self.n0]

        X = odeint(self.dXdt, initial_conditions, self.time)
        values = X[:, 0]

         # Add info needed by certain spiking features and efel features
        info = {"stimulus_start": self.time[0], "stimulus_end": self.time[-1]}

        return self.time, values, info 

def run(self, **parameters):
        self.set_parameters(**parameters)

        self.h0 = self.h_inf(self.V_rest)
        self.m0 = self.m_inf(self.V_rest)
        self.n0 = self.n_inf(self.V_rest)

        initial_conditions = [self.V_rest, self.h0, self.m0, self.n0]

        X = odeint(self.dXdt, initial_conditions, self.time)
        values = X[:, 0]

         # Add info needed by certain spiking features and efel features
        info = {"stimulus_start": self.time[0], "stimulus_end": self.time[-1]}

        return self.time, values, info 

def step(self, action):
        u = self._bound(action[0], -self._max_u, self._max_u)
        new_state = odeint(self._dynamics, self._state, [0, self._dt],
                           (u,))

        self._state = np.array(new_state[-1])
        self._state[0] = normalize_angle(self._state[0])
        self._state[1] = self._bound(self._state[1], -self._max_omega,
                                     self._max_omega)

        reward = np.cos(self._state[0])

        self._last_u = u.item()

        return self._state, reward, False, {} 

def integrate_numerical(self, t, x0=None):
        if x0 is None:
            x0 = self.x0
        else:
            self.x0 = x0
        sol = odeint(self.ode_func, x0, t)
        return sol